Parametric model order reduction for
large-scale and complex thermal systems

arXiv:1803.05240v1 [cs.SY] 14 Mar 2018

Daming Lou and Siep Weiland1
also called the one-sided moment-matching technique [12]
was proposed to match r moments of a transfer function by
a rth-order approximation. Meanwhile, two-sided moment
matching, known as Lanczos method [6], was introduced to
perform matching of 2r moments by a rth-order approximations. Then the block Arnoldi [7] and block Lanczos method
[4] were proposed to solve the same question of multi-input
multi-output (MIMO) systems. However, for large-scale systems which are parameter dependent, the aforementioned
methods are not suitable anymore for obtaining a low-order
model while maintaining the physical interpretation after the
MOR procedure. Therefore, there is a need to develop such
MOR techniques that allows the reduction of parametervarying systems. Hence, a new branch of MOR is known
as parametric model order reduction (pMOR).
The first work of pMOR using Krylov subspace methods
was introduced in [13], where the state evaluation matrix
linearly depends on a single parameter. In consideration of
the physical systems that can be described with multiple
parameters, an extensive pMOR method which matches the
coefficients of multivariate Taylor series is generalized in [5].
This method replaces the product of complex independent
variables s and the parameter p in the transfer function
G(s, p) by a set of redefined parameters, where the moments
of the reduced model are equivalent to the corresponding
moments of the original model with respect to the new
expansion points. Contrary to singular value decomposition
(SVD) based truncation, the pMOR techniques based on
Krylov subspaces have no prior error bound or global error
bound.
In this paper we explore MOR techniques for parametric
models with the purpose to match multiple moments of
the transfer function. We introduce a separation between
the physical parameters p and the frequency points s such
that, unlike previous work, the physical interpretation of the
physical parameters in the reduced model is maintained.
Also, simulation results demonstrate the benefits of this
separation. Moreover, an upper bound of the approximation
error between the reduced model and the original model
is derived, which allows prior estimating the misfit without
calculating the exact reduced model.
The remainder of the paper is structured as follows: In
Section II, a condensed problem description for parametric
model order reduction is presented. The proposed method for
parametric model order reduction with moment matching of
multiple frequencies and multiple parameters is introduced
in Section III. The analytic expression of the upper bound
of this method is given in section IV. Simulation results are

Abstract— In this paper, a parametric model order reduction
(pMOR) technique is proposed to find a simplified system
representation of a large-scale and complex thermal system. The
main principle behind this technique is that any change of the
physical parameters in the high-fidelity model can be updated
directly in the simplified model. For deriving the parametric
reduced model, a Krylov subspace method is employed which
yields the relevant subspaces of the projected state. With the
help of the projection operator, first moments of the low-rank
model are set identical to the correspondent moments of the
original model. Additionally, a prior upper bound of the error
induced by the approximation is derived.

I. I NTRODUCTION
Numerical simulations of complex dynamical systems are
an indispensable tool in studying thermodynamic phenomena. However, for complex thermal systems where ultra-high
precision simulations are required, the finite element method
(FEM) commonly yields large-scale models. These models
demand considerable computational resources. Therefore,
model order reduction techniques are employed to reduce
the computational complexity by replacing the high-order
dynamic model with a low-order one. For the class of linear
time invariant (LTI) systems, many reduction techniques have
reached a relatively high level of maturity [2]. For systems
with uncertain or time-varying parameters there is a persistent need of novel reduction techniques. Since creating a new
reduced model for every parameter value is inefficient and
computationally costly, there is a strong need for parametric
model order reduction (pMOR) techniques as introduced by
[3]. We distinguish among reduction techniques in which
both the state dimension and parameter vector dimension are
reduced, and techniques in which state dimension is reduced
while the parameter vector keeps its physical relevance.
In addition, we distinguish among time-varying and timeinvariant parametric dependence of high-fidelity models. In
either of these cases, the objective is to find a low-cost, but
accurate, parameterized model.
In the last decades, model order reduction (MOR) techniques based on Padé approximation [1] have been used as
powerful tools for large-scale system simulation, in particular
for very-large-scale system integration (VLSI) simulation [8]
and microelectromechanical (MEMS) simulation tasks [9].
The first algorithm was an asymptotic waveform evaluation
in [10], which uses explicit moment-matching. For a singleinput single-output (SISO) system, the Arnoldi algorithm,
1 Daming Lou, Siep Weiland are with the Department of Electrical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands.
Emails: {D.Lou, S.Weiland}@tue.nl

