A Perturbation Scheme for Passivity Verification and Enforcement of Parameterized Macromodels

A Perturbation Scheme for Passivity Verification and Enforcement of Parameterized Macromodels

Stefano Grivet-Talocia,  Submitted to the IEEE Transactions on components, packaging and manufacturing technology on 13-Apr-2017S. Grivet-Talocia is with the Department of Electronics and Telecommunications, Politecnico di Torino, Torino 10129, Italy (e-mail: stefano.grivet@polito.it).

This paper presents an algorithm for checking and enforcing passivity of behavioral reduced-order macromodels of LTI systems, whose frequency-domain (scattering) responses depend on external parameters. Such models, which are typically extracted from sampled input-output responses obtained from numerical solution of first-principle physical models, usually expressed as Partial Differential Equations, prove extremely useful in design flows, since they allow optimization, what-if or sensitivity analyses, and design centering.

Starting from an implicit parameterization of both poles and residues of the model, as resulting from well-known model identification schemes based on the Generalized Sanathanan-Koerner iteration, we construct a parameter-dependent Skew-Hamiltonian/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 frequency-parameter 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.

Macromodeling, parameterized modeling, transmission lines, rational approximation, reduced order modeling, Hamiltonian matrices, singular values, perturbation.

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 reduced-order, compact macromodels that can be simulated efficiently either in frequency or time-domain 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, mixed-signal and RF systems including all parasitic effects of electrical interconnects at chip, package and board level, reduced-order modeling of circuit blocks for speeding up simulation tasks, or Computer Aided Design via virtual prototyping, including sensitivity, what-if and design centering processes.

Behavioral reduced-order macromodels are typically derived either through model order reduction of large circuit equations cast in state-space 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 state-space 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 non-rational transfer functions (e.g., in the common case of distributed systems with lossy and dispersive materials). Fitting a (finite-order) rational model to represent such systems cannot lead to exact identification. Such approximation error can be sufficient to make a model non-passive, 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 system-level simulations [14, 15]. Therefore, any passivity violation must be carefully identified and removed.

Passivity verification and enforcement of univariate models is well-established. The common tool that provides an algebraic test for passivity is the so-called Hamiltonian matrix associated to a state-space 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 non-smooth (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 passivity-preserving interpolation scheme is used, also the multivariate model will be uniformly passive [37, 38, 39, 40]. This approach has two main problems. First, passivity-preserving 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 frequency-dependent and parameter-dependent 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 Sanathanan-Koerner 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 well-defined 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 first-principle 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


Since a closed-form 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


The frequency-dependent basis functions are partial fractions based on a prescribed set of distinct real poles and complex pole pairs , where ,




for , where denotes the complex conjugate. The (real-valued) numerator and denominator coefficients are further expressed as


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 a-priori 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 non-smooth 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


in least squares sense. This is achieved here through a linear relaxation of (6) known as (Generalized) Sanathanan-Koerner (GSK) iteration [1, 32, 41]. This approach transforms the nonlinear least-squares 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 above-mentioned 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]:

  1. is regular for ,

  2. ,

  3. 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 parameter-dependent 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 ,


which in turn is equivalent to


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 non-parameterized systems. We start from some initial non-passive model , and we perturb its numerator coefficients as


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 frequency-parameter 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 V-A) and then extends the check to the multivariate case (Section V-B). Section VI presents the main passivity enforcement scheme. In particular, Section VI-A constructs the algebraic constraints that, when iteratively enforced, lead to removal of all passivity violations; Section VI-B constructs a cost function whose minimization will ensure preservation of model accuracy during perturbation; and Section VI-C presents the main iterative passivity enforcement scheme. All these developments require the construction of a parameter-dependent descriptor realization of the model, which is discussed next.

Iv Descriptor Realizations

Starting from model (2), we define the following descriptor realization


where denotes internal generalized states with , and are incident and reflected scattering waves at the model ports. The descriptor matrices are constructed as


where denotes the all-zero matrix block of size . In (11), the matrices and provide regular state-space realizations of model numerator and “extended” denominator , respectively. Individual matrices in these realizations are here defined as


where is the matrix transpose, with


where is the Kronecker product, and


It is straightforward to show that, with the above definitions, the transfer matrix associated to the descriptor form (10)


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.

V-a 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 so-called Hamiltonian matrix associated to the model [23, 16, 17], which in the case of our descriptor form becomes a Skew-Hamiltonian/Hamiltonian (SHH) matrix pencil, also denoted Generalized Hamiltonian Pencil [18, 19, 20, 21]. We define the two block-matrices


