IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN OF INTEGRATED CIRCUITS AND SYSTEMS, VOL. 25, NO. 2, FEBRUARY 2006 285
Macromodel Generation for BioMEMS Components
Using a Stabilized Balanced Truncation Plus
Trajectory Piecewise-Linear Approach
Dmitry Vasilyev, Student Member, IEEE, Michał Rewie
´
nski, Member, IEEE, and
Jacob White, Associate Member, IEEE
Abstract—In this paper, we present a technique for automat-
ically extracting nonlinear macromodels of biomedical micro-
electromechanical systems devices from physical simulation. The
technique is a modification of the recently developed trajectory
piecewise-linear approach, but uses ideas from balanced trunca-
tion to produce much lower order and more accurate models. The
key result is a perturbation analysis of an instability problem with
the reduction algorithm, and a simple modification that makes the
algorithm more robust. Results are presented from examples to
demonstrate dramatic improvements in reduced model accuracy
and show the limitations of the method.
Index Terms—Biomedical microelectromechanical devices
(bioMEMS), microelectromechanical devices (MEMS), model
order reduction, nonlinear dynamical systems, perturbation
methods, piecewise linear models.
I. INTRODUCTION
T
HE application of micromachining to biological applica-
tions, such as labs-on-a-chip [1]–[3], require complicated
combinations of individual biomedical microelectromechanical
systems (bioMEMS) devices that process fluids, cells, and
molecules (e.g., mixers, separators, and pumps). In order to
simulate systems of these devices, models have been devel-
oped for common components, such as mixers and separators
[4]–[7]. The wide variety of devices currently in development,
and the need to rapidly assess the impact of candidate device
performance on system behavior, will accelerate the demand
for techniques that more automatically extract models of these
bioMEMS devices from detailed physical simulation. The re-
quired automatic techniques may include approaches similar to
the robust nonlinear model-order reduction (MOR) strategies
being developed for nonlinear circuit model reduction [8]–[12],
though bioMEMS devices can be more challenging because
Manuscript received March 5, 2005; revised June 28, 2005. This work was
supported by National Science Foundation, by the Defense Advanced Research
Projects Agency (DARPA) under the NeoCAD program, by the Semiconduc-
tor Research Corporation, and by the Singapore–Massachusetts Institute of
Technology (MIT) alliance. This paper was recommended by Associate Editor
J. Zeng.
D. Vasilyev and J. White are with the Research Laboratory of Electronics,
Massachusetts Institute of Technology, Cambridge, MA 02139 USA (e-mail:
vasilyev@mit.edu; white@mit.edu).
M. Rewie
´
nski was with the Research Laboratory of Electronics, Massachu-
setts Institute of Technology, Cambridge, MA 02139 USA. He is now with
Synopsys Inc., Mountain View, CA 94043 USA (e-mail: Michal.Rewienski@
synopsys.com).
Digital Object Identifier 10.1109/TCAD.2005.857389
they are both nonlinear and typically much less damped than
circuits.
In this short paper, we describe an effective model reduction
algorithm for bioMEMS devices that is a modification of the
trajectory piecewise-linear (TPWL) MOR algorithm [10]. In
the following section, we describe the TPWL MOR algorithm,
and then in Section III, we present an improvement on that
algorithm based on using truncated balanced realization (TBR)
[13]. In Section IV, we describe several example problems, and
in Section V, we use those examples to demonstrate both the
effectiveness of our TPWL-TBR algorithm, as well as an insta-
bility problem. Also, we demonstrate a fundamental difficulty
of the TPWL approach when applied to traveling-wave prob-
lems. In Section VI, we describe a perturbation analysis of the
instability problem, and give a second algorithm modification
that resolves this problem. Conclusions end the paper.
II. TPWL N
ONLINEAR MODEL REDUCTION
After spatial discretization of the coupled partial differential
equations (PDEs) that describe a bioMEMS component, the
dynamic behavior of the component can often be represented
using the standard state-space form
!
˙x(t)=f (x(t),u(t))
y(t)=Cx(t)
(1)
where x(t) ∈ R
N
is a vector of states (e.g., mechanical dis-
placements, fluid velocities) at time t, f : R
N
× R
M
→ R
N
is a nonlinear vector-valued function, u : R → R
M
is an input
signal, C is an N × K output matrix and y : R → R
K
is the
output signal.
We assume nonlinear function f being differentiable for all
values of x and u:
f(x, u)=f(x
0
,u
0
)+A(x − x
0
)+B(u − u
0
)+h.o.t. (2)
where matrices A and B [which are dependent on the lineariza-
tion point (x
0
,u
0
)] contain derivatives of f with respect to the
components of the state and input signals, respectively.
The goal of applying MOR to (1) is to construct a macro-
model capable of approximately simulating the input–output
behavior of the systems in (1), but at a significantly reduced
computational cost.
0278-0070/$20.00 © 2006 IEEE
286 IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN OF INTEGRATED CIRCUITS AND SYSTEMS, VOL. 25, NO. 2, FEBRUARY 2006
In order to achieve this goal, one needs to reduce the dimen-
sionality of the state-space vector, which is usually achieved by
employing projections. However, only projecting the nonlinear
system (1) is not a complete solution to the nonlinear model
reduction problem, because direct evaluation of the projected
function is still directly proportional to the size of the unreduced
system, and is too computationally expensive.
1
To reduce the cost of the projected-function evaluation, con-
sider the following generalized quasi-piecewise-linear approxi-
mate representation of the nonlinear function f, which has been
proposed, in a slightly simpler form, in [10]
f(x, u) ≈
s
"
i=1
˜w
i
(x, u)(f(x
i
,u
i
)
+ A
i
(x − x
i
)+B
i
(u − u
i
)) (3)
where x
i
’s and u
i
’s (i =1,...,s) are selected linearization
points (samples of state and input values), A
i
and B
i
are
derivatives of f with respect to x and u , evaluated at (x
i
,u
i
),
and finally, ˜w
i
(x, u)’s are state- and input-dependent weights
that satisfy
s
"
i=1
˜w
i
(x, u)=1 ∀(x, u),
˜w
i
(x, u) → 1 as (x, u) → (x
i
,u
i
). (4)
Equation (4) implies that the TPWL approximation in (3)
is simply a convex combination of samples of f and f’s
derivatives.
Projecting the piecewise-linear approximation in (3) using
biorthonormal projection bases V and W yields the following
reduced-order nonlinear dynamical system:
˙z = γ · w(z, u)+(
&
s
i=1
w
i
(z, u)A
ir
) z
+(
&
s
i=1
w
i
(z, u)B
ir
) u
y = C
r
z
(5)
where z(t) ∈ R
q
is the q-dimensional vector of states
γ =
'
W
T
(f(x
1
,u
1
) − A
1
x
1
− B
1
u
1
) ...
W
T
(f(x
s
) − A
s
x
s
− B
s
u
s
)
(
.
Here, w(z, u)=[w
1
(z, u) ···w
s
(z, u)]
T
is a vector of weights,
A
ir
= W
T
A
i
V , B
ir
= W
T
B
i
, and C
r
= CV . One should
note that
&
s
i=1
w
i
(z, u)=1for all (z, u), w
i
→ 1 as (z, u) →
(W
T
x
i
,u), and that the evaluation of the right-hand side of (5)
requires at most O(sq
2
) operations, where s is the number of
linearization points.
As proposed in [10] and [12], linearization points (x
i
,u
i
)
used in system (5) are usually selected from a “training tra-
jectory” of the initial nonlinear system, corresponding to some
appropriately determined “training input.” The choice of the
training input is an important aspect of the reduction procedure,
since this choice directly influences accuracy. As the general
rule, the training signal should be as close as possible to the
1
This issue is discussed in detail in [10].
signals for which the reduced system will be used. Additionally,
this input signal should be rich enough to collect all “important”
states in the set of linearization points (x
i
,u
i
) [12].
In order to obtain a reduced system in the form of (5),
biorthonormal projection bases V and W must also be deter-
mined. This issue is addressed below.
III. C
HOICE OF LINEAR REDUCTION METHOD
Consider a simple linearization of (1) about the initial state
(x
0
,u
0
)
!
˙x = f(x
0
) −
ˆ
Ax
0
−
ˆ
Bu
0
+
ˆ
Ax +
ˆ
Bu
y =
ˆ
Cx
. (6)
For the system in (6), a projection basis can be obtained using
one of the many projection-based linear MOR procedures. One
common choice is to reduce using the projection basis spanning
the Krylov subspace [10], [14], [15]
span(V )=span{
ˆ
A
−1
ˆ
B, . . . ,
ˆ
A
−q
ˆ
B}.
Reduction of (6) using the Krylov-subspace projection is not
guaranteed to provide a stable reduced model, even in this
linearized case [16], [17]. Therefore, TPWL macromodels ob-
tained using Krylov projection are not guaranteed to be stable
even if the original system is nearly linear.
Alternatively, one can apply a balanced truncation model
reduction (TBR) procedure [18]–[20]:
Algorithm 1: TBR
Input: System matrices
ˆ
A,
ˆ
B, and
ˆ
C.
Output: Projection bases V and W .
(1) Find observability gramian P :
ˆ
AP + P
ˆ
A
T
= −
ˆ
B
ˆ
B
T
;
(2) Find controllability gramian Q:
ˆ
A
T
Q+ Q
ˆ
A = −
ˆ
C
T
ˆ
C;
(3) Compute q dominant eigenvectors of PQ: (PQ)V =
V Σ
2
, where Σ
2
= diag(Λ
dom
q
(PQ));
(4) Compute q dominant eigenvectors of QP :
(QP )W = W Σ
2
and scale columns of W such that
W
T
V = I
q ×q
.
The projection bases V and W obtained using algorithm 1
can then be used to compute the reduced TPWL approximation
in (5).
TBR reduction can be more accurate than Krylov-subspace
reduction as it possesses a uniform frequency error bound [21],
and TBR preserves the stability of the linearized model. This
superior performance for the linear cases suggests that TPWL
approximation models obtained using TBR will be stable and
accurate as well. This is not necessarily the case, as will be
shown below.
IV. E
XAMPLES OF NONLINEAR SYSTEMS
In this section, we consider two examples of nonlinear sys-
tems that arise in the modeling of bioMEMS devices that have
nonlinear dynamical behaviors, which make good test cases for
reduction algorithms.
VASILYEV et al.: MACROMODEL GENERATION FOR bioMEMS USING TBR-BASED TPWL 287
Fig. 1. Micropump example (following Hung et al. [22]).
The first example is a fixed–fixed beam structure, which
might be used as part of a micropump or valve, shown in
Fig. 1. Following Hung et al. [22], the dynamical behavior
of this coupled electromechanical-fluid system can be mod-
eled with a one-dimensional (1-D) Euler’s beam equation and
the two-dimensional (2-D) Reynolds’ squeeze film damping
equation [22]:
ˆ
EI
∂
4
w
∂x
4
− S
∂
2
w
∂x
2
= F
elec
+
*
d
0
(p − p
0
)dy − ρ
∂
2
w
∂t
2
∇ · ((1 + 6K)w
3
p∇p) = 12µ
∂(pw)
∂t
.
(7)
Here, the axes x, y, and z are as shown on Fig. 1,
ˆ
E is the
Young’s modulus, I is the moment of inertia of the beam, S
is the stress coefficient, K is the Knudsen number, d is the
width of the beam in the y-direction, w = w(x, t) is the height
of the beam above the substrate, and p(x, y, t) is the pressure
distribution in the fluid below the beam. The electrostatic force
is approximated assuming nearly parallel plates and is given by
F
elec
= $
0
dv
2
/2w
2
, where v is the applied voltage.
Spatial discretization of (7) using a standard finite-difference
scheme leads to a nonlinear dynamical system in the form of
(1), with N = 880 states. After discretization, the state vector
x consists of the concatenation of: heights of the beam above
the substrate w, values of ∂(w
3
)/∂t, and values of the pressure
below the beam. For the considered example, the output y(t)
was selected to be the deflection of the center of the beam from
the equilibrium point (cf. Fig. 1).
The remarkable feature of this example is that the system is
strongly nonlinear, and no feasible Taylor expansion made at
the initial state can correctly represent the nonlinear function f,
especially in the so-called pull-in region
2
[22].
The exact actuation mechanism of the real micropumps may
be quite different than the above simple structure, but this ex-
ample is illustrative in that it combines electrical actuation with
the structural dynamics and is coupled to fluid compression.
We expect model reduction methods that are effective for this
example problem to be extendable to realistic micropumps.
The second example, suggested in [23], is the injection of a
(marker) fluid into a U-shaped three-dimensional microfluidic
channel. The fluid is driven electrokinetically as depicted in
Fig. 2, and the channel has a rectangular cross section of height
d and width w. In this example, the electrokinetically driven
2
If the beam is deflected by more than ≈ 1/3 of the initial gap, the beam will
be pulled- in to the substrate.
Fig. 2. Microfluidic channel.
flow of a buffer (carrier) fluid was considered to be steady, with
the fluid velocity directly proportional to the electric field as in
%v(x, y, z
+
,- .
!r
)=−µ∇Φ(%r)
where µ is the electroosmotic mobility of the fluid. The electric
field can be determined from Laplace’s equation
∇
2
Φ(%r)=0
with Neumann boundary conditions on the channel walls [24].
If the concentration of the marker is not small, the electroos-
motic mobility can become dependent on the concentration,
i.e., µ ≡ µ(C(%r, t)), where C(%r, t) is the concentration of the
marker fluid. Finally, the marker can diffuse from the areas
with high concentration to the areas with low concentration.
The total flux of the marker, therefore, is
%
J = %vC − D∇C (8)
where D is the diffusion coefficient of the marker. Again, as
the concentration of the marker grows, the diffusion will be
governed not only by the properties of the carrying fluid, but
also by the properties of the marker fluid, therefore, D can
depend on concentration. Conservation applied to the flux (8)
yields a convection-diffusion equation [25]
∂C
∂t
= −∇ ·
%
J = ∇Φ · (C∇µ(C)+µ(C)∇C)
+ ∇D(C) · ∇C + D(C)∇
2
C. (9)
The standard approach is to enforce zero normal flux at the
channel-wall boundaries, but since %v has a zero normal com-
ponent at the walls, zero normal flux is equivalent to enforcing
zero normal derivative in C. The concentration at the inlet was
determined by the input, and the normal derivative of C was
assumed 0 at the outlet.
Note that (9) is nonlinear with respect to marker concentra-
tion unless both electroosmotic mobility and diffusion coeffi-
cient are concentration-dependent.
A state-space system was generated from (9) by ap-
plying a second-order three-dimensional coordinate-mapped
288 IEEE TRANSACTIONS ON COMPUTER-AIDED DESIGN OF INTEGRATED CIRCUITS AND SYSTEMS, VOL. 25, NO. 2, FEBRUARY 2006
Fig. 3. Propagation of the square impulse of concentration of the marker.
Due to the difference in lengths of the inner and outer arc, the marker reaches
different points at the outlet with different delay.
finite-difference spatial discretization to (9) on the half-ring
domain in Fig. 2. The states were chosen to be concentrations
of the marker fluid at the spatial locations inside the channel.
The concentration of the marker at the inlet of the channel is
the input signal, and there are three output signals: the first
being the average concentration at the outlet, and the second
and third signals being the concentrations at the inner and outer
radii of the outlet of the channel, respectively.
Fig. 3 illustrates the way an impulse of concentration at
the inlet propagates through the channel: Diffusion spreads the
pulse, and due to the curvature of the channel, the front of
the impulse becomes tilted with respect to the channel’s cross
section. That is, the marker first reaches the points at the inner
radius (point 1).
V. C
OMPUTATIONAL RESULTS
In this section, results are first presented for the linear mi-
crofluid channel model, in order to emphasize the efficiency of
the TBR linear reduction. Then, results are presented for the
micromachined-pump model. The most challenging example
was a nonlinear microfluidic channel.
A. Microchannel—Linear Model
First, in order to demonstrate the effectiveness of TBR linear
reduction, consider applying a balanced-truncation algorithm
to the linearized microchannel model. This corresponds to
the problem of a very diluted solution of a marker in the
carrier liquid (a widely used approximation in the literature).
The values used for the electroosmotic mobility and diffusion
coefficients are from [23]: µ =2.8 × 10
−8
m
2
· V
−1
· s
−1
, D =
5.5 × 10
−10
m
2
· s
−1
. Physical dimensions of the channel were
chosen to be r
1
= 500 µm, w = 300 µm, and d = 300 µm.
Finite-difference discretization led to a linear time-invariant
system (A, B, C) of order N = 2842 (49 discretization points
by angle, 29 by radius, and 2 by height). Since algorithm 1
requires O(n
3
) computation, the discretized system was too
costly to reduce using the original TBR algorithm. Instead,
we used a fast-to-compute approximation to the TBR called
modified approximate implicit subspace iteration with alternate
directions (AISIAD) [26], [27].
Fig. 4. H-infinity errors (maximal discrepancy over all frequencies between
transfer functions of original and reduced models) for the Krylov, TBR, and
modified AISIAD reduction algorithms.
Fig. 5. Transient responses of the original linear (dashed lines) and reduced
(solid lines) model (order q = 13). Input signal: unit pulse with duration
0.1 s. The maximum error between these transients is ≈ 1 × 10
−4
, therefore,
the difference is barely visible. The different outputs correspond to the different
locations along the channel’s outlet (from left to right: innermost point, middle
point, outermost point).
As shown in Figs. 4 and 5, applying reduction to the spatial
discretization of (9) demonstrates the excellent efficiency of
the TBR reduction algorithm. The reduction error decreases
exponentially with increasing reduced model order, both in
frequency- and in time-domain measurements. For example, in
the time-domain simulations, the maximum error in the unit
step response for the reduced model of order q = 20 (over a
100-times reduction) was lower than 10
−6
for all three output
signals.
The modified AISIAD method was compared with Krylov-
subspace-based reduction (Arnoldi method [14]) and the origi-
nal TBR method in both time and frequency domains. As shown
in Fig. 4, TBR and modified AISIAD are much more accurate
than the Krylov method, and are nearly indistinguishable. How-
ever, the modified AISIAD model is much faster to compute.
To demonstrate the time-domain accuracy of the reduced
model, we first redefined the outputs of the model as concen-
trations at the points 1–3 on Fig. 3, and then performed approx-
imate TBR reduction using the modified AISIAD method.
VASILYEV et al.: MACROMODEL GENERATION FOR bioMEMS USING TBR-BASED TPWL 289
Fig. 6. Errors in output computed by TPWL models generated with differ-
ent MOR procedures (micromachined-pump example); N = 880; 5.5-V step
testing and training input voltage.
In Fig. 5, the output produced by a 0.1-s unit pulse is shown.
The results for the 2842-state model and modified AISIAD
reduced model of order 13 are compared in Fig. 5. One can
clearly see that the reduced model nearly perfectly represents
different delay values and the spread of the outputs.
B. Micromachined-Pump Example
The TBR TPWL MOR strategy was applied to generate
macromodels for the micromachined-pump example described
in Section IV. The reduced basis was generated using the
linearized model of system (1) only at the initial state, and the
initial state was included in the bases V and W .
Surprisingly, unlike in several nonlinear circuit examples
[13], the output error did not decrease monotonically as the
order q of the reduced system grew. Instead, macromodels with
odd orders behaved very differently than macromodels with
even orders. Models of even orders were substantially more
accurate than models of the same order generated by Krylov
reduction—cf., Fig. 6. However, if q was odd, inaccurate and
unstable reduced-order models were obtained. This phenom-
enon is reflected in the error plot shown in Fig. 6. Fig. 7
illustrates that a fourth-order (even) reduced model accurately
reproduces transient behavior.
This “even–odd” phenomenon was observed in [28] and
explained in the very general sense in [29]. The main result
in [29] is described in the following section. However, there
is also an insightful but less general way of looking at this
phenomenon.
The “even–odd” phenomenon can be viewed by examining
eigenvalues of the reduced-order Jacobians from different lin-
earization points. For the pump example, the initial nonlinear
system is stable and Jacobians of f at all linearization points
are also stable. Nevertheless, in this example, the generated
reduced-order basis provides a truncated balancing transforma-
tion only for the linearized system from the initial state x
0
.
Therefore, only the reduced Jacobian from x
0
is guaranteed to
Fig. 7. Comparison of system response (micromachined-pump example)
computed with both nonlinear and linear full-order models, and TBR TPWL
reduced-order model (seven models of order q =4); 5.5-V step testing and
training input voltage. Note: solid and dashed lines almost fully overlap.
be stable. Other Jacobians, reduced with the same projection
bases, may develop eigenvalues with positive real parts.
Fig. 8 shows spectra of the reduced-order Jacobians for
models of order q =7and q =8. One may note that, for q =8,
the spectra of the Jacobians from a few first linearization points
are very similar. They also follow the same pattern: two of the
eigenvalues are real, and the rest form complex-conjugate pairs.
Increasing or decreasing the order of the model by 2 creates
or eliminates a complex-conjugate pair of stable eigenvalues
from the spectra of the Jacobians. If the order of the model is
increased or decreased by 1 [cf., Fig. 8 (left)], the situation is
very different. A complex-conjugate pair will be broken, and
a real eigenvalue will form. At the first linearization point,
this eigenvalue is a relatively small negative number. At the
next linearization point, the corresponding eigenvalue shifts
significantly to the right half plane to form an unstable mode
of the system. An obvious workaround for this problem in
the considered example is to generate models of even order.
Nevertheless, a true solution to this problem would involve
investigating how perturbations in the model affect the balanced
reduction, and this is examined in Section VI.
C. Nonlinear Microfluidic Example
Consider introducing a mild nonlinearity into the mobility
and diffusion coefficients in (9) as follows:
µ(C) = (28 + C · 5.6) × 10
−9
m
2
· V
−1
· s
−1
D(C) = (5.5+C · 1.1) × 10
−10
m
2
· s
−1
.
Our experiments showed that even such a small nonlinearity
creates a challenging problem for the TPWL algorithm. For this
problem, the choice of training input significantly affects the
set of the input signals for which the reduced model produces
accurate outputs. For the case of a pulsed marker, this example
has, in effect, a traveling-wave solution. Therefore, linearizing
![Fig. 1. Micropump example (following Hung et al. [22]).](/figures/fig-1-micropump-example-following-hung-et-al-22-2kf8wxd3.webp)