1

shown in Section V. Finally, Section VI concludes this work.

and the state space representation of the reduced model is
then given by

Er (p)ẋr (t) = Ar (p)xr (t) + Br (p)u(t)
Σr (p) :=
(6)
y(t) = Cr (p)xr (t)

II. P ROBLEM D ESCRIPTION
The main goal of the reduction techniques proposed
in this paper is to preserve parameters in the system as
symbolic quantities in the reduced order model. In fact,
there are many key applications where this parametric model
reduction technique is of crucial interest. These include
design optimizations (where the parameter vector represents
design parameters) in which the design-loop can substantially
benefit from simplified parametric models, calibration of
parameters through simplified models or control design for
parameter dependent models. In either of these cases, a
change in parameters does not require to repeat the reduction
procedure, but simply the evaluation of the parametric model.
The aim is to infer a reduced order model with smaller state
dimension, but which remains explicit in the parameter p.
Consider a linear parameter-varying model with parameter
p(t) = [p1 (t), · · · , p` (t)]T ∈ R` described as

E(p)ẋ(t) = A(p)x(t) + B(p)u(t)
Σ(p) :=
(1)
y(t) = C(p)x(t)

The transfer function of the parametric reduced model is

−1
Gr (s, p) = Cr (p) sEr (p) − Ar (p) Br (p)
(7)
If W = V , this is called an (ordinary) Galerkin projection.
If W =
6 V , this is a Petrov-Galerkin projection.
III. PARAMETRIC M ODEL ORDER REDUCTION
In this section, the moment-matching method is elaborated first for a non-parametric model. In the second part,
we present a parametric moment-matching method which
separates the frequency variable s and the parameter p.
A. Moment-matching
Suppose that a non-parametric model is given

E ẋ(t) = Ax(t) + Bu(t)
y(t) = Cx(t)

with x(t) ∈ Rn , u(t) ∈ Rm and y(t) ∈ Rq . Now, let s0 ∈ C
be fixed and expand the Taylor series of the transfer function
G(s) of (8) at s0 . That is

where x(t) ∈ Rn , u(t) ∈ Rm and y(t) ∈ Rq denote,
respectively, the state vector, the input and the output. Here
n is the cardinality of the mesh of a spatially distributed
configuration space Ω (usually of dimension 3). The statespace matrices are functions A : R` → Rn×n , E : R` →
Rn×n etc. Before explaining the model reduction, we first
state several assumptions.
Assumption 1: i) [E(p), A(p), B(p), C(p)] system matrices are affine functions of p ∈ R` . e.g., E(p) = E0 +
Σ`i=1 Ei pi ii) p ∈ P ⊂ R` is not varying with time. iii)
For all p ∈ P, all generalized eigenvalues of (λE(p) − A(p))
have negative real parts.
With the time-invariant p, the transfer function of (1) is
meaningful and is given as follows:

−1
G(s, p) = C(p) sE(p) − A(p) B(p).
(2)

G(s)

V ∈ Rn×r

with xr ∈ Rr .

W := im(W );

where s = s0 + σ and the complex matrices mi (s0 ) are
known as the moments of G(s) at expansion point s0 . Define
the block Krylov subspace of order r > 0 by
Kr (M, f )

T

f
M

Cr (p) = C(p)V.

(10)

=
=

(s0 E − A)−1 B ∈ Rn×m
−1

(s0 E − A)

E∈R

n×n

(11)

.

Then Kr (M, f ) has dimension r and the columns of the
matrix V in (3) are constructed as a basis of Kr (M, f ).
Matrix W is chosen such that W T AV is nonsingular. Then
the transfer function Gr (s) of the resulting projection-based
reduced model is obtained by (7) and the reduced system
satisfies the following condition

(3)

(4)

mr,j (s0 ) = mj (s0 ),

j = 0, . . . , r − 1

(12)

Here, the mr,j (s0 ) and mj (s0 ) denote the jth moment of
the reduced model at s0 and the jth moment of the original
model, respectively.
The aforementioned method describes the one-sided moment matching. By choosing the columns of W as a basis
of the Krylov subspace

Ar (p) = W T A(p)V

Br (p) = W B(p),

= colspan{f, M f, · · · , M r−1 f }

where

i.e. two r-dimensional subspaces of Rn . Let ΠV and ΠW
denote the (canonical) projections of Rn onto V and W,
respectively. That is ΠV = V (V T V )−1 V T and ΠW =
W (W T W )−1 W T and consider the system matrices of the
reduced order parametric model defined by
Er (p) = W T E(p)V,

C((s0 + σ)E − A)−1 B
∞
h
ii
X
C −(s0 E − A)−1 E (s0 E − A)−1 B σ i (9)
{z
}
i=0 |
:=mi (s0 )

Define the projection spaces
V := im(V );

=
=

We consider a projection-based method for generating the
reduced order parametric model. That is, let V ∈ Rn×r and
W ∈ Rn×r be full rank matrices with r  n and let
x = V xr ,

(8)

Kr (M2 , l) = colspan{l, M2 l, · · · , M2r−1 l}

(5)
2

(13)

where

below
l

=

M2

=

(s0 E − A)−T C T ∈ Rn×q

(s0 E − A)−T E T ∈ Rn×n

G(s, p) = G(si , pj )
∂G
∂G
+
(si , pj )(s − si ) +
(si , pj )(p − pj )
∂s
∂p
1 ∂2G
1 ∂2G
+
(si , pj )(s − si )2 +
(si , pj )(p − pj )2
2
2∂ s
2 ∂2p


1 ∂2G
∂2G
+
+
(si , pj )(s − si )(p − pj )
2 ∂s∂p ∂p∂s
+ ···
(16)

(14)

another r moments can be matched. The proof can be found
in [11].
Theorem 1: Let the projection matrices V and W be
such that V and W in (4) coincide with (10) and (13),
respectively. Then, the first 2r moments mi (s0 ) and mr,i (s0 ),
i = 0, · · · , 2r − 1 of G(s) and Gr (s) defined in (2) and (7)
coincide.
It is worth to mention that the projection matrices V
and W can have complex numbers with complex expansion
points, as a result, the system matrices of the reduced model
may contain complex numbers which is often physically
unreasonable. One solution is to take s0 real or to have a
pair of complex conjugate expansion points (s0 , s¯0 ) and the
reduced model results in real-valued matrices.

Here we present the Taylor series for a single point.
Applying this for S × P yields multiple expansion points as
depicted in Figure 1. The parametric moments for multiple
expansion points are defined as

i
mi (s, p) : = C(p) − (sE(p) − A(p))E(p) ·
(sE(p) − A(p))−1 B(p)

where (s, p) ∈ S × P.

B. Multi-parameter and multi-frequency moment-matching
G(s 0

The main ingredient of the multi-parameter and multifrequency moment-matching method is to expand the transfer
function into a Taylor series at both a desired frequency and
a desired parameter (s0 , p0 ). This work has been generalized
in [5].
Using (Assumption 1), we obtain
G(s, p)

(17)

, pj)

··
·
·
·
· · · · ·· ·
· · ·· · · · ·
· ·· ··

G( s

P

i, p
0)

s

Fig. 1.

Multi-frequency and multi-parameter expansion

The parametric projection matrices V and W are now
defined as follows.

= C(p)[sE0 + Σ`i=1 sEi pi − A(p)]−1 B(p) (15)