where has Hamiltonian structure and is skew-Hamiltonian, 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


for some vector , with . These eigenvalues will be disregarded in the following, and we will consider only the finite eigenvalues of the pencil as


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 parameter-dependent


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 complex-valued 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.

Fig. 1: Graphical illustration of Theorem 1. Left panel: finite SHH eigenvalues (empty dots highlight purely imaginary eigs); right panel: singular value trajectories, which cross the unit threshold at frequencies corresponding to purely imaginary SHH eigenvalues. Local singular value maxima occur at frequencies . For this example, , , and .

As proposed in [23], we define the two sets


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


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 worst-case passivity violation for each non-passive 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.

V-B Uniform Passivity Check of Parameterized Models

The procedure discussed in Section V-A 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


where is the spectral radius of the SHH pencil (16), computed considering only the finite eigenvalues. This function vanishes when the set is non-empty, 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.

Fig. 2: Adaptive parameter sampling, see main text for details.

The proposed check performs an adaptive sampling of within the range . We start with an initial uniform partition into subintervals with endpoints


and we apply the univariate passivity check by computing the SHH eigenvalues at each , as discussed in Section V-A. 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 frequency-parameter plane where the model is uniformly non-passive (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

  • if and , or conversely and , we have a transition from a passive to a non-passive model while sweeping (see Figure 2c). Therefore, we need to refine through (24) in order to track the particular where the onset of a passivity violation occurs.

  • 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 first-order interpolation error


    Refinement is applied if


    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 entire-domain 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 .

1:frequency basis for ;
2:parameter basis for ;
3:model coefficients , in (2) and (5);
4:control parameters , , , , ;
5:set and number of initial samples ;
6:set initial samples as in (23);
8:     for  do
9:         construct SHH pencil ;
10:         find imaginary SHH eigenvalues ;
11:         extract local singular value maxima ;
12:     end for
14:     determine new samples ;
15:until  or
16:return passivity violations .
Algorithm 1 Passivity check of parameterized models

Figure 3 illustrates the proposed passivity check on a practical example (discussed in more detail in Section VII-A). 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 frequency-dependent singular values of model no. 1 for m, where the frequency of the unique imaginary SHH eigenvalue is highlighted by a yellow dot.

Fig. 3: Checking uniform passivity of two PCB link macromodels parameterized by the via antipad radius , see Section VII-A). The two models correspond to two different iterations of the passivity enforcement loop, discussed in Section VI. Panel (a): adaptively computed samples of for the two models. Panels (b) and (c): for all computed , each non-passive frequency band is highligthed with a red line; the yellow dots highlight the frequencies of imaginary SHH eigenvalues; the black squares are the local singular value maxima; the solid black line is an approximation of the contour line for of the singular value trajectories (not required by the algorithm and depicted only for illustration purpose). Panel (d): singular value plot of model at Iteration 1.

Vi Enforcing Passivity of Parameterized Models

Based on the multivariate passivity characterization discussed in Section V-B, we are now ready to formulate our proposed passivity enforcement algorithm. A perturbed model based on (9) is defined as




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




collects the entries of all matrices with a suitable ordering. We have and , with .

Vi-a Building Algebraic Passivity Constraints

Let us consider a single local singular value maximum occurring at frequency , as resulting from the passivity check of Section V-B. Model perturbation leads to a perturbation of this singular value, which under first-order approximation reads [49]


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


Defining now


with the same ordering as in (30), where


we can cast (32) in algebraic form




The constraint (32) operates on an individual singular value maximum , attempting to reduce its value to be less than one (since based on a first-order 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 .

Vi-B 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




denotes the -th entry of the squared model perturbation at frequency and parameter , weighted by the possibly frequency-, parameter-, and entry-dependent weight . Defining now




we can cast the elementwise cost function (38) as




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


where is obtained through an “economy-size” QR factorization of , where . Finally (37) can be cast as




Vi-C Iterative Passivity Enforcement

A first-order singular value perturbation with minimum induced model error is achieved by minimizing (44) while enforcing (35), resulting in the following constrained minimization problem


This problem is convex thus straightforward to solve, e.g., through a Primal-Dual Interior Point method [50]. Since based on a first-order 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 parameter-dependent passivity violations of the model, as discussed in Section V.

1:frequency basis for ;
2:parameter basis for ;
3:model coefficients , in (2) and (5);
4:control parameters , , , , ;
5:find passivity violations via Algorithm 1;
6:while  do
7:     build constraint (35) for each element in ;
8:     compute matrix in (45);
9:     solve convex optimization problem (46);
10:     update model coefficients ;
11:     find passivity violations via Algorithm 1;
12:end while
13:return passive model .
Algorithm 2 Passivity enforcement of parameterized models

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.

Vii-a A Printed Circuit Board Interconnect

The first example we consider is a high-speed 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 frequency-domain scattering responses (Courtesy of Prof. Christian Schuster and Dr. Jan Preibisch, Technische Universität Hamburg-Harburg, Hamburg, Germany) are obtained through a combination of a full-wave field solver (for the connector), lossy transmission-line 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 even-numbered for a-posteriori model validation purposes. The worst-case 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 low-frequency passivity violations throughout the parameter range.

Fig. 4: Samples of for the PCB link macromodel after passivity enforcement.

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 worst-case RMS errors (absolute) and (relative). These errors are only marginally worse than the errors of the original non-passive 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 worst-case passivity violation corresponded to a maximum singular value , which was larger than one by a very small amount). Finally, we depict the resulting implicitly-parameterized macromodel poles in Fig. 6, computed by instantiating the model at the original parameter samples .

