1

A Perturbation Scheme for Passivity Verification
and Enforcement of Parameterized Macromodels

arXiv:1706.06395v1 [] 20 Jun 2017

Stefano Grivet-Talocia, Senior Member, IEEE

Abstract—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 SanathananKoerner iteration, we construct a parameter-dependent 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 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.
Index Terms—Macromodeling, parameterized modeling, transmission lines, rational approximation, reduced order modeling,
Hamiltonian matrices, singular values, perturbation.

I. I NTRODUCTION

AND

M OTIVATION

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],
[2], [5], [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],
Submitted to the IEEE Transactions on components, packaging and manufacturing technology on 13-Apr-2017
S. Grivet-Talocia is with the Department of Electronics and Telecommunications, Politecnico di Torino, Torino 10129, Italy (e-mail: stefano.grivet@polito.it).

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]–[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 [16],
[23]. 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 [22], [23], singular value/eigenvalue
perturbation [24], [25], Positive or Bounded Real Lemma [17]
enforcement based on semidefinite programming [27], H∞
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],

2

[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]–[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 n̄, the multivariate model
will have a larger order Ln̄, where L 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 n̄ 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 parameter-dependent basis functions [32]–[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. P RELIMINARIES

AND

N OTATION

Let us consider a generic LTI multiport structure with
transfer matrix H̆(s; ϑ), where s is the Laplace variable, and
ϑ ∈ Θ = [ϑmin , ϑmax ] is some external parameter, e.g.,
some material or geometrical characteristic of the underlying
physical system. Throughout this work, it is assumed that

H̆(s; ϑ) ∈ CP ×P 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
first-principle model (e.g., the numerical solution of Maxwell’s
equations for an electromagnetic system) at a set of fixed
frequencies fk , k = 1, . . . , k̄ over a given frequency band
[fmin , fmax ], with f1 = fmin and fk̄ = fmax , and for a set of
fixed parameter values ϑm , m = 1, . . . , m̄ spanning the range
[ϑmin , ϑmax ]. We will denote this characterization as
H̆k,m = H̆(j2πfk ; ϑm ),

k = 1, . . . , k̄, m = 1, . . . , m̄. (1)

Since a closed-form dependence of the true system response
on s 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
Pn̄
Rn (ϑ) ϕn (s)
N(s, ϑ)
= Pn=0
.
(2)
H(s; ϑ) =
n̄
D(s, ϑ)
n=0 rn (ϑ) ϕn (s)

The frequency-dependent basis functions ϕn (s) are partial
fractions based on a prescribed set of distinct n̄r real poles
qn ∈ R− and n̄c complex pole pairs pn = p′n ± jp′′n ∈ C− ,
where ϕ0 (s) = 1,
ϕn (s) = (s − qn )−1 ,
and

for n = 1, . . . , n̄r ,

ϕn̄r +2ν−1 (s) = (s − pν )−1 + (s − p∗ν )−1
ϕn̄r +2ν (s) = j(s − pν )−1 − j(s − p∗ν )−1

(3)

(4)

for ν = 1, . . . , n̄c , where ∗ denotes the complex conjugate.
The (real-valued) numerator and denominator coefficients are
further expressed as
Rn (ϑ) =

ℓ̄
X

Rn,ℓ ξℓ (ϑ),

ℓ=1

rn (ϑ) =

ℓ̄
X

rn,ℓ ξℓ (ϑ)

(5)

ℓ=1

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], [36], [37], [39], [40],
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 H(s; ϑ) is a rational
function of s, 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 N(s, ϑ) and denominator D(s, ϑ) on the
respective coefficients Rn,ℓ and rn,ℓ , the latter can be easily
evaluated during a model identification step by enforcing the
fitting condition
H(j2πfk ; ϑm ) ≈ H̆k,m ,

k = 1, . . . , k̄, m = 1, . . . , m̄. (6)

3

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 r0,1 = 1.
III. P ROBLEM 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 H(s; ϑ) is Bounded
Real [12], [13]:
1) H(s; ϑ) is regular for Re {s} > 0,
2) H∗ (s; ϑ) = H(s∗ ; ϑ),
3) IP − HH (s; ϑ)H(s; ϑ) ≥ 0 for Re {s} > 0,
where H is the Hermitian transpose, and IP is the identity
matrix of size P . 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
s = jω,
IP − HH (jω; ϑ)H(jω; ϑ) ≥ 0