V = Kr (M 1 , F )

where the product of the frequency s and the parameter pi
can be redefined as ŝi = s · pi . This has been proposed in
[5]. Potentially, combining the s with pi as a single variable
ŝi can decrease the complexity of the power series. In such
a way, the computational burden of the projection matrices
V and W is reduced. However, with such conversion the
reduced model matches the first moments of the original
model with respect to ŝi and the physical interpretation of the
parameters p in the reduced model are no longer preserved.
Moreover, the Taylor expansion of the new parametric trans2
fer function includes terms such as ( ∂ ŝi∂∂ ŝGi+1 ) which are
partial derivatives of composed functions rather than the
original one.
We define two sets S := {s1 , . . . , sk } ⊂ C and P :=
{p1 , . . . , p` } ⊂ R which denote the frequency variables and
design parameters of interest, respectively. The projection
matrices V and W can be constructed along these two sets.
The following details the parametric moments which separate the frequency s from the parameters p. The Taylor series
of (15) around the expansion point (si , pj ) ∈ S × P is given

W = Kr (M 2 , L)

(18)
(19)

where
n
ok,`
F = rowspan (si E(pj ) − A(pj ))−1 B(pj )
i=1,j=1
n
ok,`
L = rowspan (si E(pj ) − A(pj ))−T C(pj )T
M 1 = rowspan{(si E(pj ) −

M 2 = rowspan{(si E(pj ) −

i=1,j=1
k,`
A(pj )) E(pj )}i=1,j=1
A(pj ))−T E(pj )T }k,`
i=1,j=1 .
−1

(20)

As in (7) with (V,W) satisfies (4), this yields a reduced
order model whose transfer function
Gr (s, p) = Cr (p)[sE(r,0) + Σ`i=1 sE(r,i) pi − Ar (p)]−1 Br (p)
(21)
coincides in as many moments of the Taylor series of the
original model about (s, p) as possible for a given order of
the reduced model. That is, we claim that
mr,i (s, p) = mi (s, p)
3