Fig. 5: Scattering responses of the PCB link macromodel after passivity enforcement (dashed red lines), compared to the raw data used for model identification (solid blue lines). Both fitting and validation points are displayed; the arrows denote the increasing direction for via antipad radius .

Fig. 6: Parameterized poles (zoomed view) of the PCB link macromodel, computed and superimposed for all original parameter samples .

Vii-B An integrated inductor

We consider here a 2-port 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 full-wave 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 odd-indexed parameter values, leaving the even-numbered responses for model validation purposes, with poles and Chebychev polynomial basis functions for numerator and denominator, respectively.

Fig. 7: Localization of non-passive bands (red lines) through imaginary SHH eigenvalues (yellow dots) of original non-passive inductor model. Top and bottom panels zoom on two different regions of the frequency-parameter plane.

Fig. 8: Normalized SHH spectral distance from imaginary axis for the integrated inductor example, before (Iteration 1) and after (Iteration 6) passivity enforcement.

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.

Fig. 9: Scattering response of the integrated inductor macromodel after passivity enforcement (dashed red lines), compared to the raw data used for model identification (solid blue lines). Both fitting and validation points are displayed; the arrows denote the increasing direction for sidelength .
Error Validation points Fitting points
abs rel abs rel
TABLE I: Absolute and relative RMS errors of non-passive () and passive () inductor models at fitting and validation points.

Vii-C A filter

Fig. 10: Normalized SHH spectral distance from imaginary axis for the double-folded microstrip filter example, before (Iteration 1) and after (Iteration 6) passivity enforcement.

Fig. 11: Localization of non-passive bands (red lines) through imaginary SHH eigenvalues (yellow dots) of filter model at iteration 1 (top) and 2 (bottom) of the passivity enforcement loop.

The last example we propose is a double-folded 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 linearly-spaced stub length values were computed through a field solver at samples spanning the frequency band  GHz. Model identification based on odd-indexed 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).

Fig. 12: Scattering responses of the filter macromodel after passivity enforcement (dashed red lines), compared to the raw data used for model identification (solid blue lines). Only validation points are displayed; the arrows denote the increasing direction for stub length .