∀ω ∈ R,

(7)

which in turn is equivalent to
σmax {H(jω; ϑ)} ≤ 1,

∀ω ∈ R,

∀n, ℓ

IV. D ESCRIPTOR R EALIZATIONS
Starting from model (2), we define the following descriptor
realization
(
Eẋ = A(ϑ)x + Bu
(10)
y = C(ϑ)x
where x ∈ RN +P denotes internal generalized states with
N = n̄P , and u, y ∈ RP are incident and reflected scattering waves at the model ports. The descriptor matrices are
constructed as




IN
0N,P
A0
B0
E=
,
A(ϑ) =
,
0P,N 0P,P
C2 (ϑ) D2 (ϑ)

(8)

where σmax {·} 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
H(s; ϑ) throughout the parameter range, so that the above
Boudned Realness conditions, in particular (8), will hold
∀ϑ ∈ [ϑmin , ϑmax ].
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 H(s; ϑ), and we perturb its
numerator coefficients as
b n,ℓ = Rn,ℓ + ∆Rn,ℓ
R

by determining the correction ∆Rn,ℓ that is required to ensure
that the perturbed scattering response of the parameterized
b ϑ) is uniformly Bounded Real in the prescribed
model H(s;
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 D(s, ϑ). 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 parameterdependent descriptor realization of the model, which is discussed next.

(9)



C(ϑ) = C1 (ϑ)



D1 (ϑ) ,




0N,P
B(ϑ) =
,
−IP (ϑ)

(11)
where 0J,K denotes the all-zero matrix block of size
J × K. In (11), the matrices {A0 , B0 , C1 (ϑ), D1 (ϑ)} and
{A0 , B0 , C2 (ϑ), D2 (ϑ)} provide regular state-space realizations of model numerator N(s, ϑ) and “extended” denominator
D(s, ϑ)IP , respectively. Individual matrices in these realizations are here defined as

where

T

A0 = blkdiag{A0r , A0c }
 T

T
BT
0 = B0r , B0c

(12)

is the matrix transpose, with

r
A0r = blkdiag{qn IP }n̄n=1
 ′
n̄c
pn IP
p′′n IP
A0c = blkdiag
−p′′n IP p′n IP n=1

T
B0r = 1, . . . , 1 ⊗ IP

T
B0c = 2, 0, . . . , 2, 0 ⊗ IP

(13)

4

where ⊗ is the Kronecker product, and


C1 (ϑ) = R1 (ϑ) · · · Rn̄ (ϑ)


C2 (ϑ) = IP r1 (ϑ) · · · IP rn̄ (ϑ)
D1 (ϑ) = R0 (ϑ)
D2 (ϑ) = IP r0 (ϑ)

Im {λ} = ω
σ

(14)

H(s; ϑ) = C(ϑ)(sE − A(ϑ))−1 B

Re {λ}

(15)

matches (2). We leave this tedious verification to the Reader,
see also [32].
OF

PARAMETERIZED M ODELS

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.
A. Passivity Check of Univariate Models

M(ϑ)v 0 = jω0 Kv 0
The proof is straightforward and follows the same flow as the
proof of Theorem 1 in [16]. See also [1], [19], [21].
Matrix K is singular, since E is singular, with
dim ker(E) = P and dim ker(K) = 2P . Therefore, we expect
that the eigenspectrum of the pencil (M(ϑ), K) includes at
least 2P infinite eigenvalues. We recall that infinite eigenvalues are characterized as pairs (α, β) such that
αM(ϑ)v = βKv

(17)

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

−ω2

0

ω1 ω̄1 ω2

(18)

ω

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 ωi corresponding
to purely imaginary SHH eigenvalues. Local singular value maxima σ̄i occur
at frequencies ω̄i . For this example, P = 2, ν = 2, Ip = {0, 2} and
Inp = {1}.

that satisfy (17) for some v 6= 0. Of course, all these
eigenvalues depend on the considered parameter value ϑ.
Therefore, the (finite part of the) SHH eigenspectrum is
parameter-dependent
Λ(ϑ) = {λ(ϑ) = β(ϑ)/α(ϑ) : α(ϑ) 6= 0}.

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




E 0
A(ϑ)
BBT
M(ϑ) =
,
K
=
−CT (ϑ)C(ϑ) −AT (ϑ)
0 ET
(16)
where M(ϑ) has Hamiltonian structure and K is skewHamiltonian, forming the SHH matrix pencil (M(ϑ), K). We
have the following result:
Theorem 1. Assume that the pencil (A(ϑ), E) has no purely
imaginary eigenvalues. Then, the following conditions are
equivalent:
• σ = 1 is a singular value of H(jω0 ; ϑ) with |ω0 | < +∞;
• jω0 is an eigenvalue of the pencil (M(ϑ), K), i.e., there
exist some vector v 0 6= 0 such that

λ = β/α,

σ=1

ω1
−ω1

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

V. C HECKING PASSIVITY

σ̄1

ω2

(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 complex-valued with
a nonvanishing real and imaginary part [46]–[48].
Based on Theorem 1, if Λ(ϑ) includes a purely imaginary
eigenvalue jω0 , then the transfer matrix of the model has a
singular value σ(jω0 , ϑ) = 1. We conclude that the singular
value trajectory for fixed ϑ in a neighborhood of jω0 crosses
the unit threshold, which is the critical condition leading to a
passivity violation, see (8). This violation will be for ω > ω0
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
ν(ϑ)

χ(ϑ) = {ωi (ϑ)}i=1 ,

χ̄(ϑ) = {0} ∪ χ(ϑ) ∪ {+∞} (20)

collecting the frequencies ωi (ϑ) > 0 corresponding to the the
ν(ϑ) purely imaginary SHH eigenvalues λi (ϑ) = jωi (ϑ) with
positive (finite) imaginary part, sorted in ascending order. The
augmented set χ̄(ϑ) includes also the DC point ω0 = 0 and
the infinte frequency ω = +∞. The elements of χ̄(ϑ) thus
induce a subdivision of the positive frequency axis [0, +∞)
into ν(ϑ) + 1 disjoint subbands
Ωi (ϑ) = (ωi (ϑ), ωi+1 (ϑ)),

i = 0, . . . , ν(ϑ),

(21)

where we defined ω0 (ϑ) = 0 and ων(ϑ)+1 (ϑ) = +∞. Since
all intersections of some singular value trajectory with the
threshold σ = 1 is captured in the set χ(ϑ), each subband
Ωi (ϑ) 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
i ∈ Ip (ϑ) and i ∈ Inp (ϑ), respectively. The worst-case

5

passivity violation for each non-passive band is defined as the
maximum singular value σ̄i (ϑ), occurring at some frequency
ω̄i (ϑ) ∈ Ωi (ϑ) with i ∈ Inp (ϑ). Numerical estimates of such
local maxima (ω̄i (ϑ), σ̄i (ϑ)) are easily determined through
local sampling. See Figure 1 for a graphical illustration.
B. Uniform Passivity Check of Parameterized Models

σ̄i (ϑ)

ω̄i (θµ+1 )

ωi (θµ )

(22)

where ρ(ϑ) = maxλ(ϑ)∈Λ(ϑ) |λ(ϑ)| 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 ϑ ∈ Θ = [ϑmin , ϑmax ]. We start with an initial
uniform partition into µ̄ subintervals with endpoints
µ
(23)
θµ = ϑmin + (ϑmax − ϑmin ), µ = 0, . . . , µ̄,
µ̄
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
Θµ = [θµ , θµ+1 ],
• if ψ(ϑ) = 0, ∀ϑ ∈ Θµ , 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 Ωi (ϑ) where at least one singular
value of the model H(jω; ϑ) exceeds one, with a local
maximum σ̄i (ϑ) occurring at frequency ω̄i (ϑ).
• if ψ(ϑ) > 0, ∀ϑ ∈ Θµ , the model is uniformly passive
in Θµ .
Assuming now that only information at the edges {θµ , θµ+1 }
of Θµ is available, we have the following subcases (Fig. 2):
• if ψ(θµ ) = ψ(θµ+1 ) = 0 and ν(θµ ) = ν(θµ+1 ),
we can infer that there exist at least one region Γµ,i
of the frequency-parameter plane Γ ⊆ [0, +∞) × Θµ
where the model is uniformly non-passive (see Figure 2a).
For each region Γµ,i , we know the local singular value
maxima σ̄i (θµ ), σ̄i (θµ+1 ), and the corresponding localization frequencies ω̄i (θµ ), ω̄i (θµ+1 ). In addition, we

ωi (θµ+1 )

σmax (jω, ϑ) > 1

ϑ
ψ(θµ+1 ) = 0
ν(θµ+1 ) = 2
ψ(θµ ) = 0
ν(θµ ) = 2

θµ+1

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
| Re {λ(ϑ)} |
ψ(ϑ) = min
ρ(ϑ)
λ(ϑ)∈Λ(ϑ)

ω̄i (θµ )

Γµ,i
θµ

ω

(a)
ϑ

ψ(θµ+1 ) = 0
ν(θµ+1 ) = 2
ψ(θµ ) = 0
ν(θµ ) = 4

θµ+1
Γµ,i
θµ

ω

(b)
ϑ

ψ(θµ+1 ) > 0
ν(θµ+1 ) = 0
ψ(θµ ) = 0
ν(θµ ) = 2

θµ+1
θµ

ω

(c)
ϑ

ψ(θµ+1 ) > 0
ν(θµ+1 ) = 0
ψ(θµ ) > 0
ν(θµ ) = 0

θµ+1
θµ

ω

(d)
ϑ
θµ+1

ψ(θµ+1 ) > 0
ν(θµ+1 ) = 0

θµ

ψ(θµ ) > 0
ν(θµ ) = 0

(e)
Fig. 2. Adaptive parameter sampling, see main text for details.

ω

6

•

•

•

know all the edges of Γµ,i for θµ and θµ+1 from the
sets χ(θµ ), χ(θµ+1 ). Since local passivity violations have
been identified, there is no need to refine subinterval Θµ .
if ψ(θµ ) = ψ(θµ+1 ) = 0 and ν(θµ ) 6= ν(θµ+1 ),
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
θµ+1/2 = 1/2(θµ + θµ+1 ).
(24)
if ψ(θµ ) = 0 and ψ(θµ+1 ) > 0, or conversely ψ(θµ ) > 0
and ψ(θµ+1 ) = 0, 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 ψ(θµ ) > 0 and ψ(θµ+1 ) > 0, we can have two
subcases:
– ψ(ϑ) > 0 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
ψ(θµ+1/2 ), and by constructing the first-order interpolation error
εµ+1/2 = ψ(θµ+1/2 ) − 1/2[ψ(θµ ) + ψ(θµ+1 )] .

(25)

Refinement is applied if
εµ+1/2 > γ ψ(θµ+1/2 ) ,

(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 γ = 0.2.
We perform a total number M of refinement passes. At each
pass, the new points θµ+1/2 are added to the previous subset
{θµ , µ = 0, . . . , µ̄}, 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
θµm at each refinement pass m, according to the above adaptive
process). Throughout this work, we use M = 10.
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 κ > 1. For all documented examples, we set κ = 4.

Algorithm 1 Passivity check of parameterized models
Require: frequency basis ϕn for n = 0, . . . , n̄;
Require: parameter basis ξℓ for ℓ = 1, . . . , ℓ̄;
Require: model coefficients Rn,ℓ , rn,ℓ in (2) and (5);
Require: control parameters ϑmin , ϑmax , γ, κ, M ;
1: set m = 0 and number of initial samples µ̄0 = κ ℓ̄;
0
2: set initial samples S0 = {θµ
, µ = 0, . . . , µ̄0 } as in (23);
3: repeat
4:
for µ = 1, . . . , µ̄m do
5:
construct SHH pencil (M(θµm ), K);
6:
find imaginary SHH eigenvalues ωi (θµm );
m
m
7:
extract local singular value maxima (ω̄µ,i
, σ̄µ,i
);
8:
end for
9:
m← m+1
10:
determine new samples Sm = {θµm , µ = 1, . . . , µ̄m };
11: until Sm = ∅ or m = M
m
m
12: return passivity violations V = ∪m {(ω̄µ,i
, σ̄µ,i
), ∀µ, i}.

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 ψ(θµ ) = 0 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 ϑ ∈ Θ∗ ≈ [437, 454] µ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 ψ(ϑ) = 0 and ψ(ϑ) > 0, which define
the endpoints of Θ∗ . Finally, Figure 3(d) shows the frequencydependent singular values of model no. 1 for ϑ = 400 µm,
where the frequency of the unique imaginary SHH eigenvalue
is highlighted by a yellow dot.
VI. E NFORCING PASSIVITY

OF

PARAMETERIZED M ODELS

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

with
∆H(s; ϑ) =

b ϑ) = H(s; ϑ) + ∆H(s; ϑ)
H(s;
∆N(s, ϑ)
=
D(s, ϑ)

Pn̄

(27)

Pℓ̄
n=0
ℓ=1 ∆Rn,ℓ ξℓ (ϑ) ϕn (s)
Pn̄ Pℓ̄
n=0
ℓ=1 rn,ℓ ξℓ (ϑ) ϕn (s)

(28)
The coefficient perturbations ∆Rn,ℓ will be computed such
that the singular values of the perturbed model that are larger

7

3

×10

A. Building Algebraic Passivity Constraints

(a) Normalized SHH spectral distance from imaginary axis

-5

Iteration 1
Iteration 2

2

1

0
400

420

440

460

480

500

520

540

560

580

600

ϑ (µm)

where uµ,i , v µ,i are the left and right singular vectors of
H(jω̄µ,i ; θµ ) associated to σ̄µ,i . Forcing this singular value
to be less than one leads to the inequality constraint

Re uH
(32)
µ,i ∆H(jω̄µ,i ; θµ ) v µ,i ≤ 1 − σ̄µ,i .

(b) Localization of non-passive bands, Iteration 1
600

ϑ (µm)

550
500
450
400

Let us consider a single local singular value maximum
σ̄µ,i = σ̄i (θµ ) > 1 occurring at frequency ω̄µ,i = ω̄i (θµ ),
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]

(31)
σ
bµ,i ≈ σ̄µ,i + Re uH
µ,i ∆H(jω̄µ,i ; θµ ) v µ,i ,

Defining now

0

0.5

1

1.5

2


T
aµ,i = aµ,i;0,1 , · · · , aµ,i;n,ℓ , · · · , aµ,i;n̄,ℓ̄

2.5

Frequency (Hz)

×10

8

(33)

with the same ordering as in (30), where

(c) Localization of non-passive bands, Iteration 2

ϑ (µm)

460

aµ,i;n,ℓ =

450

ξℓ (θµ ) ϕn (jω̄µ,i )
,
D(jω̄µ,i , θµ )

(34)

we can cast (32) in algebraic form

440

pT
µ,i x ≤ 1 − σ̄µ,i ,

430
0

0.5

1

1.5

2

2.5

Frequency (Hz)

×10

4

where

(d) Singular values at Iteration 1, ϑ=400 µm
1

0.998

0.996
0

1

2

3

4

5

6

Frequency (Hz)

7

8

(35)

3

9

10
×10 8

Fig. 3. Checking uniform passivity of two PCB link macromodels parameterized by the via antipad radius ϑ = r, 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 σ = 1 of the singular value trajectories σ(jω; ϑ) (not
required by the algorithm and depicted only for illustration purpose). Panel
(d): singular value plot of model at Iteration 1.

 T
H
T
pT
µ,i = Re (v µ,i ⊗ uµ,i ) ⊗ aµ,i .

(36)

The constraint (32) operates on an individual singular value
maximum σ̄µ,i , 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 µ, i.
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
E2 =

P
X

2
Ei,j
,

(37)

i,j=1

where
than one will be displaced below the unit threshold. Such
coefficients will be collected in a global vector of decision
varibles

T
T
T
x = xT
,
(29)
1,1 , · · · , xi,j , · · · , xP,P
where


T
xi,j = (∆R0,1 )i,j , · · · , (∆Rn,ℓ )i,j , · · · , (∆Rn̄,ℓ̄ )i,j
(30)
collects the (i, j) entries of all matrices ∆Rn,ℓ with a suitable
ordering. We have xi,j ∈ R(n̄+1)ℓ̄ and x ∈ RQ , with Q =
P 2 (n̄ + 1)ℓ̄.

2
Ei,j
=

k̄ X
m̄
X

2
wi,j;k,m
|∆Hi,j (j2πfk , ϑm )|2

(38)

k=1 m=1

denotes the (i, j)-th entry of the squared model perturbation
at frequency fk and parameter ϑm , weighted by the possibly
frequency-, parameter-, and entry-dependent weight wi,j;k,m .
Defining now
bk,m;n,ℓ =

ξℓ (ϑm ) ϕn (j2πfk )
D(j2πfk , ϑm )

(39)

and


bT
k,m = bk,m;0,1 , · · · , bk,m;n,ℓ , · · · , bk,m;n̄,ℓ̄

(40)

8

we can cast the elementwise cost function (38) as
2

2
Ei,j
= kFi,j xi,j k2 ,

(41)

where

Fi,j



o
n
e i,j
Re F
o ,
= n
e i,j
Im F



e i,j
F


wi,j;1,1 bT
1,1


..


.



T 
= wi,j;k,m bk,m 


..




.
wi,j;k̄,m̄ bT
k̄,m̄

(42)

Note that Fi,j ∈ R2k̄m̄×(n̄+1)ℓ̄ collects as many rows as
available frequency and parameter samples. The row size can
thus be large, since usually 2k̄ m̄ ≫ (n̄ + 1)ℓ̄. An equivalent
compressed form of (41) reads
2

2
= kΨi,j xi,j k2 ,
Ei,j

(43)

where Ψi,j ∈ R(n̄+1)ℓ̄×(n̄+1)ℓ̄ is obtained through an
“economy-size” QR factorization of Fi,j = Qi,j Ψi,j , where
QT
i,j Qi,j = I. Finally (37) can be cast as
2

E 2 = kΨ xk2 ,

(44)

Ψ = blkdiag{Ψi,j }P
i,j=1

(45)

where

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
2

minx kΨ xk2

subject to pT
µ,i x ≤ 1 − σ̄µ,i ,

∀µ, i. (46)

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.
VII. E XAMPLES
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.

Algorithm 2 Passivity enforcement of parameterized models
Require: frequency basis ϕn for n = 0, . . . , n̄;
Require: parameter basis ξℓ for ℓ = 1, . . . , ℓ̄;
Require: model coefficients Rn,ℓ , rn,ℓ in (2) and (5);
Require: control parameters ϑmin , ϑmax , γ, κ, M ;
1: find passivity violations V via Algorithm 1;
2: while V 6= ∅ do
3:
build constraint (35) for each element in V;
4:
compute matrix Ψ in (45);
5:
solve convex optimization problem (46);
6:
update model coefficients Rn,ℓ ← Rn,ℓ + ∆Rn,ℓ ;
7:
find passivity violations V via Algorithm 1;
8: end while
b ϑ).
9: return passive model H(s;
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 2 × 2 scattering
responses of the link are parameterized by the via antipad
radius ϑ = r within the range Θ = [400, 600] µm. The
raw frequency-domain scattering responses (Courtesy of Prof.
Christian Schuster and Dr. Jan Preibisch, Technische Universitä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 k̄ = 500 frequency samples up to 10 GHz, combined
with m̄ = 9 parameter samples were available for model
identification.
The parameterized GSK iteration of Section II was applied
with n̄ = 44 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 ℓ̄ = 3, corresponding to quadratic
polynomials. Only one half of the parameter samples (indices
m = 1, 3, 5, 7, 9) were used to fit the model, leaving the evennumbered m = 2, 4, 6, 8 for a-posteriori model validation purposes. The worst-case absolute RMS error among all scattering
matrix elements at both fitting and validation points resulted
8.788 × 10−4 , whereas the relative error was 3.367 × 10−3 .
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.
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 ψ(ϑ) > 0
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 8.792 × 10−4
(absolute) and 3.368 × 10−3 (relative). These errors are only
marginally worse than the errors of the original non-passive

9

×10

Localization of non-passive bands

Normalized SHH spectral distance from imaginary axis

-4

1.52

Sidelength (mm)

2.5

2

1.5

1.515
1.51
1.505
1.5
0

1
400

420

440

460

480

500

520

540

560

580

2

4

8

10

12

14

16

Frequency (Hz)

18
×10

9

Localization of non-passive bands

Scattering response S 11, magnitude
S(1,1), data
S(1,1), model

1.05

Sidelength (mm)

Fig. 4. Samples of ψ(ϑ) for the PCB link macromodel after passivity
enforcement.

0.6

6

600

ϑ (µm)

1.04

1.03

1.02
0

ϑ

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Frequency (Hz)

2
×10

10

0.4

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.

0.2

0

0

1

2

3

4

5

6

7

8

9

1

10
×10 9

Frequency (Hz)
Scattering response S , magnitude
12

Normalized SHH spectral distance from imaginary axis

0.04
0.03
0.02

0.9

0.01

0.8

Iteration 1
Iteration 6

ϑ
S(1,2), data
S(1,2), model

0

1

1.1

1.2

1.3

1.4

1.5

Sidelength (mm)

0.7
0

1

2

3

4

5

6

7

8

9

Frequency (Hz)

10
×10

9

Fig. 5. Scattering responses Sij (jω; ϑ) 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. 8. Normalized SHH spectral distance from imaginary axis for the
integrated inductor example, before (Iteration 1) and after (Iteration 6)
passivity enforcement.

B. An integrated inductor

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 σmax = 1.000346, 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 ϑm .

Real part (Hz)

-3

Parameterized model poles

×10 8

-3.5
-4
-4.5
-5
-5.5
0

1

2

3

4

5

Imag part (Hz)

6

7

8

9

10
×10 9

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

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
m̄ = 11 different values of the sidelength L = ϑ ranging
from 1.02 to 1.52 mm, with k̄ = 477 frequency samples up to
12 GHz. Model identification via GSK iteration was performed
using only odd-indexed parameter values, leaving the evennumbered responses for model validation purposes, with n̄ = 8
poles and l̄ = 4, 3 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.

10

S12, magnitude

1

Localization of non-passive bands, Iteration 1
Stub length (mm)

ϑ
0.9
0.8
0.7
0.6

S(1,2), data
S(1,2), model

2.2
2.15
2.1

0

2

4

6

8

10

12

Frequency (Hz)

×10

0

9

0.5

Stub length (mm)

ϑ
0
-100
0

2

4

6

8

10

12

Frequency (Hz)

2.5

3
×10 10

×10

2.2
2.15
2.1
0

Error

-40

points
rel
1.90 × 10−3
1.95 × 10−3
2.24 × 10−3
2.12 × 10−3
1.93 × 10−3
2.72 × 10−3

1

1.5

2

2.5

3
×10 10

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.

-20

Fitting
abs
7.09 × 10−4
1.74 × 10−3
8.91 × 10−4
8.58 × 10−4
1.72 × 10−3
1.10 × 10−3

0.5

Frequency (Hz)

TABLE I
A BSOLUTE AND RELATIVE RMS ERRORS OF NON - PASSIVE (H) AND
b INDUCTOR MODELS AT FITTING AND VALIDATION POINTS .
PASSIVE ( H)

H11
H12
H22
b 11
H
b 12
H
b 22
H

2

2.25

9

Fig. 9. Scattering response S12 (jω; ϑ) 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 ϑ.

Validation points
abs
rel
1.84 × 10−3
4.64 × 10−3
2.60 × 10−3
2.96 × 10−3
2.04 × 10−3
5.13 × 10−3
1.80 × 10−3
4.54 × 10−3
2.56 × 10−3
2.92 × 10−3
1.94 × 10−3
4.87 × 10−3

1.5

Localization of non-passive bands, Iteration 2

S(1,2), data
S(1,2), model

100

1

Frequency (Hz)

S12, phase (degrees)

200

-200

2.25

0

Scattering responses, magnitude (dB)

ϑ

S(2,1), data
S(2,1), model

-60
0.5

1

1.5

200

2
×10 10

Frequency (Hz)
Scattering responses, phase (degrees)

S(2,1), data
S(2,1), model

100

ϑ
0
-100
-200
0.5

C. A filter

1

1.5

Frequency (Hz)

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 ϑ = L of a microstrip stub
within the range 2.08–2.28 mm. Scattering responses for m̄ =
21 linearly-spaced stub length values were computed through
a field solver at k̄ = 300 samples spanning the frequency
band [5, 20] GHz. Model identification based on odd-indexed
parameter samples required n̄ = 10 poles and l̄ = 3 (quadratic)
polynomial basis functions to capture parameter variations.
The resulting model was not passive, as Fig. 10 shows
by depicting ψ(ϑ) = 0 at the first passivity iteration loop.

8

Normalized SHH spectral distance from imaginary axis

×10 -3

Iteration 1
Iteration 6

6
4
2
0
2.08

2.1

2.12

2.14

2.16

2.18

2.2

2.22

2.24

2.26

2.28

Stub length (mm)

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.

2
×10 10

Fig. 12. Scattering responses S21 (jω; ϑ) 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 ϑ.

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 ψ(ϑ) > 0 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 ϑ ∈ [2.08, 2.13] 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 (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
3.16 × 10−3 (absolute) and 4.64 × 10−3 (relative).

11

VIII. C ONCLUSION
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 Universitä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.
R EFERENCES
[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. 288310, 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.
863886, Oct 2009.
[16] S. Boyd, V. Balakrishnan, P. Kabamba, “A bisection method for computing the H∞ 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. 29903003, 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. 13281341, 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.

12

[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. 920933, 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 VectorFitting 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. 19051925.
[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. 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 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. Brüns, and C. Schuster, “Circular
ports in parallel-plate waveguide analysis with isotropic excitations,”
IEEE Transactions on Electromagnetic Compatibility, vol. 54, pp.
603612, 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. 393396, September 27-28 2010.