(22)

the interval (a, b) ∈ I, and that f, f 0 , · · · , f N −1 all extend
continuously to the interval [a, b]. If in addition

for all (s, p) ∈ (S × P) and i = 0, · · · , 2r − 1.
Theorem 2: Let the projection matrices V and W be such
that V and W in (4) coincide with (18) and (19), respectively.
Then, the first 2r moments mr,i (s, p) and mi (s, p), i =
0, · · · , 2r − 1 of G(s, p) and Gr (s, p) defined in (16) and
(21) coincide for all (s, p) ∈ S × P.

f (a) = f 0 (a) = · · · = f N −1 (a) = 0 = f (b)
Then there exists such c ∈ (a, b) that

f N (c) = 0.
(29)
The Theorem 3 can be proved by applying Proposition 1
and the proof is omitted due to lack of space.
Furthermore, we can find the upper bound of the remainder
term in Theorem 3 using Taylor’s inequality.
Theorem 4: Under the conditions of Theorem 3, let Ia ⊆
I be any interval that has a in its interior. We define

IV. E RROR ESTIMATION
Guarantees on the fidelity of the reduced model remain
an important question for the previously described moment
matching procedure. One major issue is establishing a prior
estimation of the approximation error. In the following we
will derive an error estimate between G(s, p) and Gr (s, p)
which provide a local accuracy at points (s0 , p0 ) ∈ S ×P and
their neighborhood. That is, the error estimate of G(s, p) −
Gr (s, p) will have local validity in (s0 , p0 ) and is based on
the truncation order of the Taylor series.

M := sup

ξ∈Ia

A. The remainder in Taylor series

f (b) − f (a)
.
(23)
b−a
The above equation is the Mean Value Theorem (MVT)
which is a fundamental theorem of calculus. Taylor formula
can be viewed as a generalization of MVT. In particular, let
f be a function such that f and its first N + 1 derivatives
are continuous at a, then the Taylor series is
f 0 (a)
f (N ) (a)
(x − a) + · · · +
(x − a)N + · · ·
1!
N!
(24)
and the N th-degree Taylor polynomial is obtained

i=1

(x − a)i

(25)
fN (x)

= f (a) +

Thus we can write Taylor’s theorem
f (x) = fN (x) + RN (x)

(31)

holds for all x ∈ Ia .
Hence an upper bound of the error between the N th
degree Taylor approximation fN and f is obtained and is a
prior bound, only depending on f itself and the considered
interval I. In the following, we generalize Theorem 4 to
multivariable function f : C` → Cm×p for MIMO systems
with ` expansion points, m inputs and q outputs.
Suppose f : C` → Cm×p is N + 1 times differentiable on
a set I ⊂ C` that contains a vector a ∈ C` in its interior.
Let Ia ⊂ I be a closed and bounded subset of I. Then the
truncated N th order Taylor polynomial around a ∈ C` at
order N is

f (x) = f (a) +

i!

(30)

|f (x) − fN (x)| = |RN (x)| ≤ M|x − a|N +1

f 0 (c) =

N
X
f i (a)

f (N +1) (ξ)
(N + 1)!

Then

Let f : R → R be differentiable on the open interval
(a, b). Then a point c with a < c < b exists such that