Finally, Fig. 12 compares the model responses to the original scattering responses at the validation points (even-indexed parameter samples). Also for this example the accuracy is excellent, with a worst-case 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 Skew-Hamiltonian/Hamiltonian matrix pencil associated to the model, we are able to detect and localize all frequency- and parameter-dependent 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 Hamburg-Harburg, Hamburg, Germany) for sharing the PCB interconnect link data, and to Prof. Piero Triverio (Univ. Toronto, Canada) for sharing the filter data.


  • [1] S. Grivet-Talocia and B. Gustavsen, Passive Macromodeling: Theory and Applications. New York: John Wiley and Sons, 2016.
  • [2] M. Swaminathan, D. Chung, S. Grivet-Talocia, 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 large-scale 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 High-Speed 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. 1052-1061, 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. 2-3, pp. 634–662, 2007.
  • [11] S. Lefteriu and A. C. Antoulas, “A new approach to modeling multi-port 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. Prentice-Hall, 1973.
  • [13] M. R. Wohlers. Lumped and Distributed Passive Networks. Academic press, 1969.
  • [14] P. Triverio, S. Grivet-Talocia, 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. 795-808, 2007.
  • [15] S. Grivet-Talocia, “On driving non-passive 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 skew-Hamiltonian/Hamiltonian matrix pencil transformations.” IEEE Transactions on Circuits and Systems I: Regular Papers vol. 55, no. 2, 2008, pp. 635-643.
  • [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. 61-65, Jan. 2010.
  • [20] Z. Zhang and N. Wong, “An Efficient Projector-Based Passivity Test for Descriptor Systems,” in IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 29, no. 8, pp. 1203-1214, Aug. 2010.
  • [21] Z. Zhang and N. Wong, “Passivity Check of S -Parameter Descriptor Systems via S-Parameter Generalized Hamiltonian Methods,” in IEEE Transactions on Advanced Packaging, vol. 33, no. 4, pp. 1034-1042, 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 Computer-Aided Design of Integrated Circuits and Systems, vol. 31, no. 4, pp. 532-545, April 2012.
  • [23] S. Grivet-Talocia, “Passivity Enforcement via Perturbation of Hamiltonian Matrices” , in IEEE Trans. Circuits and Systems I: Fundamental Theory and Applications, pp. 1755-1769, 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. Grivet-Talocia 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 Frequency-Data,” IEEE Trans. Computed-Aided Design of Integrated Circuits and Systems, vol. 23, no. 2, pp. 293–301, Feb. 2004.
  • [28] G. Calafiore, A. Chinea, and S. Grivet-Talocia, “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. Grivet-Talocia, A. Chinea, G. Calafiore, and L. Daniel, “Efficient localization methods for passivity enforcement of linear dynamical models,” IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 33, pp. 1328Ð1341, Sept 2014.
  • [30] P. Triverio, “Self Consistent, Efficient and Parametric Macromodels for High-speed Interconnects Design”, Ph.D. Thesis, Politecnico di Torino, Italy, April 2009.
  • [31] S. Grivet-Talocia, 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. Grivet-Talocia, and M.S. Nakhla, “A Parameterized Macromodeling Strategy with Uniform Stability Test,” IEEE Transactions on Advanced Packaging, vol. 32, no. 1, pp. 205-215, Feb. 2009.
  • [33] P. Triverio, M.Nakhla and S. Grivet-Talocia, “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. 1-6.
  • [34] P. Triverio, M. Nakhla, and S. Grivet-Talocia, “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. Grivet-Talocia, P. Brenner, “Parameterized and DC-compliant small-signal macromodels of RF circuit blocks,” IEEE Transactions on Components, Packaging, and Manufacturing Technology, vol. 5, pp. 508-Ð522, April 2015.
  • [36] S. Grivet-Talocia and E. Fevola, ”Compact Parameterized Black-Box Modeling via Fourier-Rational Approximations,” in IEEE Transactions on Electromagnetic Compatibility, vol. 59, no. 4, pp. 1133-1142, Aug. 2017.
  • [37] E. R. Samuel, L. Knockaert, F. Ferranti, and T. Dhaene, “Guaranteed Passive Parameterized Macromodeling by Using Sylvester State-Space Realizations”, IEEE Transactions on Microwave Theory and Techniques, vol. 61, no. 4, pp. 1444-1454, 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, “Passivity-Preserving Parametric Macromodeling by Means of Scaled and Shifted State-Space 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 Admittance-Based 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. Grivet-Talocia, 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 Real-Time 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 Vector-Fitting Algorithm,” in IEEE Transactions on Microwave Theory and Techniques, vol. 61, no. 4, pp. 1435-1443, 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. 718-722, Aug. 2016.
  • [46] V. Mehrmann, D. Watkins, “Structure-preserving methods for computing eigenpairs of large sparse skew-Hamiltonian/Hamiltonian pencils.” SIAM Journal on Scientific Computing, vol. 22, no. 6, 2001, pp. 1905-1925.
  • [47] D.S. Watkins, “On Hamiltonian and symplectic Lanczos processes.” Linear algebra and its applications, vol. 385, 2004, pp. 23-45.
  • [48] P. Benner, R. Byers, V. Mehrmann, H. Xu, “Numerical computation of deflating subspaces of skew-Hamiltonian/Hamiltonian pencils,” SIAM Journal on Matrix Analysis and Applications, vol. 24, no. 1, pp. 165-190, 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 Variability-Aware Analysis Method for High-Speed Digital Link Design Using PCE”, DesignCon, Jan 31-Feb 2, 2017, Santa Clara (CA), USA.
  • [52] X. Duan, R. Rimolo-Donadio, H.-D. Brüns, and C. Schuster, “Circular ports in parallel-plate waveguide analysis with isotropic excitations,” IEEE Transactions on Electromagnetic Compatibility, vol. 54, pp. 603Ð612, June 2012.
  • [53] P. Triverio, M. Nakhla, and S. Grivet-Talocia, Ò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 27-28 2010.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description