fN (x) = f (a) +

(28)

+

(26)

where RN (x) is the remainder term which quantifies the
difference between the Taylor series and the N th-degree
Taylor polynomial in terms of the magnitude of the (N +1)st
derivative of f .
Theorem 3: Suppose that f : R → R is defined on an
interval I that has a in its interior and f (N +1) exists on the
same interval for N ∈ Z+ . Then for each x 6= a, x ∈ I,
there is a ξ between x and a so that

`
X
δi
f (x)
δxi
i=1

x=a

(xi − ai ) + · · ·

δ n1 δ n2 · · · δ n`
(x1 − a1 )n1
f (x)
δxn1 1 δxn2 2 · · · δxn` `
x=a
· · · (x`−1 − a`−1 )n`−1 (x` − a` )n`
(32)

P`
where i=1 ni = N and where n = (n1 , . . . , n` ) is a multiindex of non-negative integers ni
any such multi-index k = (k1 , · · · , k` ) we let |k| =
PFor
`
k
i i be its cardinality. To simplify notation in (32 ), the
kth partial derivative of f at a point a ∈ C` is denoted by
[D|k| f ]k :=

(N +1)

f
(ξ)
(x − a)N +1 .
(27)
(N + 1)!
Thus, there exists a ξ ∈ I, with either x ≤ ξ ≤ a or
a ≤ ξ ≤ x such that the remainder RN (x) in (26) satisfies
(27).
A key observation is that when N = 1, this reduces to the
ordinary MVT. Following this strategy, the key is to observe
the generalization of Rolle’s theorem.
Proposition 1: Suppose f : R → R is a function on
interval I. Assume that f is N ∈ Z+ times differentiable on
RN (x) =

∂ k1 · · · ∂ k`

∂xk11 · · · ∂xk` `

f (x)

.

(33)

x=a

When ranging over all multi-indices k with |k| = K,
[D|k| f ]k can be interpreted as the coefficients of a order
|k| tensor (i.e. multi-linear functional)
[D|k| f ] : C` × C` × · · · × C` → C
{z
}
|
|k| copies

that has (33) in its elements.
4

(34)

Thus, D0 f (a) is a number. D1 f = ∇f (a) is the gradient
of f at the point a. D2 f = ∇2 f (a) is the Hessian matrix of
f at the point a and
[DN f ](v1 , · · · , vN ) :=

N
X

the Theorem 5 gives a prior upper bound on the reduction
error using multi-parameter and multi-frequency moment
matching method. Furthermore, this error bound allows for
finding the optimal expansion points and truncated order.

[DN f ]k (v1 ⊗ · · · ⊗ vN ) (35)

V. E XAMPLE : A THERMAL PARAMETERIZED MODEL
In this section, we demonstrate some numerical results
of the proposed method. To illustrate the procedure and
technique, we consider the example of a linear motor (see
Figure 2). The coils of the linear motor generate a magnetic

|k|=0

With the tensor definition of [D0 f ], [D1 f ], · · · , [DN f ],
the N th order multivariable of f at a ∈ C` (32) can be
rewritten as

fN (x) = [D0 f ] + [D1 f ] [x − a] + · · ·


1
+ [D|N | f ] [x − a], · · · , [x − a] .(36)
|
{z
}
N!

Input
Output

N copies

Define the error by the residual which has the similar form
as the single variable one (26)
RN (x) := f (x) − fN (x)

(37)
Fig. 2.

and also define k = (k1 , · · · , k` )
k

(x − a) :=

`
Y

flux that operate along the track. Meanwhile, the coil also
creates a heat flux which influences the temperature and,
with that the stability of the thermal system. The aim is to
predict the transient thermal behavior of the linear motor.
Let Ω ⊂ R3 be the spatial configuration space of the motor.
Consider the heat equation on Ω as given by

kj

j=1

(xj − aj ) .

Let, as (30)
M := sup

ξ∈Ia

D(N +1) f (ξ)
(N + 1)!

(38)

cp ρ · Ṫ

T (0)

where the norm is the Frobenius norm of a multi-linear
operator and the supremal is taken over expansion points
ξ ∈ Ia ⊂ I. We claim that
N +1 copies

(39)
Theorem 5: Suppose H : C` → Cm×p is N + 1
times differentiable on a set I ⊂ C` that contains vectors
(s0 , p0 ) ∈ (S × P) ⊆ C` in its interior. Let I(s0 ,p0 ) ⊂ I be
a closed and bounded subset of I. Then




s − s0
s − s0
kH(s, p) − HN (s, p)k ≤ M
⊗ ··· ⊗
p − p0
p − p0
{z
}
|
(40)
holds for all (s, p) ∈ I(s0 ,p0 ) , where
(ξs ,ξp )∈I(s0 ,p0 )

D(N +1) H(ξs , ξp )
(N + 1)!

∆ · (κ∆T ) + q

= T0

on ∂Ω1

in Ω

(42)
(43)

where, with some abuse of notation, T (t) ∈ Rn , n = 1560
and where p denotes the heat capacity of the material and it
is characterized by an affine function E(p) = E0 + Cp E1 .
The state T represents spatial-temporal information which
includes T (t) = col[T (x1 , t), · · · , T (x1560 , t)].
The design goal for this linear motor is to find a suitable
material which satisfies the required operating temperature.
Therefore, we illustrate three different heat capacities and the
Cp varies from 200[J/kg · K] to 900[J/kg · K]. Meanwhile,
the operating frequency range of the thermal behavior is
mainly located at low frequency (near 1 ∼ 10 Hz). Under
these conditions thereupon, we choose frequencies s =
[0, 1 ± i] and parameters p = [200, 500, 900].
Following the procedure in Section III, we obtain the reduced model with r = 12 orders which includes 6 expansion

N +1

sup

=

where T is the spatial-temporal dependent temperature of
(i.e. T : Ω × [0, ∞) → R ) and where the initial temperature
T0 (w) := T (w, 0) is given on the boundary w ∈ ∂Ω of Ω. cp
is the material specific heat capacity, ρ denotes the material
density and κ describes the thermal conductivity. The heat
flux q : Ω → R is spatially distributed and is considered as
the input of the system.
Using a finite element (FE) discretization of Ω and parameterizing the material properties, a parameter-dependent
system is generated

E(p)Ṫ = AT + Bu
,
(44)
y = CT

kf (x) − fN (x)k ≤ M [x − a] ⊗ · · · ⊗ [x − a] ; ∀x ∈ Ia
|
{z
}

M :=

Linear motor with inputs and outputs

(41)

is the supremum over the induced norm of the tensor
[D|n+1| H] over I(s0 ,p0 ) .
Hence, (40) is the error bound for multivariable function for
MIMO system.
The function H(s, p) and HN (s, p) in (40) can be interpreted as the (2) and (7) at points (s0 , p0 ) and N thorder reduced model at points (s0 , p0 ), respectively. With
the desired reduced order N and expansion points (s0 , p0 ),
5

points [0, 200], [0, 500], [0, 900], [1 ± i, 200], [1 ± i, 500], [1 ±
i, 900]. We preserve two moments for each point and obtain
the reduced model as

Er (p)ẋr = Ar xr + Br u
, for xr ∈ R12
(45)
y = C r xr

and same number of moments as the proposed method).
Obviously, the proposed method has the better performance
which is very close to the reference. The previous method
also can match the most of the thermal behavior except
for the transient parts. However, the simulation we perform
here is only 10 seconds, it has been observed that the
longer simulation time and complicated input are chosen, the
performance of the proposed method will not be as good as
this one with the same expansion points. Therefore, it leads to
another challenging topic which has been left without enough
attention, how to choose the optimal expansion points with
or without prior knowledge.

here Er (p) = Er,0 + Cp Er,1 . As demonstrated in Figure (3),
the solid lines represent the full model with different values
of Cp . The dashed lines denote the reduced order model with
corresponding p. The frequency range [0.1, 10] in the Bode
plot for all three models are matched, and there are some
mismatches from 20 Hz.

VI. C ONCLUSION
This paper considers the moment matching problem for
linear parameterized system. An explicit algorithm is proposed that establishes a perfect match of moments in both
frequency points S and parameter samples P. An error bound
has been derived that represents the local accuracy of the
transfer function of the reduced order model nearby the
points S × P. Simulation results from both time domain and
frequency domain show that the proposed method delivers
good matching and outperforms the previous work.

Bode Diagram of full model and reduced with different Cp
20

Magnitude (dB)

0
-20
-40
Full model_Cp_200
Full model_Cp_500
Full model_Cp_900
Reduced model_Cp_200
Reduced model_Cp_500
Reduced model_Cp_900

-60
-80
-100
-120
-45

Phase (deg)

-90
-135
-180

R EFERENCES

-225
-270
10

-2

10

-1

0

10
10
Frequency (rad/s)

1

10

2

10

[1] A. C. Antoulas. On recursiveness and related topics in linear systems.
IEEE Transactions on Automatic Control, 31(12):1121–1135, 1986.
[2] A. C. Antoulas. Approximation of large-scale dynamical systems.
SIAM, 2005.
[3] P. Benner, S. Gugercin, and K. Willcox. A survey of projection-based
model reduction methods for parametric dynamical systems. SIAM
review, 57(4):483–531, 2015.
[4] D. L. Boley. Krylov space methods on state-space control models.
Circuits, Systems, and Signal Processing, 13(6):733–758, 1994.
[5] L. Daniel, O. C. Siong, L. S. Chay, K. H. Lee, and J. White. A multiparameter moment-matching model-reduction approach for generating
geometrically parameterized interconnect performance models. IEEE
Transactions on Computer-Aided Design of Integrated Circuits and
Systems, 23(5):678–693, 2004.
[6] P. Feldmann and R. W. Freund. Efficient linear circuit analysis
by padé approximation via the lanczos process. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems,
14(5):639–649, 1995.
[7] P. Feldmann and R. W. Freund. Reduced-order modeling of large
linear subcircuits via a block lanczos algorithm. In Proceedings of
the 32nd annual ACM/IEEE Design Automation Conference, pages
474–479. ACM, 1995.
[8] S. Gugercin, A. C. Antoulas, and C. Beattie. H2 Model Reduction
for Large-Scale Linear Dynamical Systems. SIAM Journal on Matrix
Analysis and Applications, June 2008.
[9] H. Panzer, J. Mohring, R. Eid, and B. Lohmann. Parametric model
order reduction by matrix interpolation. at-Automatisierungstechnik
Methoden und Anwendungen der Steuerungs-, Regelungs-und Informationstechnik, 58(8):475–484, 2010.
[10] L. T. Pillage and R. A. Rohrer. Asymptotic waveform evaluation for
timing analysis. IEEE Transactions on Computer-Aided Design of
Integrated Circuits and Systems, 9(4):352–366, 1990.
[11] B. Salimbahrami and B. Lohmann. Order reduction of large scale
second-order systems using krylov subspace methods. Linear Algebra
and its Applications, 415(2-3):385–405, 2006.
[12] L. M. Silveira, M. Kamon, and J. White. Efficient reduced-order
modeling of frequency-dependent coupling inductances associated
with 3-d interconnect structures. IEEE Transactions on Components,
Packaging, and Manufacturing Technology: Part B, 19(2):283–288,
1996.
[13] D. Weile, E. Michielssen, E. Grimme, and K. Gallivan. A method for
generating rational interpolant reduced order models of two-parameter
linear systems. Applied Mathematics Letters, 12(5):93–102, 1999.

3

Fig. 3.
Bode diagram of full model and reduced model with Cp =
[200, 500, 900][J/kg · K].

To show the benefits of the separation of frequencies and
parameters in the construction of the projection matrices
V and W , we also provide a comparison of time domain
simulations for the proposed method and the method in [5],
the COMSOL simulation results which are considered as
reference are also included. By applying the same input, the
Figure (4) shows the outputs of all 9 models with different
heat capacities.
Comparsion of proposed method and conventional method

Temperature [K]

310

Full model Cp 200
Full model Cp 500
Reduced model Cp 200
Reduced model Cp 500
Reduced model Cp 900
Full model Cp 900
Conv. method Cp 200
Conv. method Cp 500
Conv. method Cp 900

305

300

295

0

1

2

3

4

5

6

7

8

9

10

Time [s]

Fig. 4. Comparison of full model and reduced model with different heat
capacities Cp = [200, 500, 900][J/kg · K] in time domain simulation.

In Figure (4), the solid lines represent the reference of
all models with different p. The dash lines are the results
of the proposed method, the dot lines are the methods
from previous work (we choose the same expansion points
6

