Eecs 2012 217
Eecs 2012 217
Eecs 2012 217
Systems
Chenjie Gu
December 1, 2012
Committee in charge:
Professor Jaijeet Roychowdhury, Chair
Professor Robert Brayton
Professor Jon Wilkening
Fall 2011
Copyright 2011
by
Chenjie Gu
1
Abstract
2
The third method, NTIM, derives a macromodel that specifically captures timing/phase
responses of a nonlinear system. We rigorously define the phase response for a non-autonomous
system, and derive the dynamics of the phase response. The macromodel emerges as a scalar,
nonlinear time-varying differential equation that can be computed by performing Floquet
analysis of the full model. With the theory developed, we also present efficient numerical
methods to compute the macromodel.
The fourth method, DAE2FSM, considers a slightly different problem finite state machine abstraction of continuous dynamical systems. We present an algorithm that learns a
Mealy machine from a set of differential equations from its input-output trajectories. The
algorithm explores the state space in a smart way so that it can identify the underlying finite
state machine using very few information about input-output trajectories.
Acknowledgments
I would like to thank all those who supported me throughout my doctoral studies.
First and foremost, I would like to thank my advisor Professor Jaijeet Roychowdhury.
Without his guidance, this dissertation would have been impossible. I really appreciate
his insight into numerous aspects in numerical simulation and analog/RF computer-aided
design, as well as his sharpness, enthusiasm, care and encouragement. Learning from and
working with him was an invaluable and unique experience.
I would also like to thank dissertation and examination committee members, Professors Robert Brayton, Jon Wilkening, Andreas Kuehlmann. They have spent a lot of time
examining and guiding my research directions, and have provided valuable comments and
suggestions along my research in the last two years of my doctoral studies.
In addition, I would like to express my special thanks to several other professors: Professor
Alper Demir for many fruitful discussions on phase macromodeling, Professor Sanjit Seshia
for his insights into verification and model checking, Professor Claire Tomlin for her series of
valuable lectures in linear, nonlinear and hybrid control systems, Professor Sachin Sapatnekar
and Kia Bazargan for their early introduction to me the field of computer-aided design and
electronic design automation.
I would also like to extend my thanks to those who have been helping me in my career
planning and job applications, including but not limited to Joel Phillips, Sani Nassif, Eric
Keiter, Heidi Thornquist, Noel Menezes, Chirayu Amin, Qi Zhu, Yuan Xie, Yu Wang.
Working in Professor Jaijeet Roychowdhurys group has been a rewarding experience,
partly because of all the wonderful people I have worked with. Especially I would like
to thank Ning Dong, Xiaolue Lai, Ting Mei for their friendly sharing of knowledge and
experiences, as well as many helpful suggestions in research and life. I would also like to
thank Prateek Bhansali with whom I have worked for more than three years, and David
Amsallem who has been a post-doctoral researcher in the group.
The DOP Center (Donald O. Pederson Center for Electronics Systems Design) in Berkeley
is without doubt a legendary research center and a great place to work in. My thanks also
extend to friends in the DOP center, including but not limited to Baruch Sterin, Sayak Ray,
Dan Holcomb, Wenchao Li, Susmit Jha, Bryan Brady, Pierluigi Nuzzo, Alberto Puggelli,
Xuening Sun, Liangpeng Guo, Jia Zou.
My gratitude also goes to friends in other research groups in Berkeley and Minnesota,
including but not limited to Jiashu Chen, John Crossley, Lingkai Kong, Yue Lu, Sriramkumar
Venugopalan, Chao-Yue Lai, Humberto Gonzale, Nicholas Knight, Chunhui Gu, Zhichun
Wang, Weikang Qian, Qunzeng Liu.
ii
Contents
List of Figures
vii
List of Tables
xi
1 Introduction
1.1 Macromodels and Motivating Applications
1.1.1 Mathematical Models . . . . . . . .
1.1.2 Applications . . . . . . . . . . . . .
1.2 Problem Overview . . . . . . . . . . . . .
1.2.1 Challenges . . . . . . . . . . . . . .
1.3 Dissertation Contributions and Overview .
1.4 Notations and Abbreviations . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
MOR
. . . .
. . . .
. . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
1
2
4
5
5
6
7
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
9
9
9
10
16
17
17
18
18
19
19
20
20
21
21
21
22
23
iii
2.6
2.7
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Reduction
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
23
24
24
24
.
.
.
.
.
.
.
.
.
25
25
26
26
28
31
33
33
34
34
.
.
.
.
.
.
.
.
.
.
.
.
.
.
36
36
37
39
41
44
44
45
46
49
50
51
51
52
55
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
56
56
56
57
58
59
61
61
iv
5.4
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
61
62
62
63
65
67
67
71
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
76
76
79
80
82
84
85
87
90
91
92
92
94
95
95
97
99
100
101
101
102
102
105
.
.
.
.
.
.
.
110
110
114
115
116
119
120
121
5.5
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
v
7.5
7.6
7.7
7.8
Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . .
Numerical Methods for Computing the Phase Model . .
7.6.1 Floquet Decomposition using Monodromy Matrix
7.6.2 Equations in terms of U(t) and V (t) . . . . . . .
7.6.3 Floquet Decomposition via Harmonic Balance . .
7.6.4 Floquet Decomposition via FDTD Method . . . .
7.6.5 Forcing Bi-orthogonality . . . . . . . . . . . . . .
Optimization-Based Phase Macromodel Computation . .
Examples and Experimental Results . . . . . . . . . . . .
7.8.1 A Simple Nonlinear System . . . . . . . . . . . .
7.8.2 A Firing Neuron Model . . . . . . . . . . . . . . .
7.8.3 An Inverter Chain Circuit . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
121
121
122
123
124
126
127
128
128
129
135
136
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
155
156
156
157
158
158
159
161
vi
A Introduction of Differential Geometry
A.1 Vector Fields and Flows . . . . . . . . . . . .
A.2 Lie Bracket and its Properties . . . . . . . . .
A.3 Coordinate Changes and Lie Bracket . . . . .
A.4 Vector Fields as Differential Operators . . . .
A.5 Equivalence between Families of Vector Fields
A.6 Sub-Manifolds and Foliations . . . . . . . . .
A.7 Orbits of Families of Vector Fields . . . . . . .
A.8 Integrability of Distributions and Foliations .
B Benchmark Examples
B.1 Nonlinear Transmission Line Circuit 1 (NTL1)
B.2 Nonlinear Transmission Line 2 (NTL2) . . . .
B.3 Inverter Chain Circuit (INVC) . . . . . . . . .
B.4 Morris-Lecar Neuron Model (NEURON ML) .
B.5 FitzHugh-Nagumo Model (NEURON FN) . .
B.6 Latch Circuit 1 (LATCH1) . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
173
173
174
174
175
176
177
178
179
.
.
.
.
.
.
180
180
181
182
182
183
184
vii
List of Figures
1.1
2.1
2.2
2.3
14
15
15
3.1
27
4.1
4.2
37
38
5.1
Simulation of the ManiMOR model. (u(t) = t) Fig. 5.1(a) shows the manifold
ManiMOR identifies (yellow), and the trajectories of the full model (blue)
and ManiMOR model (green). Fig. 5.1(b) shows the waveform of x3 (t) of the
ManiMOR model (green), the full model (blue), and the projection of the full
model solution onto the manifold (black). . . . . . . . . . . . . . . . . . . .
Simulation of the ManiMOR model. (sinusoidal input) Fig. 5.2(a) shows the
DC curve (red) and the manifold ManiMOR identifies (yellow), the trajectories of the full model (blue) and ManiMOR model (green). Fig. 5.2(b) shows
the waveform of x3 (t) of the ManiMOR model (green), the full model (blue),
and the projection of the full model solution onto the manifold (black). . . .
Comparison of ManiMOR and TPWL model. (multiple-step input) Red,
green and blue trajectories represent simulation results of ManiMOR, TPWL
and full model, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Comparison of ManiMOR and TPWL model. (two-tone sinusoidal input)
Red, green and blue trajectories represent simulation results of ManiMOR,
TPWL and full model, respectively. . . . . . . . . . . . . . . . . . . . . . .
Circuit diagram of a CMOS ring mixer [1]. . . . . . . . . . . . . . . . . . . .
5.2
5.3
5.4
5.5
63
64
65
66
66
viii
5.6
5.7
5.8
5.9
5.10
5.11
5.12
5.13
5.14
5.15
6.1
6.2
6.3
6.4
6.5
6.6
Simulation of ManiMOR model and full model for the CMOS ring mixer. (step
input) Fig. 5.6(a) shows the manifold ManiMOR identifies (yellow), and the
trajectories of the full model (blue) and ManiMOR model (green). Fig. 5.6(b)
shows the waveform of x3 (t) of the ManiMOR model (green), the full model
(blue), and the projection of the full model result onto the constructed manifold (black). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Comparison of of the circuit response of ManiMOR model and TPWL model.
(Nonlinear transmission line circuit [2]) . . . . . . . . . . . . . . . . . . . .
Circuit diagram of a current-model logic (CML) buffer. . . . . . . . . . . . .
Visualization of the low-order manifold subspace generated by ManiMOR. The
axes represent three node voltages in the circuit. The red curve consists of
the DC operating points; the black dashed lines are the second local subspace
bases; the yellow manifold is generated by stitching together all the local linear
subspaces. The blue circled points are on a transient trajectory simulated
using the full model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Singular values of the matrix containing 20 equilibrium points, sorted. . . . .
Comparison of ManiMOR model with linearized reduced-order model. The
waveforms of one output using different models are plotted. . . . . . . . . . .
Comparison of ManiMOR model with TPWL model. The waveforms of one
output voltage using different models are plotted. . . . . . . . . . . . . . . .
Components of the pheromone response pathway in yeast. [3] Each arrow
represents a type of reaction. . . . . . . . . . . . . . . . . . . . . . . . . . . .
Concentrations of x1 at equilibrium corresponding to different inputs, and the
singular values for the DC matrix. . . . . . . . . . . . . . . . . . . . . . . . .
Transient simulation of the yeast pheromone pathway. The concentration of
one reactant is plotted. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
QLMOR flow (yellow), with comparison to MOR based on polynomial form
(green) and bilinear form (gray). . . . . . . . . . . . . . . . . . . . . . . . . .
A simple CMOS circuit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Polynomial approximations of 1/(1 + x). Fig. 6.3(a) shows the interval x
[0, 1]. Fig. 6.3(b) shows the interval x [0, 2]. . . . . . . . . . . . . . . . . . .
Fig. 6.4(a), Fig. 6.4(b) show the time-domain waveforms of x1 , x2 and Fourier
coefficients of x1 , respectively, when u(t) = cos(4t). Fig. 6.4(c), Fig. 6.4(d)
show the time-domain waveforms of x1 , x2 and Fourier coefficients of x1 , respectively, when u(t) = 10 cos(2t). . . . . . . . . . . . . . . . . . . . . . .
Fig. 6.5(a), Fig. 6.5(b) show the time-domain waveforms of x1 , x2 and Fourier
coefficients of x1 , respectively, when u(t) = 8 cos(3t). Fig. 6.5(c), Fig. 6.5(d)
show the time-domain waveforms of x1 , x2 and Fourier coefficients of x1 , respectively, when u(t) = 10 cos(3t). . . . . . . . . . . . . . . . . . . . . . .
Nonlinear transmission line circuit [2]. . . . . . . . . . . . . . . . . . . . . .
67
68
68
69
70
71
72
73
74
75
81
89
103
104
106
107
ix
6.7
6.8
7.1
7.2
7.3
7.4
7.5
7.6
7.7
7.8
7.9
7.10
7.11
7.12
7.13
8.1
8.2
8.3
8.4
8.5
8.6
8.7
8.8
8.9
9.1
B.1
B.2
B.3
B.4
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
141
141
148
149
150
151
152
153
154
181
181
182
184
x
B.5 An inverter and a multiplexer. . . . . . . . . . . . . . . . . . . . . . . . . . . 185
xi
List of Tables
1.1
1.2
1.3
1.4
models .
. . . . .
. . . . .
. . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2.1
2.2
2.3
2.4
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2
4
8
8
. . . . . . . . . . . . .
. . . . . . . . . . . . .
optimization problem
. . . . . . . . . . . . .
.
.
.
.
10
16
23
24
3.1
33
4.1
55
6.1
6.2
6.3
79
83
98
8.1
8.2
. .
. .
the
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Chapter 1
Introduction
1.1
Engineers and scientists build models for many systems in various domains. Computer
simulation of these models predicts system behaviors, and thus helps us understand existing
systems and design new systems.
Typically, we have a hierarchy of models of different levels of abstraction. Higher-level
models abstract away unimportant details of lower-level models, and thus are easier to understand and more efficient to simulate. For example, in digital circuits, we have models at the
gate level describing the gates truth-table and delay/power performance, models at the RTL
level describing the register-level timing and dynamics, models at the finite state machine
level describing the functionality of the circuit, etc.; in analog circuits, we have models at
the device level describing how each electron move, models at the transistor level describing
the current-voltage (I-V) and charge-voltage (Q-V) relationship of a transistor, models at
the sub-circuit level describing sub-circuit performances such as the gain and bandwidth of
an amplifier, models at the system level describing power/functionality of the circuit, etc..
In the process of circuit design, a top-down synthesis approach is commonly employed.
That is, starting with system specifications, we map from higher-level models to lower-level
models iteratively until we reach the transistor-level implementation.
Macromodeling, or model order reduction, is the reverse procedure of synthesis. It refers
to the generation of higher-level models from lower-level models.
Traditionally, macromodeling is done manually. It usually involves a lot of domain knowledge in a specific field, and is rather ad hoc than systematic. As a result, the macromodeling
process takes long time. Moreover, while the manual macromodeling process usually captures responses of major interest, it may miss equally important corner cases, and may create
a defective macromodel.
1.1.1
Mathematical Models
y = h(x).
(1.2)
However, DAE/ODE models (1.1) and (1.2) are most commonly used in many practical
systems, and the techniques we develop easily extends to other forms of DAE/ODE models,
such as the ones listed in Table 1.2.
4
Table 1.2: Other DAE/ODE models
Name
Model
Nonlinear descriptor system dtd q(x) = f (x, u), y = h(x)
Implicit nonlinear system
f ( dtd x, x, u, y) = 0
1.1.2
Applications
The macromodels, or reduced order models, can be used in various contexts to help
understanding and designing a system.
The straightforward application is to speedup simulations of a large system. For example,
circuits nowadays can easily involve millions of transistors and even more number of parasitic
resistors and capacitors. While analog circuits used to be stand-alone and have at most
hundreds of transistors, they are now extensively coupled to digital circuits for feedback and
calibration, thus also leading to a large number of equations. With reduced-order models
of each sub-circuit block, we may replace the original sub-circuits by their reduced-order
models to obtain a smaller system of equations with little accuracy compromise. Therefore
the speedup of the simulation may be significant, and we can afford to simulate larger systems
that are otherwise impossible to simulate within a reasonable amount of time.
Similarly, model order reduction enables a bottom-up approach to create a hierarchy of
models, and to finally estimate higher-level performance metrics of a large system. For example, in aggressive memory designs, device-level PDE simulations are necessary to capture
detailed nonlinear behaviors of transistors. This can involve a large number of equations
by applying finite element methods on the underlying PDEs. Model order reduction can
therefore automatically generate a smaller model for each transistor which can then be used
to estimate circuit performances such as read/write access time of a memory cell.
Reduced order models can also be used in design optimization (e.g., transistor sizing).
In particular, parameterized reduced models [49] capture important system response with
dependence on parameters. Therefore, they simplify the optimization problem and make
the problem computationally more tractable. Some reduced models also naturally carries
sensitivity information with respect to inputs/parameters, and this is very useful in gradientbased optimization methods.
Reduced order models may also help understanding nonlinear systems. To the best of our
knowledge, there does not exist techniques to analyze behaviors of a system modeled by a
large number of differential equations, except for extensive numerical simulation. Bifurcation
theory and systems theory are developed but they are usually applied in very low dimensional
systems (e.g., 2-dimensional). Therefore, rather than a full-size model, a reduced order model
may be systematically analyzed, and can provide useful intuition of how the large system
works.
Since the MOR techniques are developed over the underlying mathematical models, they
are generally applicable to many disciplines, including electronic circuits (interconnects [5,7,
5
1014], electromagnetics [1519], power grids [2022], RF circuits [2328], oscillators [2932]),
micro-electro-mechanical systems [3339] fluid dynamics [4046], biological systems [41, 47
51], neuroscience [5259], chemical engineering [6067], mechanical engineering [6871], etc..
1.2
Problem Overview
Loosely speaking, model order reduction aims to generating simpler models that capture important dynamics of original full models.
Here, model refers to a mathematical model that describes system behaviors. The
model can be linear or nonlinear, continuous or discrete value, continuous or discrete time,
deterministic or stochastic, as shown in Table 1.1.
Simple typically means computationally efficient, but may also represent other properties of the model, such as sparsity of matrices in the model, number of parameters, number
of state variables, correlation among state variables, etc..
What important dynamics mean depends highly on applications. For examples, it
can be transfer functions over a certain frequency range (for LTI systems), the time-domain
waveforms, stability, passivity, etc..
For DAEs (1.1), model order reduction typically means deriving a system of reduced
DAEs
d
qr (xr ) + fr (xr ) + Br u(t) = 0, y = hr (xr ),
(1.3)
dt
where xr RNr are the state variables of the reduced model; qr (xr ), fr (xr ) : RNr RNr ,
hr (xr ) : RNr RNo are functions; B RNr Ni ; Nr is the order of the reduced model.
1.2.1
Challenges
While we will have more detailed discussions about challenges in nonlinear MOR in
Chapter 2 and other chapters, we highlight the major ones in this section.
1. Difficulties intrinsic in the current projection framework. LTI MOR techniques generally are formulated within a linear projection framework. However, it may not be
efficient/reasonable to use linear projection for nonlinear models, as it leads to less
computational speedup and lacks enough theoretical support, compared to LTI cases.
It is worthwhile to investigate how nonlinear projections may work better than linear projections, and that casts even more questions than answers regarding nonlinear
model order reduction.
2. Lack of canonical representation of a nonlinear model. LTI models (linear differential
equations) C dtd x + Gx + Bu = 0 have a very special and compact matrix representation (C, G, B). There are also other equivalent canonical representations such as
controllable/observable canonical forms [72, 73]. However, there have been hardly any
6
literature on canonical forms of general nonlinear systems. The nonlinear functions
can be arbitrary and complicated, making the system hard to analyze and manipulate,
and even harder to reduce.
3. Lack of characterization of system responses. LTI systems can be characterized by
their transfer functions which are extensively used as a performance metric to preserve
in linear model order reduction. However, for nonlinear systems, there is no simple
way to derive analytical expressions of system responses such as harmonic distortion,
output frequency, etc.. This further complicates the reduction problem.
1.3
In this dissertation, we formulate the model order reduction problem, explore and develop nonlinear model order reduction techniques, and apply these techniques in practical
examples. In particular, the main contributions are as follows:
1. We formulate the nonlinear model order reduction problem as an optimization problem
and present a general nonlinear projection framework that encompasses previous linear
projection-based techniques as well as the techniques developed in this dissertation.
2. We develop a nonlinear model order reduction technique, ManiMOR, which is a direct implementation of the nonlinear projection framework. It generates a nonlinear
reduced model by projection on a general-purpose nonlinear manifold. The proposed
manifold can be proven to capture important system dynamics such as DC and AC
responses. We develop numerical methods that alleviates the computational cost of
the reduced model which is otherwise too expensive to make the reduced order model
of any value compared to the full model.
3. We develop a nonlinear model order reduction technique, QLMOR, which transforms
the full model to a canonical QLDAE representation and performs Volterra analysis to
derive a reduced model. We develop an algorithm that can mechanically transform a
set of nonlinear differential equations to another set of equivalent nonlinear differential
equations that involve only quadratic terms of state variables, and therefore it avoids
any problem brought by previous Taylor-expansion-based methods. With the QLDAE
representation, we develop the corresponding model order reduction algorithm that
extends and generalizes previously-developed Volterra-based technique.
4. We develop a nonlinear phase/timing macromodeling technique, NTIM, that derives
a macromodel that specifically capture timing/phase responses of a nonlinear system.
We rigorously define the phase response for a non-autonomous system, and derive the
dynamics of the phase response. The macromodel emerges as a scalar, nonlinear timevarying differential equation that can be computed by performing Floquet analysis of
7
the full model. With the theory developed, we also present efficient numerical methods
to compute the macromodel.
5. We develop a finite state machine abstraction technique, DAE2FSM, that learns a
Mealy machine from a set of differential equations from its input-output trajectories.
The algorithm explores the state space in a smart way so that it can identify the underlying finite state machine using very few information about input-output trajectories.
The remainder of the dissertation is organized as follows. Chapter 2 formulates the MOR
problem, reviews previous work and summarizes challenges in nonlinear MOR. Chapter 3
provides a deeper view of the linear projection framework that is used by most previous
MOR methods. Chapter 4 generalizes the linear projection framework to a nonlinear one,
and summarizes the skeleton of an algorithm within this framework. Chapter 5 presents
the method ManiMOR that is a direct implementation of the nonlinear projection framework. Chapter 6 presents the method QLMOR that generates a quadratic-linear reduced
order model. Chapter 7 presents the method NTIM that generates a phase/timing reduced
model. Both QLMOR and NTIM can be interpreted in the nonlinear projection framework.
Chapter 8 looks at a slightly different problem of deriving a finite state machine model from
a set of differential equations, and presents the method DAE2FSM. Chapter 9 concludes the
dissertation and discusses the related future work.
1.4
Except for special mention, we will use use the following notations and abbreviations
throughout the dissertation.
f ()
Jacobian of function f
r
residual of the differential equations
ei
the column vector [0, , 0, 1, 0, , 0]
with the i-th element being 1 and other elements being 0
transfer function of an LTI system
H(s)
span(A)
column space of matrix A
span(v1 , , vN ) span of vectors v1 , , vN
Kronecker product
M N
Chapter 2
Background and Previous Work
In this chapter, we formally formulate the model order reduction problem, trying to
encompass all previous methods under the same formulation. We review previous linear and
nonlinear model order reduction techniques, and summarize a few challenges in nonlinear
MOR, some of which are addressed in this dissertation.
2.1
Problem Formulation
Suppose that we are given a full model M M, where M is the set of a specific type of
model, (e.g., ordinary differential equations). The goal of model order reduction is to derive
a reduced model Mr Mr , where Mr is the set of another type of models (e.g., ordinary
differential equations), so that Mr approximates M in some sense.
Mathematically, we can define the problem of computing Mr from M as the following
optimization problem
minimize F (M, Mr )
Mr
subject to G(M, Mr ) = 0,
H(M, Mr ) 0,
(2.1)
where F (, ), G(, ), H(, ) are functions on models M and Mr . F (, ) specifies how close
two models are. G(, ) and H(, ) specify the equality and inequality constraints.
Choices of M, Mr , F (, ), G(, ) and H(, ) lead to different problems. Usually, if Mr
has less number of state variables (i.e., the order of the system) than M, it is called a reduced
order model, and Mr is expected to be computationally more efficient than M.
2.1.1
Types of Models
As listed in Table 1.1, there are many types of mathematical models. However, in the
MOR problem, not all pairs of models (M, Mr ) are reasonable. For example, it may be
10
unreasonable to reduce an LTI model to a nonlinear model because they may exhibit fundamentally different behaviors; it may also be unreasonable to convert a deterministic model
to a non-deterministic one and vise versa.
In this dissertation, unless otherwise mentioned, we focus on the problem where M and
Mr are both linear ODEs/DAEs or nonlinear ODEs/DAEs.
We list possible model reduction problems with other choices of M and Mr in Table 2.1.
Table 2.1: Choices of mathematical models in MOR
Original Model M Reduced Model Mr
Applications
ODEs/DAEs
ODEs/DAEs
Traditional MOR
ODEs/DAEs
FSMs
FSM abstraction/learning
ODEs/DAEs
Hybrid Systems
Hybridization
Hybrid Systems
FSMs
FSM abstraction/learning
ODEs/DAEs
Difference Equations
Time Discretization
PDE
ODE/DAE
Time/Space Discretization
References
[7478]
[7880]
[79, 81, 82]
[83, 84]
[85]
[85]
Remark In Table 2.1, the mappings M Mr can be one-to-one, one-to-many, many-toone or many-to-many. For example, for continuous ODEs, we can apply different integration/discretization schemes to obtain different difference equation models. Another example,
widely used in IC design, is to map a finite state machine to a circuit implementation which
corresponds to a set of ODEs/DAEs. With different circuit topologies, we obtain different
ODEs/DAEs that realize the same FSM functionality. It is sometimes also interesting to
look at the reverse problem to map Mr to M. For example, the mapping from FSMs to
ODEs may correspond to the synthesis problem [8688].
2.1.2
Model Fidelity
In (2.1), F (, ) defines a measure describing how close two models are, or how well Mr
approximates M.
For an input-output system whose inputs are u and outputs are y, a straightforward
choice of F (, ) is
X
F (M, Mr ) ,
||y yr ||,
(2.2)
uU
where U is a given set of input trajectories/traces; y and yr are outputs of M and Mr under
inputs u, respectively. In general, y and u may be described in different forms, depending
on applications.
For example, since LTI systems are typically described by their transfer functions, y can
be chosen as the transfer function and u can vary over a range of frequencies. For nonlinear
systems, often u is chosen as a set of (training) input signals, and y is chosen to be the
11
corresponding time-domain responses, or the first few harmonic distortions when the system
is settled at the periodic steady state. In specific applications where other special responses
or higher-level performance metrics (such as phase/timing responses or bit error rate) are of
interest, y can also be set to these quantities.
There have also been many techniques in reducing dynamical systems with no inputs
(e.g., the ODE model x = f (x)). In this case, we may still view it as an input-output
system, with the virtual input to the system being the set of all possible impulse functions
that set the initial condition of the system when applied. These methods basically compute
reduced order models that approximate outputs of full models for arbitrary initial conditions.
G(, ) and H(, ) in (2.1) impose equality and inequality constraints on the reduced
order model, respectively. They are usually used to guarantee certain important behaviors
of reduced models.
For example, in moment-matching methods [89] for LTI systems, G(, ) is defined to be
mi (M) = mr,i (Mr ),
(2.3)
where mi (M) is the i-th moment of the transfer function of model M (see Section 2.3.2 for
details). This guarantees that the first few moments of the reduced model mr,i match those
of the full model mi .
In stability-guaranteed reduced-order modeling techniques [90,91], H(, ) may be defined
as a linear matrix inequality
Pr Ar + ATr Pr 0
(2.4)
where Ar is the matrix in the reduced model dtd x = Ar x and Pr is a positive definite matrix.
To further guarantee the passivity of the reduced model, we may simply expand G(, ) by
one more equality constraint [91].
While there are many fidelity metrics that prove to be useful in different applications, we
try to summarize the most important ones that reduced models should satisfy, as follows:
1. Attractor preservation. Attractors, such as equilibria (DC states), limit cycles, or even
strange attractors, are system responses as time t goes to infinity. Many engineering
systems are designed to work in a region close to its DC states or limit cycles, and
therefore the attractors are the most important system response to preserve in the
reduced model.
Nonlinear systems present much richer attractors than linear systems, such as multiple
DC equilibria and limit cycles. They correspond to important circuits such as memory
cells and oscillators.
2. Stability preservation. Stability describes how the state of a dynamical system moves
under small perturbation at the input. It is a key system property of interest, and
has many variants of definitions [73, 90, 91]. The stability (e.g., of equilibria) should
12
be preserved in the reduced model since it ensures that the qualitative behavior of the
system is correct, given that the attractors are preserved in the model.
3. Passivity preservation. Passivity of a system refers to the property that the system
only consumes energy and never produces energy. This is commonly seen in practical
systems, such as RLC interconnects. This is another qualitative behavior of a system,
and is desirably preserved in reduced models.
4. Structure preservation. Reduced models are often written as a set of dense (in terms
of Jacobian matrices) and less structured differential equations. However, full models
usually have certain structures. For example, the system matrices may be sparse,
symmetric or have block patterns; the system may be composed of a cascade of several
sub-systems, or a feedback; the system may be be a second order system; the system
may obey certain conservation laws; the system may correspond to a realistic circuit
implementation. Therefore, it is sometimes desirable for reduced models to preserve
these structures, even with certain accuracy loss.
5. Linearity preservation. This refers to the property that the reduced model of a linear
system should remain linear. This is rather a necessary condition that should be
appreciated in practice, especially for nonlinear MOR techniques. Linear systems have
limited behaviors than nonlinear systems, and should not exhibit any nonlinear effects
induced by model reduction methods.
Example 2.1.1 The TPWL [2] reduced model, reviewed in Section 2.4.3, for the ODE
model (1.2) is
X
d
xr =
wi (xr )(V T f (xi ) + V T G(xi )V (xr xr,i )).
dt
i
(2.5)
When f (x) = Gx, i.e., the full model is linear, then (2.5) becomes
X
d
xr =
wi (xr )(V T Gxi + V T GV (xr xr,i )).
dt
i
(2.6)
Because V xr,i = xi does not always hold, the TPWL model (2.6) of an LTI model
becomes nonlinear.
6. Composability. One application of reduced models is to replace complex system blocks
by their reduced models so that simulation/validation of larger systems can be computationally more efficient. When the full models are replaced by reduced models and
these reduced models are composed together, we hope to retain the coupling effects,
and important behaviors (e.g., stability) of the composed system. Composability
13
also allows one to generate reduced models of sub-systems separately and obtain reduced models of larger sub-systems by composition. It can lead to a hierarchical model
order reduction framework.
7. Identification of simpler models under coordinate transformation. A nonlinear system,
in some cases, can be converted to a simpler model (e.g., with less nonlinearity), or
even to an LTI model. In such cases, the nonlinear MOR problem may be reduced to
a linear MOR problem, and therefore becomes much easier. One such example is as
follows:
Example 2.1.2 (True Linearization of Nonlinear Models) Consider the nonlinear system
x1 = x1 + u,
x2
x2 = u,
x1
(2.7)
(2.8)
This simple example shows the possibility of finding appropriate coordinate transformations to simplify nonlinear models.
Such capability is desirable in a nonlinear MOR technique.
8. Separation of independent/repetitive sub-systems. There are many cases where the
full model is composed of several independent sub-systems. Some sub-systems may
be redundant in terms of controllability, observability, repetition, etc.. It is therefore
another desirable property for MOR techniques to separate these sub-systems automatically.
Example 2.1.3 Fig. 2.1, Fig. 2.2 and Fig. 2.3 show three special examples where the
full model is composed of two sub-systems.
(a) Shown in Fig. 2.1, the output y depends solely on the sub-system A. Therefore subsystem B is redundant from an input-output perspective, and should be identified
and removed in model order reduction.
14
15
16
reduced due to sparsity. Also, a sparse projection may lead to reduced models whose
state variables are associated with physical quantities.
Summarizing the above discussion, in Table 2.21 , we give a relatively complete list of the
model fidelity that has been, or potentially can be, used in model order reduction.
Table 2.2: Model fidelity and examples
Model Fidelity
Example in terms of F (, ), G(, ) and H(, )
Moment matching
G(M, Mr ) = m
Pi (M) mr,i (M) = 0
Minimum output-MSE
F (M, Mr ) = uU ||y(t) yr (t)||2
G(M, Mr ) = limt (y(t) yr (t)) = 0
Attractor preservation
Stability preservation
H(M, Mr ) = Pr Ar + ATr Pr 0
Passivity preservation
G(M, Mr ) = Pr B D = 0 with stability condition
Structure preservation
G(M, Mr ): Mr corresponds to a realizable RLC circuit.
Linearity preservation
G(M, Mr ): Mr is linear, given M is linear.
Composability
G(M, Mr ) and/or H(M, Mr ):
suppose M = (MA , MB ), Mr = (Mr,A , Mr,B ),
then Mr preserves certain properties of M.
Redundancy removal
Mr removes non-observable and non-controllable sub-systems.
Lower computational cost H(M, Mr ): Computational cost of Mr is less than that of M.
Sparsity preservation
F (M, Mr ): number of non-zero elements in Jacobians in Mr .
In Table 2.3 in Section 2.5, we also summarize a selected set of previous model order
reduction techniques in terms of the optimization problem (2.1), i.e., in terms of the functions
F (, ), G(, ) and H(, ).
2.1.3
Model Verification
In (2.1), if F (, ) is null, the problem is called the feasibility problem in literature. It can
be viewed as a verification problem to check whether two models M and Mr are equivalent
(under the constraints imposed by G(, ) and H(, )).
The solution of this problem renders a model Mr that satisfies all the constraints G(, )
and H(, ). If the (perfect) algorithm cannot find such a model, then such an equivalent
model Mr Mr does not exist.
This problem has much wider engineering applications. For example, in the digital circuit
synthesis and verification flow, abstraction/models of different levels are verified against
each other (equivalence checking) to make sure that the abstraction/synthesized circuit is
1
There are many choices of constraints that may lead to the same fidelity criteria, but we only provide
one example for each to illustrate the idea. The symbols used in the table are straightforward and therefore
without detailed explanation.
17
functionally correct. Solving such problems is therefore the key enabler of the digital circuit
design flow.
However, in model order reduction literature, this equivalence checking procedure is
typically skipped, or is done by heuristics. We believe that this is an important problem to
solve for employing MOR in practice, and some initial studies [79, 9296] in analog/mixedsignal verification may be used.
2.2
MOR Overview
According to the description of the full system, model order reduction methods can be
classified into black-box methods and white-box methods.
In black-box methods, the internal structure or representation of the model is not accessible. Rather, we have restricted access to only a few controllable inputs (e.g., stimulations)
and observable outputs (e.g., from measurements or simulations). This leads to various
fitting algorithms, including regression, neural network, machine learning, etc.. These
methods are often also called system identification methods. Due to the nature of such
methods, we normally cannot assert any deterministic relationship between the full model
and the reduced model beyond the training data.
In white-box methods, we have full access to the model. This gives us a complete view
of the full model, and thus opportunities to guarantee certain properties can be preserved
from full models. Depending on the type of the full model, the behavior can be qualitatively
different, and therefore different methods have been proposed to handle different types of
models, including LTI, LTV and nonlinear models. The most common framework encompassing most of these methods is the linear projection framework, to be introduced in this
section, and detailed in Chapter 3.
2.2.1
System Identification
System identification methods [97,98] view the system as a black/grey box. Usually only
limited number of input-output data/trajectory is given, either by simulation or measurement. Sometimes, partial information of the system is known (such as the connectivity of
a graph). These methods can be classified into two categories: parametric methods and
non-parametric methods.
In parametric methods, we make assumptions of the model structure for the system, and
parameterize the model by certain key parameters. With the parameterized model, the data
are used to fit the model by solving for parameters that minimizes the error between the fitted
model and the data. For example, we may use methods such as least squares, maximum
likelihood, maximum a posteriori probability estimation, compressed sensing, etc..
In non-parametric methods, we construct the model directly from data. For example, to
find an LTI model of a system, we can measure directly the impulse response or step response
18
to reconstruct the LTI system. To compute the transfer function at specific frequency points,
we can measure directly the response to a sinusoidal input. Often the experiments are
carefully designed for the identification purpose in non-parametric methods.
2.2.2
(2.9)
(2.10)
by projection.
That is, the nonlinear functions in (2.10) are defined by projection, i.e.,
qr (xr ) = W T q(V xr ),
fr (xr ) = W T f (V xr ),
(2.11)
br (t) = W T b(t),
where V, W RN Nr are two projection matrices.
Therefore, the problem of MOR amounts to finding two appropriate projection matrices.
2.3
For an LTI system, q(x) = Cx and f (x) = Gx are linear functions of x, and often the
input is written in the form of b(t) = Bu(t), i.e.
C
d
x + Gx + Bu = 0,
dt
y = D T x,
(2.12)
d
xr + Gr xr + Br u = 0,
dt
y = DrT xr ,
(2.13)
19
reduced model.
2.3.1
Since LTI models can be equivalently defined by their transfer functions H(s), a straightforward model order reduction method is to fit a reduced transfer function Hr (s) that gives
a good approximation to H(s).
The famous approach is the Pad`e approximation [99], i.e., to find a rational reduced
transfer function of the form
Hr (s) =
a0 + a1 s + a2 s2 + + am sm
b0 + b1 s + b2 s2 + + bn sn
(2.14)
which satisfies
Hr (0) = H(0),
k
d
dk
H
(0)
=
H(0),
r
dsk
dsk
k = 1, 2, , m + n,
(2.15)
d
where ds
k H(0) is also referred to as the k-th moment of the transfer function. Therefore, any
LTI MOR method that leads to a reduced model satisfying (2.15) is also called a momentmatching method.
Intuitively, Pad`e approximation is good for MOR because a rational function is a good
fit for LTI transfer functions the original transfer function H(s) is typically a rational
function given that there is no repeating poles. In fact, we may analytically write out H(s)
of (2.12) as
H(s) = D T (sC + G)1 B,
(2.16)
2.3.2
Krylov-Subspace Methods
The explicit moment calculation in the Pad`e approximation methods can have serious
numerical inaccuracy, and this inaccuracy will be inherited to the reduced model. Even with
PVL, since it uses a transfer function representation of the LTI system, it is not directly
clear what the state-space ODE/DAE model of the reduced model is.
20
The source of error in explicit moment computation comes from the fact that higherorder moments tend to decay/grow at an exponential rate. To see that, we may derive the
expression of the i-th moment to be
mi = D T (G1 C)i B.
(2.17)
span(W ) = KNr ((G1 C)T , D) = span((G1 C)T D, ((G1C)T )2 B, , ((G1C)T )Nr D),
(2.18)
then the projected reduced models with Cr = W T CV , Gr = W T GV , Br = W T B, Dr = V T D
have the property that the first 2Nr moments of the reduced transfer function match those
of the full transfer function.
Moreover, the projection matrices V and W can be efficiently and accurately computed
(e.g., Arnoldi and Lanczos algorithms [102, 103]). These algorithms are computationally
efficient and numerically stable even for large scale problems (e.g., sparse matrices of dimension on the order of millions), and therefore enables the moment-matching technique to be
applicable to large-scale problems.
2.3.3
The method of truncated balanced realization (TBR) originated from control theory,
based on the idea of controllability and observability [104]. Loosely speaking, certain states
in a dynamical system are hard to control and some others are hard to observe. The balanced
realization obtains a dynamical system whose state variables have equal controllability and
observability. Then, the truncation of the states that are hard to control and observe leads
to a reduced order model.
It turns out that the TBR method also falls in the linear projection framework. Ultimately, it defines two matrices V and W , and the reduced model is a projected model.
2.3.4
Proper orthogonal decomposition (POD) is also a method that derives reduced models
by linear projection. It takes several data/trajectories of the state variables, and constructs
a linear subspace that minimizes certain error between the projected reduced model and the
21
full model.
On the one hand, POD is general because it can take any trajectory of state variables.
Therefore, prior information of input signals and system dynamics can help generating a
better reduced order model. In can also be shown that if the trajectory is chosen to be the
impulse response of an LTI system, then one version of the POD is equivalent to the TBR
method.
On the other hand, POD is less systematic than other methods in the sense that there
is no universal guideline in choosing a good trajectory. Therefore, it is possible to have
bad trajectories that lead to bad reduced order models.
2.4
Nonlinear MOR methods are mostly built upon the LTI MOR methods. For so-called
weakly-nonlinear systems, we may approximate the nonlinear model by a local model, such as
LTI, LTV, Volterra models. The reduction is then performed on the approximated model to
obtain an approximated reduced model. For general nonlinear systems, we may approximate
the nonlinear model around several different regions in the state space, and reduce each
linearized model using LTI MOR methods.
2.4.1
There are many practical systems that may be well-approximated by linear time-varying
systems, partly because these systems are engineered to work under time-varying operating
points. For example, in RF circuits such as mixers, the local frequency input is fixed,
rendering an approximate periodic time-varying system. Similarly, in switched-capacitor
circuits, the clock input is fixed, and the circuit is approximately working periodically under
two modes (clock high and clock low).
The MOR methods in [25, 28] extends the notion of transfer function for LTI system to
time-varying transfer functions. Then it can be shown that using a Pad`e-like approximation, or an appropriate projection, we can obtain reduced order models whose time-varying
transfer functions approximate the original one.
2.4.2
Volterra-Based Methods
Generally, for a nonlinear system, we may derive a Volterra series approximation [105]
which approximates its local behavior. Therefore, instead of matching the moments of
transfer functions, we may match the moments of Volterra kernels (higher-order transfer
functions), and this leads to the Volterra-based MOR methods such as NORM [106] and [107].
Since the Volterra series only converges when input is small, Volterra-based models are
small-signal or local models.
22
Similar to Krylov subspace methods, in this case, we may also choose projection matrices
V and W that cover a series of Krylov subspaces so that the projected models match the
specified moments of Volterra kernels.
2.4.3
Trajectory-Based Methods
Unlike Volterra-based methods, trajectory-based MOR methods generate global reduced models, instead of local ones.
Among several variants of trajectory-based methods, the representative one is the TPWL
method [2] whose algorithm is shown in Algorithm 1. The TPWL method identifies two linear
subspaces V and W , and generates the reduced model by linear projection. The subspaces
V and W are the aggregated subspaces of many reduced models of linearized models (e.g.,
Krylov subspaces of linearized models), therefore making sure that the moments of transfer
functions of reduced linearized models match those of the full linearized models. In order to
achieve the computational efficiency of the reduced model, i.e., the computation of the term
W T f (V z), it uses a weighted-summation of linear functions to approximate W T f (V z).
Algorithm 1 TPWL (Trajectory Piece-Wise Linear) Model Order Reduction
Inputs: Differential equations of the full model dtd q(x) + f (x) + Bu = 0.
Outputs: Two linear subspaces defined by the column span of V and W . Differential equations of the reduced model dtd qr (xr ) + fr (xr ) + Br u = 0.
where
PNs
P Ns
T
T
T
qr (xr ) = i=1 wi (xr )(W q(xi ) + W C(xi )V (xr xr,i )), fr (xr ) = i=1 wi (xr )(W f (xi ) +
W T G(xi )V (xr xr,i )), Br = W T B, wi (xr ) are weighting functions described in Section 2.4.3,
Ns is the number of sample points.
1: Simulate the full model using a set of training input signals, and obtain a set of
trajectories.
2: Sample Ns points {xi }, i = 1, , Ns on these trajectories as expansion points of the
reduced model.
3: for all i = 1, , Ns do
4:
Linearize the full model around each sample point xi
5:
Compute the reduced linearized models, and obtain projection matrices Vi and Wi for
each reduced model.
6: end for
7: Compute the projection matrix for the nonlinear model V = [V1 , x1 , , Vk , xk ], W =
[W1 , x1 , , Wk , xk ].
8: Compute W T f (xi ), W T G(xi )V , W T q(xi ), W T C(xi )V and W T B, and store them in the
reduced model.
23
2.5
We summarize a selected set of methods in Table 2.3, relating them with our formulation
(2.1).
Table 2.3: Summary of previous MOR methods in terms of the optimization problem
Methods
F (, )
G(, )
H(, )
AWE [100]
null
moment-matching mi (M) = mr,i (Mr )
null
PVL [101]
null
moment-matching mi (M) = mr,i (Mr )
null
PRIMA [74]
null
moment-matching mi (M) = mr,i (Mr )
null
the full model is the MNA equation and is passive
TBR [104]
null
controllability and observability are balanced
null
NORM [106] null
moment-matching of Volterra kernels
null
TPWL [2]
null
moment-matching of linearized models
null
2.6
Rethinking about the nonlinear MOR problem, we summarize the following major challenges:
1. Lack of canonical form for the reduced model. LTI systems can be completely represented by several matrices C, G, B, D. However, nonlinear systems generally do not
have a universal canonical form, making it hard to efficiently represent a nonlinear
system as well as to develop algorithms manipulating the full model.
2. Lack of analytical expressions for system responses. Unlike LTI systems which are
characterized by transfer functions, nonlinear systems generally do not have such explicit analytical formula for system responses. Volterra kernels of polynomial systems
alleviate this problem.
3. Limited guarantee of reduced order model behaviors. Often we are only able to quantify the behavior of the reduced model locally or with respect to given training data.
However, global properties such as stability and multiple equilibria are extremely important for nonlinear systems.
4. Potentially expensive computational cost. In linear projection-based methods, the
nonlinear function in the reduced model W T f (V xr ) can be as expensive to compute as
the full model, thus making the reduced model of less use compared to the full model.
24
2.7
Benchmark Examples
We have created and collected a few benchmark systems for testing MOR methods. Some
of them are copied from existing literatures. A summary of these benchmarks is shown in
Table 2.4. Detailed description of these benchmarks are presented in Appendix Chapter B.
Benchmark
NTL1
NTL2
INVC
NEURON ML
NEURON FN
LATCH1
Size
variable
variable
variable
2
variable
3
Most of these examples are from circuit applications. Some of them describe the chemical
reactions, biological pathways, or neuron dynamics. We briefly review the circuit MNA
equations and chemical rate equations.
2.7.1
In circuit simulation, the MNA equation formulation [108, 109] is most widely used.
Briefly, in the MNA equations, the state variables (unknowns) are mainly the node voltages, together with a few branch currents of non-voltage-controlled devices. The differential
equations mainly consist of KCL equations, i.e., each equation defines the sum of currents
flowing into one node to be 0. For non-voltage-controlled devices, the BCR equations are
added to describe the voltage-current relationship of these devices.
2.7.2
Chemical kinetics can be modeled by a set of differential equations, called rate equations.
They describe the dynamics/rate-of-change of concentrations of reactants. They are usually
in the form of ODEs dtd x = f (x, u), where x are the concentrations of reactants, and u are
the controllable parameters such as reaction rates, temperatures.
K
There are rate laws that determine the function f (x, u). For a generic reaction A+B
C, the f function has expressions such as K[A]m [B]n , where K is a constant, and [A], [B]
are the concentrations of A and B, respectively. Some other rate laws include models of
Henri, Michaelis & Menten, and Briggs Haldane [110]. For example, in Michaelis & Menten
x
model, the f function has expressions such as k+x
.
25
Chapter 3
Linear Projection Framework for
Model Order Reduction
While many previous MOR methods have used the linear projection framework, the major
motivation seems to be the fact that the projection of the full model on certain subspaces
leads to a reduced model that matches certain accuracy metrics (e.g., moments or Hankel
singular values). However, no literature has discussed why such such a framework is a good
one for MOR. For example, some key questions that need to be answered are: can any
reduced model be written as a projected model; what is the degree of freedom in choosing
projection matrices; what might be other choices of projections; etc..
In this chapter, we take a deeper look at the linear projection framework. We discuss
the application of the linear projection framework in LTI MOR the representation power
of the linear projection framework, the degree of freedom in choosing projection matrices,
the other potential criteria in choosing projection matrices. We also discuss the application
of the linear projection framework in nonlinear MOR, and discuss potential problems and
solutions.
3.1
The basic idea of the linear projection work has already been presented in Section 2.2.2.
Here, we focus on the algorithmic aspect of the linear projection framework. The problem
is to reduce a model M,
d
q(x) + f (x) + b(t) = 0,
(3.1)
dt
to a reduced model Mr ,
d
qr (xr ) + fr (xr ) + br (t) = 0.
(3.2)
dt
26
The algorithm is generally in the form of Algorithm 2. Obviously, the choice of V and W
plays a key role in the linear projection framework.
Algorithm 2 MOR based on Linear Projection
Inputs: dtd q(x) + f (x) + b(t) = 0
Outputs: dtd qr (xr ) + fr (xr ) + br (t) = 0
1: Construct a matrix V RN Nr ;
2: Construct a matrix W RN Nr ;
3: Define qr (xr ) = W T q(V xr ), fr (xr ) = W T f (V xr ), br (t) = W T b(t).
3.2
For simplicity, we consider the LTI system where q(x) = x, f (x) = Gx, b(t) = Bu, i.e.,
d
x + Gx + Bu = 0,
dt
y = D T x.
(3.3)
Therefore, under the linear projection framework, we obtain a reduced model where
qr (xr ) = W T V xr , fr (xr ) = W T GV xr , br (t) = W T Bu.
To further simplify the analysis and discussion, we assume that the left projection matrix
W is chosen as
W T = (V T V )1 V T , (W T V = I),
(3.4)
and therefore, the reduced model is
d
xr + Gr xr + Br u = 0,
dt
y = DrT xr
(3.5)
where
Gr = W T GV,
Br = W T B,
(3.6)
Dr = V T D.
3.2.1
Under the optimization formulation (2.1), the LTI reduced model is in the set Mr,LT I of
all differential equations in the form of (3.5), and is determined by all choices of (Gr , Br , Dr ),
27
i.e.,
d
T
Nr Nr
Nr Ni
Nr No
Mr,LT I =
xr + Gr xr + Br u = 0, y = Dr xr Gr R
, Br R
, Dr R
.
dt
(3.7)
In contrast, the reduced model in the set of projected models Mr,V , (3.5) with (3.6), is
determined by all pairs of (V, W ) where W T = (V T V )1 V T , i.e.,
d
N Nr
T
T
1 T
T
T
T
.
, W = (V V ) V
xr + W GV xr + W Bu = 0, y = D V xr V R
Mr,V =
dt
(3.8)
Furthermore, instead of choosing the projection matrix V directly, we may choose a
Nr -dimensional subspace V in RN , i.e., an element in the Grassmann manifold [111]. Mathematically, V Grass(Nr , N) = RN Nr /GLNr where GLNr is the general linear group,
i.e., the set of all Nr Nr invertible matrices. The dimension of the Grassmann manifold
Grass(Nr , N) is therefore Nr (N Nr ). Using V, the set of projected models is
d
T
T
1
T
T
T
T
Mr,V =
.
xr + W GV xr + W Bu = 0, y = D V xr V Grass(Nr , N), W = (V V ) V
dt
(3.9)
From (3.7), (3.8), (3.9), we may conclude that Mr,V is a subset of Mr,LT I , and that
there is a many-to-one mapping from Mr,V to Mr,V , as illustrated in Fig. 3.1. It is also
tempting to say that the degrees of freedom in choosing a model in Mr,LT I , Mr,V , Mr,V are
Nr2 + Nr (Ni + No ), Nr N, Nr (N Nr ), respectively, by counting the number of variables.
28
2. Are the projected models good candidates to use in MOR (i.e., do they encompass all
possible LTI reduced models)? If not, why are the projected models useful?
3.2.2
To determine the degree of freedom in choosing an LTI model in the form of (3.3),
we study its canonical forms. Mathematically, given a set and an equivalent relation, the
canonical form is a unique representation of an equivalence class, i.e., all equivalent elements
can be mapped to the same canonical form. Therefore, by studying the canonical forms, we
may determine its degree of freedom.
For simplicity, we consider SISO systems in the following, i.e., Ni = No = 1. For SISO
LTI systems, the canonical forms include transfer function, controllable canonical form,
observable canonical form, etc.. In particular, we study the transfer function and controllable
canonical form.
Definition 3.2.1 (Transfer function) The transfer function of (3.3) is
H(s) = D T (sI + G)1 B
(Let G = P P 1 be its eigen-decomposition) = D T P (sI + )1P 1 B
(Assuming is diagonal) =
N
X
i=1
ri
,
s + i
(3.10)
where
ri = dibTi R,
= [d1 , , dN ] = P T D,
D
= [bT ; ; bT ] = P 1B.
B
1
(3.11)
G B=
n1
X
i Gi B.
(3.12)
i=0
Then we can find a matrix T RN N so that by defining x = T , the LTI model in terms
of is in the controllable canonical form, i.e.,
d
= T 1 GT + T 1 Bu,
dt
y = D T T ,
(3.13)
29
where
T 1 GT =
0
0
..
.
1
0
..
.
0
1
..
.
..
.
0
0
..
.
0
0
..
.
0
0
0
1
0
0
0
0
0
1
0 1 2 N 2 N 1
0
0
..
.
T 1 B = .
0
0
1
(3.14)
= [GB, G2 B, , GN B]
0 0 0 0 0
1 0 0 0 1
..
0 1 0 0
.
= Mc .. .. . . . . ..
..
.
. .
.
. .
0 0 0 . . . 0 N 2
0 0 0 1 N 1
= Mc Q.
(3.15)
= T 1 GT and B
= T 1 B. Because the two systems have the same characteristic
Let G
polynomial, we have
M
c = M
c Q.
G
(3.16)
Consequently, choosing
c1 ,
T = Mc M
(3.17)
leads to (3.14).
To fully specify (3.13), we need to choose = [0 , , N 1 ] RN , and D T T RN ,
leading to a degree of freedom 2N.
Therefore, we conclude the following result
Theorem 3.2.2 The degree of freedom in determining an N-dimensional SISO LTI model
is 2N.
From theorem 3.2.2, the degree of freedom in determining an Nr -dimensional reduced
model is 2Nr . We also see that the degrees of freedom in choosing Mr,LT I , Mr,V , Mr,V are
misleading because in each of these sets, a lot of elements are equivalent.
30
Corollary 3.2.3 (Modal form) Let G = P P 1 where diag() are the eigenvalues of G
and columns of P are eigenvectors of G. Define x = P 1x, and replace x = P x in (3.3), we
obtain the decoupled form
d
x =
x + P 1 Bu,
dt
y = D T P x,
(3.18)
y = D1T x
(3.19)
and
d
x + G2 x + B2 u, y = D2T x
dt
are equivalent if they have the same transfer function, i.e.,
H1 (s) = D1T (sI + G1 )1 B1 = H2 (s) = D2T (sI + G2 )1 B2 .
(3.20)
(3.21)
31
3.2.3
y = D T V1 xr ,
(3.22)
d
xr + W2T GV2 xr + W2T Bu = 0,
dt
y = D T V2 xr ,
(3.23)
and
We derive conditions for (3.22) to be equivalent to (3.23), and then discuss the degree of
freedom in choosing a projected model.
Theorem 3.2.6 If span(V1 ) = span(V2 ), then (3.22) is equivalent to (3.23).
Proof If span(V1 ) and span(V2 ) are the same, then there exists an invertible matrix Q
RNr Nr , such that V2 = V1 Q.
Therefore, we have
H2 (s) = D T V1 Q(sI + (QT V1T V1 Q)1 QT V1T GV1 Q)1 (QT V1T V1 Q)1 QT V1T B
= D T V1 Q(sQT V1T V1 Q + QT V1T GV1 Q)1 QT V1T B
= D T V1 (sV1T V1 + V1T GV1 )1 V1T B
= H1 (s).
(3.24)
32
However, the projected models have limited power in representing general Nr -dimensional
models. Therefore, we have the following theorem.
Theorem 3.2.7 Projected models cannot represent arbitrary Nr -dimensional LTI SISO models. In particular, they are highly dependent on the original full model.
In the following, we show several counter examples in which projected models are not
general enough to cover all possible Nr -dimensional models:
Example 3.2.8
1. Suppose N = Nr + 1, we have only Nr degrees of freedom in choosing
the projected model, while we have 2Nr degrees of freedom in choosing an arbitrary
Nr -dimensional LTI SISO model. It is therefore not possible to cover all possible LTI
SISO models.
2. Consider deriving a one-dimensional projected LTI model by choosing V RN 1 and
V T V = 1. Let G be positive definite. Therefore, in the reduced model,
Gr = W T GV = V T GV 0.
(3.25)
Therefore, the reduced model is always stable. While this is probably rather a good result
in practice, it shows that the projected model cannot represent any non-stable model.
3. If G is identity, we have
G r = W T I N V = I Nr ,
(3.26)
and therefore, the poles of the reduced transfer function can only be 1. The projected
model is obviously not general.
In general, in order to be equivalent to a given Nr -dimensional LTI model, the projected
model must satisfy 2Nr or Nr2 + Nr constraints for SISO or SIAO systems, respectively, (i.e.,
by equating variables in the canonical forms). Therefore, the number of elements in the set of
Nr -dimensional subspaces must be larger than the degree of freedom (i.e., Nr (N Nr ) 2Nr
(N Nr + 2) or Nr (N Nr ) Nr2 + Nr (N 2Nr + 1), respectively), in order to make
the projected model general. While there is still no guarantee that the projected model is
general enough, this appears as a viable necessary condition.
While projected models are not general enough to represent all Nr -dimensional models,
they are still most widely used, due to the following reasons:
1. Projection is a natural way to map from G RN N to Gr RNr Nr ;
2. There are efficient algorithms to generate projected reduced models that are reasonably
accurate.
33
3.2.4
3.3
y = D T x,
(3.27)
y = DrT xr ,
(3.28)
34
where
fr (xr ) = W T f (V xr ),
Br = W T B,
(3.29)
Dr = V T D.
3.3.1
Similar the the linear case, the nonlinear reduced model is in the set Mr,N L of all differential equations in the form of (3.28), and is determined by all choices of (fr (), Br , Dr ),
i.e.,
d
Nr
Nr
Nr Ni
Nr No
T
Mr,N L =
, Dr R
xr + fr (xr ) + Br u = 0, y = Dr xr fr : R R , Br R
.
dt
(3.30)
In contrast, the reduced model in the set of projected models Mr,V , (3.28) with (3.29),
is determined by all pairs of (V, W ) where W T = (V T V )1 V T , i.e.,
d
T
T
T
N
N
T
T
1
T
r
.
xr + W f (V xr ) + W Bu = 0, y = D V xr V R
, W = (V V ) V
Mr,V =
dt
(3.31)
Furthermore, instead of choosing a projection matrix V , but choosing a Nr -dimensional
subspace V in RN , we have
d
T
T
1
T
T
T
T
.
Mr,V =
xr + W f (V xr ) + W Bu = 0, y = D V xr V Grass(Nr , N), W = (V V ) V
dt
(3.32)
We hope to answer the same two questions we have addressed in the linear case:
1. What is the degree of freedom in choosing a nonlinear reduced model?
2. Are the projected models good candidates to use in MOR (i.e., do they encompass all
possible nonlinear reduced models; if not what subset of models can be represented)?
However, we are not able to give exact answers, due to the lack of canonical form and
the lack of analytical response functions such as transfer functions.
3.3.2
In the full model, the nonlinear function f (x) is sparse, meaning that each element
fi (x) depends on only a few variables (K of them), and only takes at maximum L flops for
each nonlinear kernel. Therefore, the function evaluation takes NKL flops, and the Jacobian
matrix has only NK non-zero elements and is sparse.
35
In contrast, the projected nonlinear function fr (xr ) is dense, meaning that fr,i (xr )
depends on all variables in xr . Without any symbolic simplification, fr (xr ) can be computed
by computing x = V xr , f (x) and W T f (x) sequentially, resulting in NKL + 2NNr flops.
This is always larger than the that of the full model.
On the other hand, unlike the LTI case where W T GV can be pre-computed and stored in
the reduced model, there is no simple way to store a nonlinear function W T f (V ), therefore
resulting in difficulties of storing the reduced model.
We will discuss methods to handle such difficulties in Section 4.7.
36
Chapter 4
Nonlinear Projection Framework for
Model Order Reduction
While the linear projection framework works very well for reducing LTI models, it is less
so for reducing nonlinear models. In this chapter, we study how linear projections can fail
in nonlinear MOR, from which we motivate and propose the nonlinear projection framework
that provides an elegant alternative for nonlinear MOR. We analyze each component in
this framework (e.g., choosing manifolds and projections), and present general guidelines of
deriving nonlinear projected reduced models under this framework.
4.1
Motivation
The linear projection framework amounts to finding two subspaces V = span(V ) and
W = span(W ) and imposing constraints x V, r W. Loosely speaking, we want to find
two subspace V and W so that the solution of the full model stay very close to V and that
the residual of the reduced model solution is orthogonal to W.
Example 4.1.1 (A simple LTI system) Consider a simple LTI system defined by
x1
10 1
1
x1
1
d
x2 =
1 1 0
x2 + 0 u(t),
(4.1)
dt
x3
1
0 1
x3
0
where x = [x1 , x2 , x3 ]T are the three state variables, u(t) is the input.
We perform a series of simulations of (4.1) by applying several random sinusoidal inputs.
(t) are plotted in Fig. 4.1(a) and Fig. 4.1(b), respectively. In
The trajectories of x(t) and dx
dt
this case, the trajectories lie exactly on a 2-dimensional linear subspace, and we expect that
the MOR algorithm is able to identify that linear subspace.
37
0
0.5
0.1
0.2
0.3
0.3
x3
0.4
0.2
0.4
0.1
5
4
3
0
0.5
0.2
0.4
0.3
0.4
0.6
0.1
0.2
0.1
0
x2
x1
0.2
x1
0.3
1
(a) x(t).
0.4
x
(b)
dx
dt (t).
x1
10 1
1
x1
0
1
d
x2 = 1 1 0 x2 0 + 0 u(t),
dt
x3
1
0 1
x3
x21
0
(4.2)
which adds a quadratic nonlinear term to the LTI system (4.1). Simple as it is, it exhibits
very different behaviors from the LTI system (4.1).
We apply the same random inputs to (4.2) (as those to (4.1)), and the trajectories of
(t) are plotted in Fig. 4.2(a) and Fig. 4.2(b), respectively. We observe that x(t)
x(t) and dx
dt
dx
and dt (t) no longer stay on a linear subspace. Rather, there seem to be a nonlinear manifold
attracting the trajectories.
This example suggests that nonlinear projections might be better for nonlinear model
order reduction, and motivates us to consider the nonlinear projection framework.
4.2
In the nonlinear projection framework, we try to identify two manifolds V and W that
are sub-manifolds of RN . To generalize the linear projection framework, we define x = v(xr )
38
0.05
0
0.2
x3
0.05
x3
0.15
0.1
0.15
0.1
0.2
0.05
0
0
0.5
0.2
0.4
0.3
0.4
0.2
0.6
0.1
0
x
x1
0.1
(a) x(t).
(b)
0.2
0.3
0.4
dx
dt (t).
39
Algorithm 4 MOR based on Nonlinear Projection with w(r) = W T (x)r
Inputs: dtd q(x) + f (x) + b(t) = 0
Outputs: W T (v(xr )) dtd q(v(xr )) + f (v(xr )) + b(t) = 0
1: Construct a mapping v : RNr RN ;
2: Construct a mapping w(r) = W T (x)r where W T : RN RNr N ;
3: Define w(r(v(xr ))) = 0.
how to identify manifolds V and W. Similar to LTI MOR, the manifold to be used in MOR
depends on the reduction criterion F (, ), G(, ) and H(, ) in (2.1). In this section, rather
than being very general to answer this question, we study some candidate manifolds, and
try to draw insights and conclusions.
The notations and simplifications that will be used in this section include:
1. We consider reducing ODE models, instead of DAE models, i.e.,
d
x = f (x) + Bu.
dt
2. Suppose x V, i.e., x = v(xr ). We denote V (x) for all x = v(xr ) by
dv
.
V (x) =
dxr x=v(xr )
(4.4)
(4.5)
(4.6)
(4.7)
where
fr (xr ) = W T (v(xr ))f (v(xr )),
(4.8)
Therefore, W is determined by V.
4.3
The first candidate manifold we study is the one that integrates local Krylov subspaces.
For nonlinear systems, properties are often described in terms of linearized systems. For
example, TPWL [2] tries to find a linear subspace which covers the local Krylov subspaces
40
around all expansion points. The linear subspace, however, is likely to grow very quickly as
the number of expansion points increases.
In order to avoid the growth of the order of reduced models due to the inclusion of many
linear subspaces, we expect to construct a manifold V so that at each point x V, the
tangent space around x, Tx V matches the local Krylov subspace of the linearized system at
x.
In another word, the mapping v : RNr RN should satisfy
!
v
= KNr ,x ,
(4.9)
span(V (x)) = span
xr
x=v(xr )
where KNr ,x is the Nr -dimensional Krylov subspace of the linearized system around x. Such
a manifold is called an integral manifold.
To see why the manifold satisfying (4.9) is useful, we show that responses of the full
linearized model and the reduced linearized model are closely related.
By linearizing (4.4) around xi RN , (i.e., let x = xi + x, with higher-order terms of
x truncated), we obtain
d
x = G(xi )x + f (xi ) + Bu,
dt
(4.10)
where
f
(x).
x
In contrast, the linearized model of (4.7) around xr,i
G(x) =
(4.11)
d
xr = Gr (xr,i )xr + fr (xr,i ) + Br (xr,i )u,
dt
(4.12)
where1
Gr (xr ) =
T
W (v(xr ))[f (v(xr )) + Bu]
xr
W T (v(xr ))
(f (v(
x)) + Bu),
xr
(4.13)
T
Br (xr ) =
W (v(xr ))[f (v(xr )) + Bu]
u
= W T (v(xr ))B.
To compare (4.10) with (4.12), let xi = v(xr,i ), and consider xi that is an equilibrium
T
W T (v(xr ))
xr
(f (v(
x)) + Bu) denotes
41
point, i.e.,
f (xi ) + BuDC = 0,
(4.14)
4.4
The manifold defined in the last section would be useful for model reduction. However,
the key question is: can we find an Nr -dimensional manifold satisfying conditions (4.9)?
Unfortunately, generally the answer is no. To show that, we resort to a fundamental result
in differential geometry [73, 112, 113], the Frobenius theorem.
To explain the result of Frobenius theorem, we introduce a few terminologies in differential
geometry2 :
1. Let X denote a manifold.
2. A distribution on X is a map which assigns to each point p in X a subspace of
the tangent space at this point, i.e.,
p X (p) Tp X.
3. A foliation {S }A of X of dimension k is a partition
[
X=
S
(4.15)
(4.16)
of X into disjoint connected sub-manifolds S , called leaves, which has the following
property: for any p X there exists a neighborhood U of p and a diffeomorphism
: U V RN onto an open subset V such that
\
N
((U
S )cc ) = {p = (x1 , , xN ) V |xk+1 = ck+1
(4.17)
, , xN = c },
For self-contained ness, we provide a short introduction to differential geometry that is relevant to our
discussion in Appendix Chapter A.
42
foliation {S }A on X such that for any p X,
Tp S = (p),
(4.18)
g
f
(x)f (x)
(x)g(x).
x
x
(4.20)
43
In Section 4.3, we have specified the distribution to be the one assigning Krylov subspaces to each point on the manifold, i.e.,
K (x) , span(KNr ) = span(G1 (x)B, , GNr (x)B),
x = v(xr ), xr RNr .
(4.21)
x RN .
(4.22)
dx
+ x + Bu = 0,
dt
(4.23)
(4.24)
x1 x3
x3 0 x1
g
g(x) = x2 x3 , A(x) =
(x) = 0 x3 x2 ,
(4.25)
x
x1
1 0 0
1
B = 0 .
0
(4.26)
Therefore, the two-dimensional Krylov subspace is the span of {v1 = A(x)B, v2 = A(x)2 B},
where
x3
x3 + x1
.
v1 (x) = 0 , v2 (x) = x2
(4.27)
1
x3
The Lie bracket [v1 , v2 ] is
2x3
v1
v2
v1
v2 = 0 .
[v1 , v2 ] =
x
x
1
(4.28)
44
When x1 = 1, x2 = 1, x3 = 1, we have
1
2
v1 (x) = 0 , v2 (x) = 1 ,
1
1
2
[v1 , v2 ] = 0
/ span({v1 , v2 }),
1
(4.29)
4.5
While the distribution defined by (4.21) is typically not involutive, we may induce an
involutive distribution from it. Given a set of analytical vector fields {fu }, = span({fu }),
we may define
L(p) = span{g(p)|g Lie{fu }}
(4.30)
where Lie{fu } is the smallest family of vector fields that contains {fu } and is closed under
taking linear combinations and Lie bracket (or the Lie algebra of {fu }).
According to corollary A.7.2, the orbit of (p, {fu }) defines a sub-manifold whose tangent
space
Tp X = L(p).
(4.31)
In another word, L(p) defines an involutive distribution that covers the subspace defined
by at each point p. Therefore, L(p) corresponds to a well-defined sub-manifold.
In order to compute L(p), we can start with L(p) = span({fu }) = , and then iteratively
add Lie brackets of generator functions in L(p) into L(p) until L(p) is involutive.
The possible caveat of this approach is that there does not seem to have any bound of
the dimension of L(p). Therefore, even in the case where we want to match all 2-dimensional
Krylov subspaces, L(p) may end up being of a large dimension, making the manifold not so
efficient in preserving accuracy of linearized models.
4.6
We have seen that K in (4.21) is generally not involutive, and K -induced L(p) is involutive but may be of large dimension. However, the discussion gives us a better understanding
on how we may define a manifold (more precisely, a sub-manifold in RN ):
1. We may define an involutive distribution, or equivalently, its generator functions, from
which a sub-manifold can be derived;
2. We may directly define the projection function x = v(xr ) that specifies the manifold
v
V, and derive its tangent spaces V (x) = x
(v(xr )).
r
45
In this section, we consider several ways to construct integral manifolds, and consider
special kinds of distributions/manifolds that may be analyzed or computed efficiently.
In the following discussion, we assume the projections x = v(xr ) and xr = w(x) such
that xr = w(v(xr )) for all xr . Note that this definition leads to the reduced model
dxr
dw dx
dw
=
=
(v(xr ))(f (v(xr )) + Bu),
dt
dx dt
dx
(4.32)
2. xr = x21 + x22 + x23 defines equivalent points in R3 being points on spheres of radius
(xr 0).
4.6.1
xr
Involutivity Constraints
i, j, i 6= j.
(4.33)
vj
vi
vi
vj span(V (x)),
x
x
i, j, i 6= j.
(4.34)
To simplify (4.34), we note that from the proof of the Frobineous theorem, if span(V (x))
is involutive, then there exists a commuting basis of generators for span(V (x)). Therefore,
without loss of generality, we assume that vi (x) and vj (x) are commuting, leading to algebraic
equations
vj
vi
[vi , vj ] =
vi
vj = 0, i, j, i 6= j.
(4.35)
x
x
3
46
For each pair (i, j), (4.35) correspond to N equations. Since there are 12 Nr (Nr 1) (i, j)
pairs, the involutivity constraints consist of 12 Nr (Nr 1)N equations.
4.6.2
Method 1
(4.36)
which states that the projection from xr to x and then x to xr leads to the same point on
the manifold.
With (4.36), we may immediately prove that for any x = v(xr ), the distribution defined
by span(V (x)) is involutive, and that the basis given by V (x) is commuting.
Theorem 4.6.2 If xr = w(v(xr )), then for all x = v(xr ), [vi (x), vj (x)] = 0 for all i, j, i 6= j.
Proof By definition,
Therefore,
v
,
vi (x) =
xr,i xr =w(x)
v
vj (x) =
.
xr,j xr =w(x)
vj
vi
vi
vj
x
x
2 v w v
2 v w v
=
.
xr xr,j x xr,i xr xr,i x xr,j
(4.37)
[vi , vj ] =
(4.38)
w v
,
x xr
x = v(xr ),
w v
= ei .
x xr,i
(4.39)
(4.40)
47
Therefore, (4.38) simplifies to
2v
2v
ei
ej
xr xr,j
xr xr,i
v
v
=
xr,j xr,i
xr,i xr,j
=0
[vi , vj ] =
(4.41)
While we havent shown that V (x) is involutive for all x RN , we see that by definition,
we have an integral manifold V that satisfies Txr V = span(V (x)). Therefore, the projections
v() and w() lead to a well-defined integral manifold.
However, in practice, it can still be hard or impossible to construct analytical expressions
of the nonlinear functions v() and w() that satisfy (4.36), except for a few special cases.
Example 4.6.3 Consider linear functions
x = v(xr ) = x0 + V xr ,
xr = w(x) = xr,0 + W T x,
(4.42)
(4.43)
(4.44)
W T V = I.
Therefore, given V, x0 , one of the many choices of W T and xr,0 is
W T = (V T V )1 V T ,
(4.45)
xr,0 = W T x0 .
Example 4.6.4 Consider a nonlinear function v() and a linear function w():
x = v(xr ) = [V, Vg ][xr ; g(xr )] = V xr + Vg g(xr ),
xr = w(x) = W T x,
(4.46)
48
In this case, (4.36) becomes
xr = W T (V xr + Vg g(xr )).
(4.47)
(4.48)
xr = w(x) = W T x = [x1 , x2 ].
(4.49)
In (4.49), we define the manifold by v() which is the surface defined by x1 = xr,1 , x2 =
xr,2 , x3 = xr,1 xr,2 . On the other hand, we define the projection of each point x to the manifold
by w(): for any point (x1 , x2 , x3 ), the projection leads to the point v(w(x)) = (x1 , x2 , x1 x2 ).
(4.48) corresponds to V = W = [I2 ; 0], Vg = e3 . As a result, W T Vg = 0. Therefore, V (x)
defines an integral manifold.
Note also that in this example, the manifold (defined by v()) is nonlinear, but the projection (defined by w()) is linear.
Example 4.6.5 Consider a linear function v() and a nonlinear function w():
x = v(xr ) = V xr ,
(4.50)
(4.51)
(4.52)
To simplify the problem, we further assume g(x) are quadratic functions, i.e., g(x) =
x x. Then,
WgT (V xr ) (V xr ) = WgT (V V )(xr xr ) = 0,
(4.53)
for which a sufficient condition is
WgT (V V ) = 0.
4
5
(4.54)
49
For example, consider
V = [1; 0; 0],
W = [1; 0; 0].
(4.55)
We then have
and a valid choice of WgT is
V V = [1; 0; 0; 0; 0; 0; 0; 0; 0],
(4.56)
(4.57)
x = v(xr ) = [xr ; 0; 0]
xr = w(x) = x1 + x23
(4.58)
In this case, the manifold is the x1 axis (a linear subspace), but the projection of any
point to the x1 axis is the parabola x1 + x23 = 0 (nonlinear).
4.6.3
Method 2
In method 2, we start with an involutive and commuting basis V (x), from which we
induce the mapping v() by solving
dx
= V (x).
dxr
(4.59)
dx
= vi (x).
dxr,i
(4.60)
(4.61)
Since columns of V (x) commute, we can change the order of the integration without
changing the mapping. With V (x), we can also induce a W (x) by ensuring W T (x)V (x) = I.
In this method, we have not explicitly defined the projection xr = w(x) although it may
= W T (x). However, we do not have to compute w(x) for the
be induced by solving dw(x)
dx
purpose of reduced order modeling W T (x) alone is enough to construct a reduced order
model.
One way to find such a V (x) that also preserves system responses (such as moments
of linearized models) is to start with a non-involutive and non-commuting accurate basis
50
(e.g., the Krylov basis), and then try to find a closest distribution to it.
Another way is to consider specific basis functions, and derive simpler conditions for
involutivity and commutativity.
Example 4.6.6 Consider V (x) = [v1 , , vNr ] where vi (x) is in the form of
vi (x) = Ai x + bi ,
i = 1, , Nr ,
(4.62)
where Ai RN N and bi RN 1 .
Applying the involutivity constraint in (4.35) on the pair (vi , vj ), we obtain
[vi , vj ] = Aj (Ai x + bi ) Ai (Aj x + bj )
= (Aj Ai Ai Aj )x + (Aj bi Ai bj ) = 0.
(4.63)
(4.64)
Ai = Ui U 1 ,
(4.65)
(4.66)
4.6.4
(4.67)
Method 3
In method 3, we start with an arbitrary distribution defined by U(x) that may be noninvolutive or non-commuting. However, we can still define a mapping v : xr 7 x by
v : (xr,1 , , xr,Nr ) 7 exp(xr,Nr uNr ) exp(xr,1 u1 ).
(4.68)
(4.69)
51
in a specified order (i = 1, , Nr ).
From this mapping, we may then derive the distribution
dv(xr )
V (x) =
,
dxr xr =w(x)
(4.70)
4.7
(4.71)
4.7.1
Symbolic/Analytical Representation
If v(xr ), w(x) have analytical expressions, we may simply store their symbolic expression.
In particular, some special functions even have a matrix representation, such as quadratic
functions
v(xr ) = G1 xr + G2 xr xr ,
(4.72)
for which we may just store the matrices G1 and G2 .
Moreover, analytical expressions of v() and w() may help simplify (4.71). In the linear
projection framework, we have
fr (xr ) = W T GV xr ,
(4.73)
where W T GV RNr Nr can be precomputed and stored, and is N-independent. We may
do the similar trick in the nonlinear case.
Example 4.7.1 Consider the case where
v(xr ) = V1 xr + V2 xr xr ,
f (x) = G1 x + G2 x x,
W T (x) = W T .
(4.74)
52
We obtain
fr (xr ) =W T (G1 (V1 xr + V2 xr xr ) + G2 (V1 xr + V2 xr xr ) (V1 xr + V2 xr xr ))
=W T G1 V1 xr + W T G1 V2 xr xr
(4.75)
=W T G1 V1 xr + (W T G1 V2 + W T G2 (V1 V1 ))(xr xr )
+ W T G2 (V1 V2 + V2 V1 )(xr xr xr )
+ W T G2 (V2 V2 )(xr xr xr xr ),
where we only have N-independent matrices to pre-compute and store in the reduced model.
4.7.2
When no analytical expression is available, we have to use certain approximation for the
nonlinear mappings and nonlinear functions. One such approximation is PWL approximation6 , i.e., concatenating (interpolating) point-wise linearized functions to obtain a global
function.
Approximating mapping functions
In PWL approximation, we assume that, instead of the actual nonlinear function x =
v(xr ), we are given a few pairs of (xr,i , xi ), i = 1, , Ns where Ns is the number of samples
on the manifold. Therefore, we can approximate v() by
x=
Ns
X
Ki (xr )xr,i ,
(4.76)
i=1
Ns
X
Ki (xr ) 0,
Ki (xr ) = 1,
(4.77)
i=1
Ki (xr,i ) = 1.
When we also have information about the tangent space (e.g., Krylov subspaces) at each
6
While the method is traditionally called piecewise-linear in MOR community, it is also known as a
special case of kernel methods.
53
point, i.e.,
v
Vi =
,
xr x=v(xr,i )
(4.78)
Ns
X
i=1
(4.79)
Ns
X
i=1
(4.80)
fr (xr ) W (xr )
Ns
X
i=1
(4.81)
Ns
X
j=1
i=1
Ns
X
i=1
Ns
X
Kif (xr ) fi + Gi (
j=1,j6=i
Ns
X
Kif (xr ) fi + Kiv (xr )Gi Vi (xr xr,j ) + Gi (Kiv (xr )xi xi )
(4.82)
xr,j ) + xj ) xi ) .
where Kif () and Kiv () are the weighting functions for f () and v(), respectively.
Note that when xr is close to xr,i , Ki (xr ) 1 and Kj (xr ) 0, j 6= i. Therefore, we have
fr (xr ) W T (xr )(fi + Gi Vi (xr xr,i )).
(4.83)
54
Therefore, the PWL approximation of f () is often simplified to
T
fr (xr ) W (xr )
Ns
X
i=1
(4.84)
which can be viewed as approximately truncating second-order terms of Kif (xr )Kjv (xr ).
To approximate W T (xr ), we differentiate (4.80) and obtain
w
W T (xr ) =
.
(4.85)
x x=v(xr )
Similarly, a PWL approximation of W T (xr ) is
W T (xr )
Ns
X
wi (xr )WiT .
(4.86)
i=1
Combine (4.84) with (4.86), and use a simplified PWL approximation, we obtain
fr (xr )
Ns
X
i=1
(4.87)
Therefore, we can pre-compute and store WiT fi RNr and WiT Gi Vi RNr Nr in the
reduced model. As a result, the computation cost of the reduced model (i.e., (4.87)) is
depends only on Nr and Ns , and is N-independent.
Approximating Jacobians
Taking the derivative of (4.87) with respect to xr , we obtain the Jacobian
N
s
s
X
fr X
Ki (xr ) T
Wi fi + WiT Gi Vi (xr xr,i ) .
Ki (xr )WiT Gi Vi +
xr
xr
i=1
i=1
(4.88)
Note that when Kxi(xr r ) is small (which is true for xr far away from xi ), we can further
ignore the second part in (4.88), and obtain
N
s
X
fr
Ki (xr )WiT Gi Vi .
xr
i=1
(4.89)
55
4.7.3
To generate the reduced order model, we pre-compute and store the reduced matrices/vectors (e.g., those used in (4.87)) in a look-up table. They are summarized in Table
4.1.
Table 4.1: Look-up table stored in the PWL-based reduced model
Variable Expression Size
fri
WiT f (xi )
RNr 1
T
qri
Wi q(xi )
RNr 1
Bri
WiT B
RNr 1
f
Gr i
WiT x (xi )Vi RNr Nr
q
(xi )Vi RNr Nr
C ri
WiT x
Therefore, if these is Ns samples, the time and space complexity of the reduced model
(i.e., computation and storage of (4.87)) is O(Ns Nr2 ) it is N-independent and therefore is
more efficient than computing the full model given Nr N and Ns N.
Furthermore, instead of using all Ns samples for the computation at a specific point
xr , we may use k-nearest neighbors of xr in the computation of (4.87) because the weight
function wi (xr ) evaluates approximately to 0 when xr is far away from xi . Therefore, the
computational complexity is further reduced to O(kNr2 ) plus the complexity of computing
the k-nearest neighbors.
56
Chapter 5
ManiMOR: General-Purpose
Nonlinear MOR by Projection onto
Manifolds
In this chapter, we present a general-purpose nonlinear MOR technique based on the
nonlinear projection framework presented in Chapter 4. The algorithm explicitly constructs
and parameterizes a manifold in the state space, and projects the original nonlinear model
to the manifold to obtain a nonlinear reduced order model.
We describe in detail how we construct and parameterize a general-purpose manifold
that captures typical system responses (such as DC and AC responses). We validate and
illustrate the MOR technique using several examples.
5.1
Many dynamical systems derived from engineering problems work in certain highly constrained region in the state space, either because the fast dynamics die out very quickly and
are not easily observable, or because the systems are designed/engineered to work under
certain mode (e.g., amplifiers are designed to work around a DC steady state / equilibrium,
and oscillators are designed to work around a limit cycle). In another word, state variables
typically stay close to the attractors, and that part of the state space covers important dynamics. In order for the reduced model to reproduce the attractors, a necessary condition is
that the manifold must cover the attractors.
5.1.1
DC Manifold
We first consider the manifold covering equilibria (also known as DC operating points in
circuit community) of a system.
57
For the nonlinear system (1.2), the DC equilibria are the solution to
f (x) + Bu = 0,
u [umin , umax ],
(5.1)
where umin and umax are bounds of the input u according to physical constraints.
For simplicity, we consider the case where f () is invertible. We consider this assumption
because it is often satisfied in a large number of circuit examples.
As a result of this assumption, we have
x = f 1 (Bu),
(5.2)
x(u = 0) = f 1 (0),
(5.3)
where fb () is the Jacobian of f 1 (). (5.3) is obtained by taking the partial derivative of
(5.2) with respect to u.
Similarly, we may directly take the partial derivative of (5.1) with respect to u, and
obtain
x
G(x)
+ B = 0,
(5.4)
u
and therefore we can define
x
= G1 (x)B,
(5.5)
v1 (x) =
u
which associate each point on the DC manifold with its tangent space.
The RHS of (5.5) is exactly the direction of the first Krylov vector of the linearized model
around x. Therefore, as proved in [89], at each DC equilibrium, the linearized reduced model
matches the first moment of the linearized full model.
1
5.1.2
AC manifold
58
projection matrix VNr (x) that satisfies
span(VNr (x)) = span(G1 (x)B, G2 (x)B, , GNr (x)B).
(5.6)
Because the tangent space around a point x gives a good local approximation of the
manifold, the points close to x is, loosely speaking, also on the manifold. Therefore, one
may grow the manifold by integrating along the Krylov vectors to obtain the manifold.
That is, solving
x
(5.7)
= VNr (x)
xr
to obtain the projection x = v(xr ).
However, according to the discussion in Section 4.4, in general there may not exist a
solution to (5.7). Nevertheless, we can always ensure that at each point on the DC manifold,
(5.7) holds (one example is given in Section 5.1.3. We will call any manifold that satisfies
this condition an AC manifold, because such manifolds lead to reduced order models that
match up to the Nr -th moment of the linearized full models, and therefore the AC responses
are well approximated.
5.1.3
While we do not have a manifold whose tangent spaces match local Krylov subspaces, we
can construct manifolds that cover DC and AC manifolds. In the following, we provide one
such solution by constructing parametric equations of the manifold in terms of Nr variables
xr = (xr,1 , , xr,Nr ).
Let VNr (x) = [v1 (x), , vNr (x)] be the Krylov basis of the linearized system around x,
and let x0 = f 1 (0) be the DC solution for u = 0. We grow the manifold V starting from
x0 :
1. Solve for x = v(xr,1 , 0, , 0) by integrating
x
= v1 (x),
xr,1
x(0) = x0 .
(5.8)
(5.9)
59
4. Solve for x = v(xr,1 , xr,2 , , xr,Nr ) by integrating
x
= vNr (x),
xr,Nr
(5.10)
(5.11)
In another word, the manifold is composed of points that are solutions of (5.11).
The rationale behind the construction of such a manifold is that at a point x in the
state space, the state variable can move in N independent directions, corresponding to N
modes which are approximated by the N Krylov basis vectors. Ignoring N Nr fast
modes, the state variable will first move along the Nr -th direction (the fastest remaining
mode), and then the (Nr 1)-th direction (the second fastest remaining mode), and so on.
In the end, it moves along the first direction (the slowest mode). Therefore, at each point on
the Nr -dimensional manifold, we only guarantee that vNr (x) is covered by the local tangent
space.
This leads to the conclusion that at each point on the k-dimensional manifold, the local
tangent space covers span(vk (x), , vNr (x)). In particular, at each DC equilibrium, the
local tangent space covers the Krylov subspace of the linearized model. Therefore, the
resulting nonlinear model covers the DC and AC manifolds.
5.2
Parameterization
(5.12)
d
x( ) = ( )g(
x( )),
(5.13)
d
respectively, where t = ( ) is a function of , and g() is Lipschitz. Then x(t) and x( )
define the same points in the state space, i.e.t, t, such that x(t) = x(t).
60
Proof Since t = ( ), we have
dt = ( )d
(5.14)
(5.15)
d
x( ) dt
d
x( ) =
= g(x(t)) ( ) = g(
x( )) ( ),
d
dt d
(5.16)
Define
Therefore,
i.e.,
d
x( ) = ( )g(
x( )).
(5.17)
d
According to the existence and uniqueness theorem for ordinary differential equations,
x(t) and x( ) are the unique solutions to (5.12) and (5.17), respectively.
Therefore, x(t) and x(t) define the same points in the state space.
Based on theorem 5.2.1, the following corollary follows.
Corollary 5.2.2 The solutions to the normalized integral curve equation
x
vi (x)
=
xr,i
||vi (x)||
(5.18)
define the same trajectory as the solutions to the regular integral curve equation
x
= vi (x).
xr,i
(5.19)
Proof Suppose the solutions to (5.19) and (5.18) are x(t) and x(t), respectively. Let
t = (t) =
1
d.
||vi (
x())||
(5.20)
By theorem 5.2.1, x(t) and x(t) span the same state space.
From the above discussion, we conclude that different choices of (t) leads to different
parameterizations of the manifold. Nevertheless, they all define the same manifold. Therefore, we may define (t) to obtain better numerical accuracy in computing the manifold and
the projection function.
61
5.3
5.3.1
Numerical Computation
It can be very expensive to use (4.68), or (5.11), to compute the projection functions and
the reduced nonlinear models, because they themselves involve solving a series of differential
equations.
x
For example, to compute V (x) = x
, we have to solve
r
x
= vNr (x),
xr,Nr
dx
dxr,Nr 1
x
x0,xr,Nr 1
vNr 1 (x),
dx
dxr,Nr 2
x
x0,xr,Nr 2
vNr 2 (x), ,
(5.21)
where x0,xr,i = x(xr,1 , , xr,i , 0, , 0), and the computation of RHS in (5.21) involves
expensive transient sensitivity analysis.
To alleviate the computational cost of the reduced model, we employ the piece-wise linear
approximation method described in Section 4.7.2 that makes the time and space complexity
O(Ns Nr2 ).
5.4
62
5.4.1
5.5
In this section, we apply ManiMOR to several examples, including artificial simple nonlinear systems, practical circuits such as ring mixer, nonlinear transmission line and I/O
buffer circuit, as well as a bio-chemical signaling pathway example. We validate the ManiMOR model by comparing the simulation results using ManiMOR model and the full model.
63
We show comparisons of ManiMOR to existing methods, mainly against the TPWL method,
to examine the accuracy and efficiency of ManiMOR.
5.5.1
In the first example, we apply ManiMOR to the simple nonlinear system (4.2) mentioned
in Section 4.1, in order to illustrate how ManiMOR works.
We apply the ManiMOR algorithm to generate a 2-dimensional reduced order model for
this 3-dimensional system. We then perform transient analysis of both models and compare
the results.
The simulation results obtained by applying a slowly changing input u(t) = t are shown in
Fig. 5.1. We observe that the trajectories of the reduced order model follow the DC manifold
in the state space. The reduced model covers the DC manifold and therefore reproduces the
slowest dynamics.
0
reducedorder model
full model
projection on manifold (full model)
0.05
x3
0.1
0.15
0.2
0.25
0
10
20
30
40
50
time
(b) x3 (t)
Figure 5.1: Simulation of the ManiMOR model. (u(t) = t) Fig. 5.1(a) shows the manifold
ManiMOR identifies (yellow), and the trajectories of the full model (blue) and ManiMOR
model (green). Fig. 5.1(b) shows the waveform of x3 (t) of the ManiMOR model (green), the
full model (blue), and the projection of the full model solution onto the manifold (black).
The simulation results obtained by applying a sinusoidal signal u(t) = 2.5 + 2sin(0.2t)
are shown in Fig. 5.2. Because the frequency of the input signal is low, the trajectory of
the full model is nearly on the manifold. In this case, the trajectory stays close to the DC
solution, but nonlinear distortion is exhibited in steady state solution. Since the manifold
in ManiMOR is good in capturing behaviors around DC solutions, this type of distortion is
well approximated by the ManiMOR model, as shown in Fig. 5.2. The maximum absolute
and relative mean squared error of ManiMOR, compared to the full model, are 0.0136 and
0.0473, respectively.
64
0.05
reducedorder model
full model
projection on manifold (full model)
x3
0.1
0.15
0.2
0.25
0
10
time
15
20
(b) x3 (t)
Figure 5.2: Simulation of the ManiMOR model. (sinusoidal input) Fig. 5.2(a) shows the
DC curve (red) and the manifold ManiMOR identifies (yellow), the trajectories of the full
model (blue) and ManiMOR model (green). Fig. 5.2(b) shows the waveform of x3 (t) of the
ManiMOR model (green), the full model (blue), and the projection of the full model solution
onto the manifold (black).
Fig. 5.3 and Fig. 5.4 show simulation results to compare the ManiMOR model and the
TPWL model. Fig. 5.3 shows the simulation results obtained by applying a multiple-step
function
1, t < 5,
5, 5 t < 10,
u(t) =
(5.22)
2, 10 t < 15,
4, 15 t 20.
We observe that the ManiMOR model tracks the trajectory of the full model better than
TPWL. Better means that 1) the distance (error) between the trajectories of ManiMOR
model and the full model is smaller than that of the TPWL model; 2) the ManiMOR model
converges to the right DC solution, while the TPWL model fails to get back to the right
operating point after transient response dies out, which is better observed in Fig. 5.3(b).
In this example, the maximum absolute mean square error of the ManiMOR model and
TPWL model, compared to full model, are 0.0534 and 0.0627, respectively. Although the
TPWL model only leads to an error that is slightly larger than the ManiMOR model, the
trajectory of TPWL model fails to converge to the right DC operating point this can cause
serious problems when using reduced models because the attractor is not preserved. Such
models may not be acceptable in practice.
Similar results are observed when a two-tone sinusoidal input u(t) = 2.5 + sin(0.1t) +
cos(0.4t) is applied, as shown in Fig. 5.4. Although the waveforms obtained using the
TPWL model looks similar to those of the full model, there is an observable offset. In
65
0
maniMOR model
TPWL model
full model
0.05
0.1
x3
0.15
0.2
0.25
0.3
0.35
0
10
time
15
20
(b) x3 (t)
Figure 5.3: Comparison of ManiMOR and TPWL model. (multiple-step input) Red, green
and blue trajectories represent simulation results of ManiMOR, TPWL and full model,
respectively.
contrast, the ManiMOR model tracks the nonlinear dynamics correctly.
5.5.2
The second example, a CMOS ring mixer [1] (shown in Fig. 5.5), is highly nonlinear.
When the input changes, the DC operating point moves from one state to another along a
highly nonlinear curve, as shown by the red curve in Fig. 5.6(a). The 2-D manifold identified
by ManiMOR is also depicted in Fig. 5.6(a). Clearly the highly nonlinear DC curve is not on
a linear subspace, and the 2-D manifold is extremely twisted. In fact, regions that are near
two ends of the DC curve are relatively linear (the two linear spaces are almost orthogonal to
each other), while the region in the middle shows strong nonlinearity, and smoothly connects
two nearly-orthogonal linear subspaces. This again suggests that TPWL, which asserts that
the solutions will be close to a linear subspace, will be inefficient to reduce the full model.
To examine how well ManiMOR captures dynamics and nonlinearities in this example,
we apply a step input, and simulate the circuit using the full model and the ManiMOR
model. The simulation results are shown in Fig. 5.6, where the trajectory starts from one
DC operating point, goes out of the manifold, and finally converges to another DC operating
point. Again, the ManiMOR model shows close match to the full model and qualitatively
correct behaviors. The maximum absolute and relative mean square error in this case is
0.1293 and 0.0492, respectively.
66
0
maniMOR model
TPWL model
full model
0.05
x3
0.1
0.15
0.2
0.25
0
10
time
15
20
(b) x3 (t)
Figure 5.4: Comparison of ManiMOR and TPWL model. (two-tone sinusoidal input) Red,
green and blue trajectories represent simulation results of ManiMOR, TPWL and full model,
respectively.
67
4.5
reducedorder model
full model
projection on manifold (full model)
4
3.5
x3
3
2.5
2
1.5
1
0
0.5
1
time
1.5
2
3
x 10
(b) x3 (t)
Figure 5.6: Simulation of ManiMOR model and full model for the CMOS ring mixer. (step
input) Fig. 5.6(a) shows the manifold ManiMOR identifies (yellow), and the trajectories of
the full model (blue) and ManiMOR model (green). Fig. 5.6(b) shows the waveform of x3 (t)
of the ManiMOR model (green), the full model (blue), and the projection of the full model
result onto the constructed manifold (black).
5.5.3
A test example that is examined in most papers about trajectory-based methods [2, 114]
is a nonlinear transmission line circuit. To verify ManiMOR, we apply both ManiMOR and
TPWL on the nonlinear transmission line circuit with 100 nodes.
Using the same training set and model size in [2] (i.e., TPWL model is of order 10;
training input is the step input i(t) = u(t 3)), we apply a test input of i(t) = 0.8u(t 3)
to the reduced models and the full model. As shown in Fig. 5.7, since the DC solution
corresponding to i(t) = 0.8 is not trained, an observable error, compared to the full model,
presents in the result of the TPWL model.
Accordingly, a ManiMOR model of order 4 is generated, and the results are shown in
Fig. 5.7. We observe that (1) the error is less than that of TPWL model, even if the model
size is almost halved; (2) despite the initial transient error, the trajectory converges to the
exact steady state solution.
5.5.4
Another test circuit is a current-mode logic (CML) buffer chain [115], which has also
been studied in previous MOR literatures [116, 117]. The circuit diagram of the CML buffer
is shown in Fig. 5.8. We use the BSIM3 [118] model for MOSFETs in our simulations, and
the size of differential algebraic equations for the full circuit is 52.
Applying ManiMOR, a size 5 model is generated for this circuit. The nonlinear manifold
68
0.012
0.01
v1
0.008
0.006
0.004
Full model
TPWL q=10
maniMOR q=4
0.002
0
0
10
time
Figure 5.7: Comparison of of the circuit response of ManiMOR model and TPWL model.
(Nonlinear transmission line circuit [2])
Vdd
Vout
Rload
Vin
Vbias
69
constructed by ManiMOR is visualized in Fig. 5.9 (only first two dimensions are plotted)
by picking three node voltages as axes. Clearly, the manifold is nonlinear, and blue circled
points, which are sampled on a transient trajectory, do stay close to this nonlinear manifold.
This confirms that the manifold attracts system responses.
Figure 5.9: Visualization of the low-order manifold subspace generated by ManiMOR. The
axes represent three node voltages in the circuit. The red curve consists of the DC operating
points; the black dashed lines are the second local subspace bases; the yellow manifold is
generated by stitching together all the local linear subspaces. The blue circled points are on
a transient trajectory simulated using the full model.
In contrast, the TPWL model aggregates all the samples and all projection matrices into
the final projection matrix. We used 20 DC operating points for training by aggregating
only the state variables, we obtain a matrix of 20 columns. The singular values of this
matrix are shown in Fig. 5.10. Therefore, to get a set of bases which approximate all the DC
operating points with an overall accuracy of 0.01, first 8 bases must be used. This implies
that the size of the TPWL model must be at least 8, only to capture the DC response.
When generating the TPWL model, we use the same method (i.e., Arnoldi algorithm),
70
10
log()
10
20
30
40
0
10
ith singular value
15
20
Figure 5.10: Singular values of the matrix containing 20 equilibrium points, sorted.
the same linear subspace size (i.e., 5), to calculate the local linear projection matrices as
in ManiMOR. As a result, a size 12 TPWL model is generated in order to match all the
DC operating points, and the first 5 moments of the transfer functions for all the linearized
systems.
In this example, we also generate a small-signal linearized model of size 5 we linearize
the nonlinear model around its DC solution, and generate a linear reduced model from that
linear system.
We first compare the ManiMOR model to the small-signal linearized model, by applying
a small signal input superimposed on a DC bias, u(t) = 0.005 sin(109 t)+0.0025 sin(1/3
109 t). Not surprisingly, since the amplitude of the input is small, the transient trajectory
is close to the DC operating point, and the linear small-signal reduced order model almost
tracks the full model, as shown in Fig. 5.11(a). For the same reason, the result of the
ManiMOR model also matches that of the full model.
However, when the input amplitude is large, the trajectory will not be confined close to a
DC operating point, but travel among different regions. In that case, the small signal model
fails to match the waveforms of the full model. To see that, we apply a sinusoidal input with
larger amplitude u(t) = 0.25 sin( 109 t). As shown in Fig. 5.11(b), the output waveform
is highly distorted. However, the waveform of the ManiMOR model is almost overlapped
with that of the full model this is because the ManiMOR model uses different local models
when the state variable travels to different regions in the state space. On the contrary, the
small signal model, which makes the linear and small-signal assumption, still generates
a sinusoidal output waveform, and therefore fails to capture the distortion effect.
We then compare the ManiMOR model to the TPWL model by applying a step input
71
2.5
1.6
full model
linearized model
maniMOR model
1.59
1.58
2
Vout1
Vout1
1.57
1.56
full model
linearized model
maniMOR model
1.5
1.55
1.54
1.53
1.52
0
0.5
1
time
1.5
0.5
0
0.2
0.4
0.6
time
x 10
0.8
1
8
x 10
Figure 5.11: Comparison of ManiMOR model with linearized reduced-order model. The
waveforms of one output using different models are plotted.
as shown in Fig. 5.12(a), which is chosen to traverse four DC operating points. As shown in
Fig. 5.12(b), without the heuristic to handle the DC solutions, the TPWL model with model
size 10, fails to converge back to the correct DC operating points. After including the DC
operating points into the projection matrix, the TPWL model of size 10 still leads to some
error. During experiments, we find that the model size must be at least 12 in order to match
the original trajectory. This result is obtained by using a training input to be the same as
the test input potentially an even larger model will be needed.
In comparison, ManiMOR model with size 5 renders an output waveform almost indistinguishable from that of the full model it is less than half of the size of the TPWL model
with the same accuracy. We conclude that the redundancy in TPWL is because of the aggregation and the SVD step which embed all the DC solutions and all the local Krylov
subspaces into a linear subspace. Under the linear subspace constraint, it is highly possible
for TPWL to obtain less reduction than ManiMOR.
5.5.5
In this example, we show the application of ManiMOR macromodeling method on a Biochemical system. We have found that ManiMOR model reduces the system size by a large
amount for this example.
In a chemical reaction, the number of reactants is usually large. For example, yeast
pheromone signaling pathway [3] contains a chain of chemical reactions, as shown in Fig. 5.13.
For each reaction, several reactants participate. Shown in [119, 120], 304 reactants in this
pathway are modeled, rendering a dynamical system of size 304. However, questions, such
72
0.3
0.2
Vin
0.1
0
0.1
0.2
0.3
0
0.5
1
time
1.5
2
8
x 10
1.8
full model
TPWL without heuristic
TPWL q=10
TPWL q=12
maniMOR q=5
1.7
1.6
Vout1
1.5
1.4
1.3
1.2
1.1
1
0
0.5
1
time
1.5
2
8
x 10
Figure 5.12: Comparison of ManiMOR model with TPWL model. The waveforms of one
output voltage using different models are plotted.
73
as how inputs affect specified outputs of the system, are of more interest to bio-chemical
scientists, and this is where MOR can play an important role. Reduced-order models also
makes it possible to simulate large scale bio-chemical systems.
Figure 5.13: Components of the pheromone response pathway in yeast. [3] Each arrow represents a type of reaction.
In the procedure of manifold construction in ManiMOR, we carefully simulate the DC
operating points due to a wide range of input values, since it is found that DC analysis
for this specific system is hard to converge unless continuation/limiting methods are used.
Moreover, it is observed that several equilibriums exist for some inputs, although some of
them are unstable and cannot occur in reality. (The equilibrium solutions of one output are
plotted in Fig. 5.14(a). In the process of building the ManiMOR model, we also implicitly
eliminates all the useless equilibriums, which makes the DC convergence of the reduced
model much better than the full model.
In contrast, TPWL model has to generate a matrix containing all the equilibrium points.
We picked up 20 equilibrium points into a matrix. The singular values of this matrix is again
plotted in Fig. 5.14(b). Setting the threshold of the singular value to be 1010 , we have to
include 8 basis vectors into the projection matrix, which renders at least a size 8 model.
74
1
0.9
0
5
log()
0.8
0.7
10
15
0.6
20
0.5
0.4
0
25
3
input
(a) Concentrations of x1 .
5
3
x 10
30
0
10
ith singular value
15
20
Since the first few basis vectors in the local projection matrix Vi are more important than others, we
construct Vk , whose columns are the aggregation of the k-th basis vector from all Vi at different samples.
Then SVD is performed to each Vk , and the basis vectors corresponding to singular values larger than a
threshold are chosen. For different k, the threshold changes due to their importance. Then, we simply
aggregate all Vk into Vagg , and do an SVD on Vagg = U V T just for ortho-normalization. The projection
matrix is then set to U . Using this heuristic, the important basis vectors are typically not filtered in the
aggregation and SVD step.
75
0.95
0.9
0.85
0.8
0.75
0
full model
TPWL model q=144
maniMOR model q=3
maniMOR model q=30
5
10
time (hour)
15
20
Figure 5.15: Transient simulation of the yeast pheromone pathway. The concentration of
one reactant is plotted.
76
Chapter 6
QLMOR: MOR using Equivalent
Quadratic-Linear Systems
In this chapter, we present another nonlinear projection-based MOR technique, based on
an equivalent quadratic-linear system representation of the original nonlinear system. This
representation is exact in the sense that the solution of the equivalent system is exactly the
same as the original nonlinear system. Therefore, the quadratic-linear model differentiates
itself from Taylor-expansion-based quadratic/polynomial models. The reduction is then
performed on the quadratic-linear representation to obtain a quadratic-linear reduced model.
We first summarize previous work on polynomial reduced order modeling techniques, and
the major problems therein. We then present algorithms that transform almost all nonlinear
differential equations to their quadratic-linear representation. With the resulting quadraticlinear representation, we perform Volterra analysis, and derive a reduced-order modeling
technique. We show experimental results and comparisons to other techniques.
6.1
Previous Volterra-series based MOR methods [106, 107, 121] are derived from the result
that the response of a nonlinear system can be decomposed into responses of a series of homogeneous nonlinear systems, i.e., the system response x(t) can be written as the summation
of responses of all n-th order homogeneous nonlinear systems xn (t):
x(t) =
X
n=1
xn (t),
(6.1)
77
where
xn (t) =
hn (1 , , n )
(6.2)
u(t 1 ) u(t n ) d1 dn .
hn (1 , , n )
(6.3)
(s1 1 ++sn n )
e
d1 dn ,
which is the multi-variable Laplace transform of hn (1 , , n ) [105].
Previous Volterra-series based methods can be summarized in four steps:
1. Perform a Taylor expansion of the nonlinear function f (x), i.e., expand f (x) in a
series of multidimensional polynomials, and truncate the expansion to a predefined
order (e.g., 3, as shown in (6.4)1 ).
C
d
x = G1 x + G2 x x + G3 x x x + Bu,
dt
(6.4)
d
x = Gx + Dxu + Bu,
dt
(6.5)
d
xr = Gr,1 xr + Gr,2 xr xr + Gr,3 xr xr xr + Br u
dt
(6.6)
78
by defining
x = V xr ,
Cr = V T CV,
Br = V T B,
Gr,1 = V T G1 V,
(6.7)
Gr,2 = V T G2 (V V ),
Gr,3 = V T G3 (V V V ).
The common step in [106, 107, 121] is Taylor approximation. This results in three major
drawbacks that limit the wide applicability of these methods. Firstly, the convergence of
the Volterra system representation is only valid when |u| , for some small > 0, i.e., for
small inputs, or the weakly nonlinear case. This constraint stems partly from the fact
that the Taylor series only converges in a small region around the expansion point. Incorrect
behavior of the reduced model due to large inputs was reported in [107], and similar results
are reproduced in Section 6.5.2.
Secondly, Taylor expansion typically destroys important global properties of the system,
such as passivity and stability. For example [107], the quadratic approximation for ex , which
is part of the I-V curve of a diode, around x = 0 leads to an artificial diode that can
generate energy when the voltage applied is negative. As a result, the reduced model, which
is created using the Taylor approximated system as a starting point, will typically not be
passive/stable.
Thirdly, the computational cost for higher-order terms such as Gr,3 xr xr xr dominates
3
the computational cost of the reduced model. Since xr RNr , Gr,3 RNr Nr , and with
Gr,3 usually dense, the computational cost is O(Nr4 ). (It is possible to make Gr,3 of a
smaller size by combining terms in xr xr xr , but the computational complexity does not
change.) Although it is N-independent, it can potentially be too expensive to use in practice.
Moreover, it is generally difficult to accurately extract the high-order matrices G2 , G3 , ,
which also constrains applications of Taylor-expansion based methods.
Remark MOR methods based on bilinear forms and Volterra series [107] actually define
a nonlinear projection e.g., by explicitly writing out the differential equations for ~x
~x, ~x ~x ~x, the final projection is defined by xr = V T [xT , xT xT , xT xT xT ]T , which
is a polynomial projection function. In this sense, QLMOR indeed allows a richer set of
projections than those in previous methods [107].
Briefly summarizing this section, we list in Table 6.1 the main drawbacks and difficulties
that have not been solved in existing MOR methods, and how well QLMOR addresses these
drawbacks.
79
Existing methods
Volterra series
based
Trajectory based
6.2
QLMOR
Solved
Partially solved
Improved
Improved
QLMOR Overview
In this section, we outline the key ideas of QLMOR. Details are provided in Section 6.3
and Section 6.4.
From the discussion in Section 6.1, the drawbacks of previous Volterra-series based methods are mainly caused by Taylor expansion step, in which information about the system
is lost.
In contrast, QLMOR derives a polynomial system that is equivalent to the orignal system
without losing any information. By doing so, QLMOR avoids the small input assumption,
and preserves properties of the original system in the first step.
We show that any polynomial ODEs can be converted into quadratic-linear differential
algebraic equations (QLDAEs)
C
d
x = G1 x + G2 x x + D1 xu + D2 (x x)u + Bu,
dt
2
(6.8)
where C RN N , G1 RN N , G2 RN N , D1 RN N , D2 RN N , B RN 1 . That
is, the differential equations are quadratic in state variables x and linear in the input u.
Similarly, any generalized polynomial system (to be defined in Section 6.3.8) can be
converted into generalized QLDAEs
d
[C1 x + C2 x x] = G1 x + G2 x x + (D1 x + D2 x x)u + Bu,
dt
2
(6.9)
2
where C RN N , C2 RN N , G1 RN N , G2 RN N , D1 RN N , D2 RN N ,
B RN 1 .
The QLDAE system is an equivalent representation of the original system, meaning that
the solution of the QLDAE system is exactly equal to that of the original one if the initial
conditions are consistent (to be defined in Section 6.3). Therefore, no approximation is
involved in this step.
QLMOR reduces this QLDAE system to a smaller model that is also in the QLDAE
form, as opposed to the polynomial form in prior methods. The important advantage here
is that higher-order terms beyond quadratic terms are absent in QLMOR-reduced model.
Therefore, it avoids the expensive computational cost (e.g., xr xr xr ) in a high-order
80
polynomial reduced model, and quadratic terms dominates the computational complexity.
The reduced order model is obtained by projection of the QLDAE system. While it looks
like a linear projection (of the QLDAE system), it indeed defines a nonlinear projection (of
the original system). To obtain the minimum subspace in the sense of moment matching,
we analyze the transfer functions of systems in the QLDAE form, and results related to
prior methods are derived. For example, it turns out that the moments of H2 (s1 , s2 ) of a
QLDAE system are a summation of the moments of H2 (s1 , s2 ) of a quadratic system and
a bilinear system direct application of existing codes are possible. In our algorithm, we
also avoid explicit moment evaluations, and therefore the algorithm has better numerical
stability properties.
Another point to note is that QLMOR serves as a new core method that can be combined
with any trajectory-based methods. For example, in TPWL method, to generate reduced
models along a trajectory, previous Krylov-subspace methods or TBR methods can be replaced with QLMOR straightforwardly, potentially leading to a more accurate model and
less number of regions along the trajectory. For another example, POD method can be directly applied to the transformed QLDAE system, which results in a reduced model whose
computational cost is lower than the usual POD reduced model if the reduced model size
remains the same.
To summarize, QLMOR is composed of three steps:
1. construct an equivalent polynomial system from the orignal system;
2. construct an equivalent quadratic-linear system from the polynomial system;
3. reduced the quadratic-linear system to a smaller quadratic-linear system.
We chart in Fig. 6.1 the basic flow of QLMOR, as well as the flow of MOR based on polynomial systems and bilinear systems for comparison. In what follows, Section 6.3 presents
the polynomialization and the quadratic-linearization procedure, and Section 6.4 discusses
the reduction step.
6.3
Quadratic-Linearization Procedure
In this section, we develop an algorithm to convert systems of ordinary differential equations to an equivalent system which has functions that are quadratic in state variables and
linear in inputs. We call it quadratic-linearization procedure, and we call the resulting
system QLDAEs (quadratic-linear differential-algebraic equations). We also discuss extensions to quadratic-linearize DAEs in the end of this section.
Considering nonlinear ODEs, the quadratic-linearization procedure is composed of two
steps:
81
Figure 6.1: QLMOR flow (yellow), with comparison to MOR based on polynomial form
(green) and bilinear form (gray).
82
1. convert the original nonlinear ODEs
d
x = f (x) + Bu
dt
(6.10)
6.3.1
To be more specific, we consider semi-general nonlinear ODEs (6.10) where the nonlinear functions f (x) = [f1 (x), , fN (x)] can be written as a linear combination of a set of
elementary functions {gi (x)}, i = 1, , k, i.e., the differential equation for xi is
d
xi = pTi x + ai,1 g1 (x) + ai,2 g2 (x) + + ai,k gk (x) + bi u, i = 1, , N,
dt
(6.11)
(6.12)
where P T = [p1 , , pN ], A = [ai,j ], g(x) = [g1 (x); , gk (x)] and B = [b1 ; ; bN ]. Note
that A is extremely sparse due to the fact that each differential equation typically contains
just a few nonlinear functions.
As mentioned in Section 2.7, this form of differential equations is prevalent in circuit
simulation, as well as other engineering problems. For example, in circuit MNA equations [108, 109], each gi (x) represents the current flowing into a node from a specific device;
in chemical rate equations, each gi (x) represents the reaction rate according to a specific
reaction; in mechanical applications, gi (x) is the force to the system.
The elementary functions gi (x) include, but not limited to, the ones listed in Table 6.2,
and compositions of these elementary functions. Because of the function composition, the
nonlinear functions in Table 6.2 in fact cover a very rich set of nonlinear functions that are
almost adequate in engineering problems.
Specifically, in circuit simulation, two important classes of models use functions that are
a subset of composition of these elementary functions.
The first class is the most commonly used transistor models, such as BSIM [118] or PSP
models [123]. These models are written in C/C++, and contain if-else statements to
implement piecewise-linear functions. For example, the C code
83
Table 6.2: Commonly used nonlinear functions
Function
Applications
1
N
x1 xN
chemical rate equations
MOSFET in saturation mode (x2 , x [122])
x
chemical rate equations, smoothing functions
k+x
x
e
diodes, BJTs, ion-channel models, smoothing functions
sin(x), cos(x) control systems (e.g., x is the angle to be steered)
if (y<1)
return f1(x,y);
else
return f2(x,y);
implements the function
f (x, y) = f1 (x, y)(1 s(y 1)) + f2 (x, y)s(y 1),
(6.13)
where s() is the step function, Because that the step function can be well-approximated and
smoothed by
s(t) 0.5 tanh(Kt) + 0.5,
(6.14)
where K is a large constant number, and that tanh(t) is a composition of elementary functions, we can represent all piecewise-linear functions as compositions of elementary functions
in Table 6.2.
The second class is models composed of look-up tables (LUT). The LUT model can be
viewed to implement a function whose outputs of given input values are stored in the LUT.
Therefore, one way to interpret the function that LUT implements is through the kernel
method (interpolation), i.e.,
X
f (x) =
f (xi )K(x xi ),
(6.15)
i
where K() is a kernel function [124]. By choosing K() appropriately, we can make f (x)
a piecewise-linear function, a polynomial function, or more generally the composition of
elementary functions.
Given these nonlinear ODEs, in the first step, we convert them into a polynomial system
which is defined as follows:
Definition 6.3.1 (Polynomial system of order M) A polynomial system (6.16) is said
to be of order M, if the highest order of monomial of x is M, i.e.,
84
M = maxk,i
X
P
ci,k xi =
N
l=1
X
k
i,k,l ,
l=1 i,k,l .
PN
X
k
u,
i = 1, , N.
(6.16)
In (6.16),
ci,k R, i,k R, RN , i,k,j 0, i,k,j 0, i,k,j Z, i,k,j Z, k, j.
(6.17)
Note that the nonlinear function on the RHS of (6.16) is linear in terms of the input
u which is an important property we will preserve in the process of polynomialization and
quadratic-linearization.
6.3.2
1
+ u.
k+x
(6.18)
1
,
k+x
and rewrite y =
1
k+x
as a
(6.19)
This approach preserves the linearity to the input u which is shown by the following
lemma:
Lemma 6.3.2 By adding polynomial algebraic equations through variable change, the resulting polynomial system is linear in terms of the input u.
85
Proof The new variable introduced is a function of x, i.e., yi = gi (x), and therefore, the
new polynomial equation introduced is independent of the input u.
In original differential equations, we only replace gi (x) with yi , and therefore, the ODEs
are still linear in terms of the input u.
The assumption for this approach to work is that the variable change can lead to a
polynomial algebraic equation. For example, any variable change to a rational function
leads to a polynomial algebraic equation. Constrained by this assumption, however, this
approach cannot deal with other nonlinear functions such as y = ex .
6.3.3
dgi
dx
consists of
(6.20)
i = 1, , k
i (x)
where gi (x) = dgdx
. Hence, if gi (x) consists of polynomial functions in x and yi , (6.20) is a
polynomial system.
We now examine some functions yi = gi (x) and show that in all cases we encounter, gi (x)
consists of polynomial functions in x and yi .
Considering uni-variable function (i.e., functions of only one variable), we have
y = ex (ex ) = ex = y
y = 1/(k + x) (1/(k + x)) = 1/(x + k)2 = y 2
y = x (x ) = x1 = yx1
y = ln(x) (ln(x)) = x1 .
2
(6.21)
The new equation is obtained by applying chain rule, and is the Lie derivative of yi with respect to
d
dt x.
86
Note that when y = x , y = yx1, where another new variable z = x1 can be further
introduced to eliminate the negative powers of x. The same reasoning applies to y = ln(x).
Therefore our claim that gi (~x) are polynomial functions in x and yi is correct for these
functions.
Some other uni-variable functions need to be handled by introducing two new variables.
For example, if gi (x) = sin(x), then we have
y1 = sin(x), y2 = cos(x)
(sin(x)) = cos(x) = y2
(cos(x)) = sin(x) = y1 .
(6.22)
(6.23)
which also shows that gi (x) are polynomial functions in x and yi is correct for these functions.
Furthermore, if g(x) is a composition of several functions, e.g., g(x) = (g2 g1 )(x) =
g1 (g2 (x)), we can convert it into polynomial nonlinearities by a similar procedure:
1. Introduce new variables y2 = g2 (x) and y1 = g1 (y2 );
2. Replace g1 (g2 (x)) by y1 in original equations;
3. Add y1 =
y1
y
y2 2
and y2 =
y2
x
x
2
1
and y
are polynomial functions of x, y1 , y2 , we conclude that
Since we can show that y
y2
~
x
we can polynomialize nonlinear ODEs which have nonlinear functions that are compositions
of nonlinear kernels in Table 6.2.
Similar to the previous approach, we have the following lemma regarding the linearity to
the input u:
Lemma 6.3.3 By adding polynomial differential equations through variable change, the resulting polynomial system is linear in terms of the input u.
Proof The new variable introduced is a function of x, i.e., yi = gi (x), and the new differential
equation is dtd yi = gi (x) dtd x.
Given that dtd ~x is a function that is linear in terms of the input u, and gi (x) is a function
of x and yi , we have dtd yi is also linear in terms of the input u.
Combining the above derivations, Theorem 6.3.4 follows:
87
Theorem 6.3.4 By iteratively applying polynomialization by adding polynomial algebraic
equations and taking Lie derivatives, a nonlinear system with nonlinear functions being compositions of functions in Table 6.2 can be polynomialized into a polynomial system in the
form of (6.16).
Moreover, the size of the equivalent polynomial system is linear with respect to the number
of elementary functions in original ODEs/DAEs.
Proof The proof for the first part is constructive the procedure described in previous
subsections (also summarized in Algorithm 6 in Section 6.3.4) converts ODEs/DAEs to an
equivalent polynomial system.
Regarding the size of the equivalent polynomial system, we note that for an elementary
function, we can always polynomialize it by adding constant number (O(1)) of new variables
and equations (normally 1 or 2); For a nonlinear function that is a composition of m elementary functions, i.e., g(x) = g1 (g2 ( gm (x) )), we can always polynomialize it by adding
O(m) new variables and equations. Therefore, the size of the equivalent polynomial system
is linear with respect to the number of elementary functions in original ODEs/DAEs.
Remark In circuit applications, the number of elementary functions is linear with respect
to the number of nonlinear devices in the circuit. Therefore, the polynomialization procedure
scales well with the size of the circuit. Similar results can be obtained for other engineering
problems.
6.3.4
Polynomialization Algorithm
Let the original nonlinear system be given in the form of (6.11), we summarize the polynomialization algorithm in Algorithm 63 . Because the size of the equivalent polynomial system
is linear with respect to the number of elementary functions in original ODEs/DAEs, the
computational complexity of Algorithm 6 is linear with respect to the number of elementary
functions in original ODEs/DAEs.
Example 6.3.5 Consider an ODE
x =
1
.
1 + ex
(6.24)
1
.
1 + y2
(6.25)
In this algorithm, we perform polynomialization by only taking Lie derivatives, because it is more general.
88
Algorithm 6 Polynomialization of Nonlinear ODEs
Inputs: E = [x1 , , xn ], the list of symbolic expressions of the ODEs.
Outputs: Yvar , the set of new variables;
Yexpr , the set of expressions of new variables;
E,
the
list
of
symbolic
expressions
of
the
polynomial
ODEs.
1: Initialize Yvar {}, Yexpr {};
2: repeat
3:
From E, the list of symbolic expressions of the nonlinear ODEs, pick a nonlinear
function g(x, v) that is not a polynomial function of x and v Yvar .
4:
Define a new state variable y = g(x, v).
5:
Add y into Yvar , and g(x, v) into Yexpr .
6:
Compute the symbolic expression of y = g (x, v)[x;
v].
7:
Add the symbolic expression of y to E.
8:
In E, replace occurrences of expressions in Yexpr by corresponding variables in Yvar .
9: until All expressions in E are polynomial functions of x and variables in Yvar .
By adding differential equations for y1 and y2 , we obtain a polynomial system
x = y1
y1 = y12 y2 = y2 y13
y2 = y2 x = y2 y1 .
(6.26)
Example 6.3.6 (A MOS Circuit) Consider a MOSFET circuit shown in Fig. 6.2, where
x1 and x2 are two node voltages, u1 and u2 are inputs, the resistance of resistors is 1 and
the capacitance of capacitors is 1. We can write the differential equations in terms of node
voltage x1 and x2 as follows:
x1 = x1 + u1
x2 = x2 + u2 IDS (x1 , x2 ),
(6.27)
2
(VGS VT )
if VDS VGS VT
(6.28)
2
IGS = 0,
where and VT are constants. The smoothed version of the Schichman-Hodges equations
89
u2
+
x2
x1
u1
+
VDS
VT )
VDS s(VGS VT VGS )
2
(6.29)
(6.30)
where K is a constant.
Therefore, the circuit equations (6.27) are
x1 = x1 + u1
h
e2K(x1 VT x2 )
x2 i
x2 2K(x1 V x2 )
x2 = x2 + u2 (x1 VT )
T
2
e
+1
1
.
+ (x1 VT )2 2K(x1 V x2 )
T
2
e
+1
(6.31)
90
To polynomialize (6.31), we introduce two new variables
y1 = e2K(x1 VT x2 )
1
,
y2 =
y1 + 1
(6.32)
and add differential equations for y1 and y2 . We then obtain a polynomial system
x1 = x1 + u1
h
x2 i
2
x2 y1 y2 + (x1 VT ) y2 .
x2 = x2 + u2 (x1 VT )
2
2
y1 =2Ke2K(x1 VT x2 ) (x1 x2 ) = 2Ky1 (x1 x2 )
1
y1 = y22 y1 .
y2 =
(1 + y1 )2
(6.33)
Remark The polynomialization of a nonlinear system is not unique. In fact, there should
exist a minimum-size polynomial system that corresponds to the original system. The
algorithm to find this minimum-size polynomial system could be devised by symbolic computational tools. Methods in logic synthesis might also be applicable.
Despite the non-uniqueness, however, as long as the existence and uniqueness conditions
of solutions to differential equations are satisfied, the polynomial system is equivalent to
the original system since they have the same solution.
Formally, suppose the original DAEs have state variables x. We define new variables
y = g(x), and expand the DAEs by adding differential equations for y. Let x(0) be the initial condition of the original DAEs, and (x(0), y(0)) be the initial condition of the expanded
DAEs. If y(0) = g(x(0)) and x(t) is the solution to the original DAEs, then (x(t), g(x(t))) is
the solution to the expanded DAEs. Moreover, if the expanded DAEs satisfy the existence
and uniqueness condition of solutions, (x(t), g(x(t))) is also the unique solution.
After the polynomialization, we obtain (6.16), a polynomial system of order M. Now we
present two approaches to convert a polynomial system (6.16) to a QLDAE system (6.8).
6.3.5
91
state variables up to that iteration. We then replace x1 1 xNN in original equations by y x,
and we add (6.34) into the original system. The equations in the form of (6.34) are then
re-written as quadratic algebraic equations, if possible, and are added to the original system.
Example 6.3.7 Consider an ODE
x = x3 + u(t).
(6.35)
(6.36)
Similar to the polynomialization procedure, quadratic-linearization procedure also preserves the linearity to the inputs u. Besides, we prove in the following theorem that a
polynomial system can always be quadratic-linearized by adding quadratic algebraic equations. As a byproduct, the theorem also provides an upper bound on the size of the QLDAE
system derived from a polynomial system.
Theorem 6.3.8 For a polynomial system of order M, the maximum number of state variables in the resulting QLDAE system does not exceed the number of state variables of
3 , , (xT )( M2 ) ]T , by adding quadratic algebraic equations,
2 , (xT )
[xT , (xT )
p ]T , p Z+ , then QLDAEs in y can generate any
2 , , (xT )
Proof Let y = [xT , (xT )
nonlinear functions of the form G1 y+G2 yy, which is essentially all the polynomial functions
in x of order not greater than 2p.
So, M 2p, i.e., p = M2 .
6.3.6
Similarly, taking Lie derivatives of polynomial functions preserves the linearity to the
input u, and it can also lead to a QLDAE system, as proved in Theorem 6.3.9. However,
the theoretical upper bound on the size of the resulting QLDAE system is more than that
by adding quadratic algebraic equations.
Theorem 6.3.9 For a polynomial system of order M, the maximum number of state variables in the resulting QLDAE system does not exceed the number of state variables of
3 , , (xT )M
2 , (xT )
]T , by adding differential equations that are derived by taking
[xT , (xT )
Lie derivatives of changed variables.
92
P
3 , , (xT )M
2 , (xT )
]T .
T
T
Proof Let y = xi11 xiNN , where N
k=1 ik M. Then, y [x , (x )
The differential equation of y is derived by taking its Lie derivative
d i1
(x1 xiNN ) = Lx (xi11 xiNN )
dt
N
X
d
ik 1
i1
xi xiNk .
x1 ik xi
=
dt
(6.37)
k=1
6.3.7
Quadratic-Linearization Algorithm
6.3.8
More generally, systems may be modeled by differential algebraic equations in the form
of
d
(q(x)) = f (x) + Bu.
(6.38)
dt
Similarly, we deal with equations where each element in f (x), denoted by fi (x), and each
element in q(x), denoted by qi (x), are summations of elementary functions, i.e.,
fi (x) =gi,1 (x) + + gi,ki (x),
qi (x) =hi,1 (x) + + hi,li (x),
(6.39)
where ki , li are constant integers denoting the number of nonlinear functions in fi and qi ,
respectively.
93
Algorithm 7 Quadratic-Linearization of Polynomial ODEs
Inputs: E = [x1 , , xn ], the list of symbolic expressions of the polynomial ODEs.
Outputs: Yvar , the set of new variables;
Yexpr , the set of expressions of new variables;
E,
the
list
of
symbolic
expressions
of
the
QLDAEs.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
Polynomialization
The polynomialization procedure is similar to the ODE case. We first introduce new
variables
yi,1 = gi,1 (x), , yi,ki = gi,ki (x),
zi,1 = hi,1 (x), , zi,li = hi,li (x).
(6.40)
i = 1, , N.
(6.41)
If equations in (6.40) can be written in the polynomial form, we can simply add them as
polynomial equations.
If equations in (6.40) cannot be written in the polynomial form, we add the differential
94
equations
d
d
d
d
yi,1 = gi,1
(x) x, , yi,ki = gi,k
(x) x,
i
dt
dt
dt
dt
d
d
d
d
zi,1 = hi,1 (x) x, , zi,li = hi,li (x) x.
dt
dt
dt
dt
(6.42)
However, we may not have explicit expressions for dtd x, and we do not reach a polynomial
system in the form of (6.16). We therefore extend the definition of a polynomial system to
be the following:
Definition 6.3.2 (Generalized polynomial system) The generalized polynomial system
is defined to be of the form (6.43), and it is said to be of order M, if the
highest order of
PN
PN
PN
monomial of x is M, i.e., M = maxk,i
j=1 i,k,l,j ,
l=1 i,k,l ,
l=1 i,k,l .
XX
k
X
k
i,k,lx1i,k,l,1
xNi,k,l,N
xl
X
k
(6.43)
(6.44)
Similar to the proof of Theorems 6.3.8 and 6.3.9, it is easy to show that by adding
quadratic algebraic equations and differential equations, we can convert a generalized polynomial system (6.43) to a generalized QLDAE system (6.44).
6.4
QLMOR
In this section, we focus on the reduction step for QLDAEs (6.8). Similar, but more
complicated derivations can be performed for generalized QLDAEs (6.44).
95
We start by performing a variational analysis to show the validity of applying Krylovsubspace methods. We then derive the moments of transfer functions of (6.8), which provides
a guideline of how to choose the proper Krylov subspace for generating a smallest reduced
model in the sense of moment-matching. The derivations in this section unifies previous
works on polynomial systems and bilinear systems [106, 107].
6.4.1
Variational Analysis
The variational analysis [105] starts by assuming the QLDAE system is a combination
of a series of homogeneous nonlinear sub-systems, whose responses to u(t) are x1 (t), x2 (t),
etc.. That is, when input is u(t), the response x(t) is
x(t) = x1 (t) + 2 x2 (t) + 3 x3 (t) + 4 x4 (t) + .
(6.45)
Therefore, by inserting (6.45) into (6.8), and equating the terms corresponding to i , i =
1, 2, , we obtain
d
x1
dt
d
C x2
dt
d
C x3
dt
d
C x4
dt
=G1 x1 + Bu,
=G1 x2 + G2 x1 x1 + D1 x1 u,
=G1 x3 + G2 (x1 x2 + x2 x1 ) + D1 x2 u + D2 (x1 x1 )u,
(6.46)
=G1 x4 + G2 (x1 x3 + x2 x2 + x3 x1 )
+ D1 x3 u + D2 (x1 x2 + x2 x1 )u.
From (6.46), if in the i-th set of differential equations (i.e., equations where LHS is dtd xi ),
we lump the terms without xi to be pseudo inputs, then the original system can be viewed as
a series of linear systems. Therefore, the Krylov subspace that guarantees that the projected
1
system matches the moments of these linear systems is K(G1
1 C, G1 ([B, G2 , D1 , D2 ])), if G1
is non-singular. However, the size of this subspace can easily exceed the size of the original
system. Therefore, care must be taken to choose a useful subspace.
6.4.2
Following [105], we can show that transfer functions for (6.8) can be written as follows.
(Only first three transfer functions are listed.)
H1 (s) = (sC G1 )1 B,
(6.47)
96
H2 (s1 , s2 ) =((s1 + s2 )C G1 )1
1
G2 (H1 (s1 ) H1 (s2 ) + H1 (s2 ) H1 (s1 ))
2
+ D1 (H1 (s1 ) + H2 (s2 )) ,
(6.48)
H3 (s1 , s2 , s3 )
=((s1 + s2 + s3 )C G1 )1
(
1
G2 H1 (s1 ) H2 (s2 , s3 ) + H1 (s2 ) H2 (s1 , s3 )
6
+ H1 (s3 ) H2 (s1 , s2 ) + H2 (s2 , s3 ) H1 (s1 )
(6.49)
Based on (6.47), (6.48) and (6.49), the Taylor expansion of these transfer functions can be
derived. Correspondingly, we obtain the moments of the transfer functions, i.e., coefficients
of sk for H1 (s), sk1 sl2 for H2 (s1 , s2 ), etc.. For example,
H1 (s) =
M1,k s =
k=0
H2 (s) =
Ak R1 sk
k=0
1
G2
2
+
X
Ak R1 sk1
k=0
Ak R1 sk2
k=0
+ D1
Ak R1 (s1 + s2 )k G1
1
k=0
k=0 l=0
"
X
k=0
Ak R1 sk2 +
Ak R1 sk2
k=0
k=0
Ak R1 sk1
Ak R1 sk2
k=0
!
#)
(6.50)
97
1
where A = G1
1 C, R1 = G1 B, and M1,k , M2,k,l are the moments for H1 (s) and H2 (s1 , s2 ),
respectively. Similar derivations can be performed for higher order moments.
An interesting point to note is that H2 (s1 , s2 ) is the same as the summation of the
H2 (s1 , s2 ) of a bilinear system C dtd x = G1 x + D1 xu + Bu and the H2 (s1 , s2 ) of a quadratic
system C dtd x = G1 x + G2 x x + Bu.
6.4.3
To illustrate how QLMOR generates the projection matrix, we derive the theorem and
algorithm for H2 (s1 , s2 ) in the following. These derivations can be easily extended for higherorder moment-matching.
From (6.50), we can derive the vectors that are the coefficients for an n-th order moment.
(E.g., for H2 (s1 , s2 ), the n-th order moments refers to the set of coefficients for sk1 sl2 , k+l = n.)
As an example, those vectors for H2 (s1 , s2 ) are listed in Table 6.3. We then have the following
theorem:
Theorem 6.4.1 Given a QLDAE system (6.8), suppose
[
[
j
R(V ) ={Ai R1 , i q} {Ai G1
D
A
R
,
i
+
j
q}
1
1
1
j
k
{Ai G1
1 G2 (A R1 ) (A R1 ), i + j + k q, k j}
(6.51)
Then if
x = V xr ,
Cr = V T CV,
Gr,1 = V T G1 V,
T
Dr,1 = V D1 V,
Br = V T B,
Gr,2 = V T G2 V V,
(6.52)
Dr,2 = V D2 V V,
the reduced system (6.53) matches the moments of H2 (s1 , s2 ) of the original system up to
q-th order.
Cr
d
xr = Gr,1 xr + Gr,2 xr xr + Dr,1 xr u + Dr,2 (xr xr )u + Br u
dt
(6.53)
98
Table 6.3: Moments of H2 (s1 , s2 )
Moments
Bilinear system Quadratic system
0th order
G1
G1
1 D1 R1
1 G2 R1 R1
1
1
1st order
AG1 D1 R1
AG1 G2 R1 R1
1
G1 D1 AR1
G1
1 G2 (AR1 ) R1
2 1
2nd order A2 G1
D
R
A
G
1 1
1
1 G2 R1 R1
1
1
AG1 D1 AR1
AG1 G2 (AR1 ) R1
1
2
2
G1 D1 A R1
G1
1 G2 (A R1 ) R1
1
G1 G2 (AR1 ) (AR1 )
3 1
3 1
3rd order A G1 D1 R1
A G1 G2 R1 R1
2 1
A G1 D1 AR1
A2 G1
1 G2 (AR1 ) R1
1
2
2
AG1
D
A
R
AG
1
1
1
1 G2 (A R1 ) R1
3
3
G1
G1
1 D1 A R1
1 G2 (A R1 ) R1
1
AG1 G2 (AR1 ) (AR1 )
2
G1
1 G2 (A R1 ) (AR1 )
A naive generation of the above subspace might lead to numerical stability problems
because of the explicit computation of Ai terms. For quadratic systems, [106] computes
the starting vectors directly and suggests to perform an orthogonalization between starting
vectors. However, that can still lead to numerical stability problems since the starting vectors
are computed explicitly.
Here we present a new algorithm to generate the basis for this subspace, as shown in
Algorithm 8. The elementary operation in the algorithm is to generate an orthonormal
projection matrix V = [v0 , , vq ] for a Krylov subspace Kq+1 (A, R). We employ Arnoldi
algorithm [126] in our implementation.
The computational complexity for this algorithm depends on the number of moments
to be matched. Suppose we want to match up to the p0 -th order moments of H1 , up to
the p1 -th order moments of the bilinear system, and up to the p2 -th order moments of the
quadratic system. From TheoremP
6.4.1, we know that the sizesP
of the subspaces for the three
2
1 +1
i) = O(p21) and O(1 + pi=1
i) = O(p22), respectively.
classes of moments are O(p0), O( pi=1
Therefore, the size of the final subspace is of size O(p0 + p21 + p22 ), and this is also the time
and space complexity of the Krylov basis generation algorithm, with the valid assumption
that system matrices are sparse. If we choose p0 = p1 = p2 = q, then the complexity of the
algorithm is O(q 2 ).
In this algorithm, we avoid the direct computation of Ai in the starting vector for each
1
1
i
Krylov subspace. We use G1
1 D1 vi instead of G1 D1 A R1 , and G1 G2 vi vj instead of
1
i
j
G1 G2 (A R1 ) (A R1 ), where vi is the i-th vector
the Krylov subspace
R1 ).
S i in
S K i(A,
1
1
i
j
We can prove that R1 (V ) = {A R1 , i q} {A GS
G1 G2 (Aj R1 )
1 D1 A R1 , i+j q} {A
S
1
(Ak R1 ), i+j+k q, k j} and R2 (V ) = {vi , i q} {Ai G1 D1 vj , i+j q} {Ai G1
1 G2 vj
vk , i + j + k q, k j} span the same subspace (proof by induction). Therefore, in the
99
computation of the starting vector, no Ai terms are involved.
We can also prove that vi1 vj1 is orthogonal to vi2 vj2 if i1 6= i2 or j1 6= j2 (by definition of
Kronecker product). Hence, the starting vectors have better numerical properties than [106].
1
However, it is not always the case that G1
1 G2 vi1 vj1 is orthogonal to G1 G2 vi2 vj2 , even
if vi1 vj1 is orthogonal to vi2 vj2 . Therefore, one further improvement of this algorithm
is to perform an orthonormalization of the starting vector to the basis already generated.
Proofs are omitted here due to page limits.
Algorithm 8 Generation of basis for subspace in Theorem 6.4.1
Inputs: C, G1, G2 , D1 , D2 matrices defining the QLDAEs;
q: the reduced model will match up to the q-th order moments of H1 and H2 .
Outputs: V : the projection matrix.
1
1: A = G1
1 C, R1 = G1 B.
2: Generate matrix V = [v0 , , vq ] for Kq+1 (A, R1 )
3: for i = 0 to q do
4:
Generate matrix W = [w0 , , wqi ]
for Kqi+1 (A, G1
1 D1 vi ). {subspace for moments of the bilinear system}
5:
V = qr([V, W ]).
6: end for
7: for i = 0 to q do
8:
for j = 0 to min(q i, i) do
9:
Generate matrix W = [w0 , , wqij ] for Kqij+1(A, G1
1 G2 vi vj ). {subspace
for moments of the quadratic system}
10:
V = qr([V, W ]).
11:
end for
12: end for
6.4.4
The loss of passivity in [106, 107, 121] are mainly caused by Taylor approximation. QLMOR alleviates this problem in the sense that the conversions to the polynomial system and
the quadratic-linear system are exact, and no passivity loss is introduced in that step.
In the reduction step, however, the passivity may again be lost. One may show that our
reduced model linearized at any xr,0 RNr , is a congruence transformation of the original
system linearized at x0 = V xr,0 RN . Under somewhat strict conditions such as the ones
in [11], the passivity of the linearized system can be preserved.
100
6.4.5
As mentioned in Section 6.3, different QLDAE formulations are all equivalent. However, the reduced QLDAEs are generally not equivalent. We discuss how different QLDAE
formulations lead to different reduced models.
Consider a nonlinear system,
d
x = f (x) + Bu,
dt
(6.54)
d
xq = G1 xq + G2 xq xq + D1 xq u + D2 xq xq u + Bu,
dt
(6.55)
where the state variables are xq = (x, y) where y = g(x) are nonlinear functions of x.
According to Theorem 6.4.1, by choosing appropriate projection matrices, we can match
the first few moments of (x, y). Therefore, by formulating different QLDAEs where different
ys are added, the moments of different ys are matched in different reduced models. This
result is of great importance in practice depending on the output signal, which is a nonlinear
function of state variables, we may compute reduced models that matches moments of these
output signals.
There is another subtle but important guideline of choosing new state variables only
add variables y = g(x) so that yDC = g(xDC ) = 0 where xDC is the DC solution of (6.54) for
u = uDC . (Without loss of generality, we assume that xDC = 0.)
To see why this is useful, we consider two QLDAEs that are equivalent to (6.54), and are
in terms of x and x , respectively:
d
x = G1, x + G2, x x + D1, x u + D2, x x u + B u,
dt
(6.56)
d
(6.57)
x = G1, x + G2, x x + D1, x u + D2, x x u + B u.
dt
The state variables x and x satisfy 1) x = 0 when u = uDC ; 2) x = x + x0 , where
x0 is a non-zero constant vector.
Generally, G1, 6= G1, , G2, 6= G2, , D1, 6= D1, , D2, 6= D2, , B 6= B . Therefore
the moments of two QLDAEs are not the same. This, however, does not imply that the
two QLDAEs are not equivalent. Rather, loosely speaking, it implies that the moments of
systems quadratic-linearized around different states are not the same.
To see that, we note that G1, is the Jacobian matrix of RHS of (6.56) computed at
x = 0 = x,DC , and G1, is the Jacobian matrix of RHS of (6.57) computed at x = 0.
However, when u = uDC , x = x +x0 = x0 6= 0. Therefore, in order to match moments of the
original system (6.54) around xDC , it is inappropriate to use G1, in moment computations,
101
and it is inappropriate to use (6.57) in QLMOR.
Rather, we can show that the moments of (6.56) linearized around 0 match the moments of
(6.54) linearized around xDC , thus leading to the guideline of choosing y such that yDC = 0.
The proof of this result involves a bit linear algebra, but is quite straightforward, and is
omitted here.
6.4.6
6.4.7
Multi-Point Expansion
The multi-point Krylov subspace method [127] can be directly applied in QLMOR, and
it potentially leads to a smaller and more accurate model for a specified frequency range of
interest.
Also note that the quadratic-linearization procedure in Section 6.3 might render a QLDAE system where G1 is singular. This can be a potential problem for generating Krylov
subspaces at s = 0. The current workaround for this problem is to generate Krylov subspaces
at s = s0 , s0 0.
102
6.5
In this section, we illustrate the applicability of the QLMOR method presented in the
previous sections, demonstrating its efficiency and accuracy on example circuits and systems.
6.5.1
A System with a
x
1+x
Nonlinearity
x
are common nonlinearities in biochemical rate
As noted in Section 6.3, x, x x and K+x
equations. To test how QLMOR works, a random second-order polynomial system with a
x
function was generated to mimic a biochemical reaction system4 :
1+x
10x1
d
+ Bu = 0.
x + G1 x + G2 x x + e1
dt
1 + x1
(6.58)
X
1
=
(x)i = 1 x + x2 x3 + x4 x5 + , (|x| < 1).
1+x
i=0
(6.59)
Notice that (6.59) converges only when |x| < 1 (also clearly seen in Fig. 6.3), making the
Taylor approximation irrelevant and highly inaccurate for |x| 1; as a result, any reduced
model derived from this approximation cannot be accurate for |x| 1. Indeed, the Taylor
approximation model turns out to be an unstable system for polynomial orders from 2 to
7, i.e., the polynomial model has a finite escape time [73]. As mentioned in [107], such
approximations do not preserve properties such as stability and passivity, except by chance.
To apply the QLMOR method, we perform quadratic-linearization of the original system
x1
by making the variable change y = 1+x
. Following the procedure in Section 6.3, we obtain
1
a QLDAE system equivalent to (6.58):
d
x + G1 x + G2 x x + e1 10y + Bu = 0,
dt
x1 + y + x1 y = 0.
4
(6.60)
The eigenvalues of G1 were manipulated to be positive such that the linearized system at 0 is stable. The
eigenvalues are uniformly sampled on [0, 20]. G2 was made sparse (5% non-zeros), which is approximately
x1
were made
of the same order in real rate equations and circuit equations. The coefficient in front of 1+x
1
large enough to make the system exhibit high nonlinearity.
103
1
50
0.9
40
1/(1+x)
1storder approximation
2ndorder approximation
3rdorder approximation
4thorder approximation
5thorder approximation
6thorder approximation
0.8
30
0.7
20
0.6
0.5
10
0.4
1/(1+x)
1storder approximation
2ndorder approximation
3rdorder approximation
4thorder approximation
5thorder approximation
6thorder approximation
0.3
0.2
0.1
0
0
0.1
0.2
0.3
0.4
0
10
20
0.5
x
(a)
0.6
0.7
0.8
0.9
30
0
0.2
0.4
0.6
0.8
1
x
1.2
1.4
1.6
1.8
(b)
Figure 6.3: Polynomial approximations of 1/(1 + x). Fig. 6.3(a) shows the interval x [0, 1].
Fig. 6.3(b) shows the interval x [0, 2].
Notice that the size of (6.60) is 11, only one more than the size of the original system.
Applying QLMOR on (6.60), a size 6 reduced model is generated that matches 4 moments
of H1 (s) and 2 moments of H2 (s1 , s2 ). We then apply a small and fast input signal
u(t) = cos(4t) so that small higher-order harmonic distortions are generated in the output
waveform, as plotted in Fig. 6.4(a). To assess the accuracy of the reduced order model, we
also plot the first few Fourier coefficients of x1 (t) in Fig. 6.4(b). (Because x1 (t) exhibits the
largest nonlinearities compared to the other state variables, we focus on its waveforms to
evaluate the accuracy of the reduced model.) We also excite this system by a large and
slow input signal u(t) = 10 cos(2t), so that large higher-order harmonic distortions are
generated in the output waveform, as plotted in Fig. 6.4(c). The speedup of the transient
simulation of the reduced model over the original model is about 1.5 although the model
2 matrix becomes dense, and therefore, the speedup is
size is nearly halved, the reduced G
not huge. However, considering that the original model size is very small, a speedup of 1.5
is, in fact, not inconsiderable.
We further manipulate the system (6.58) to make it more nonlinear, and see how QLMOR
behaves. In this case, we simply change the nonlinear term to be a function of x2 5 , i.e.,
d
10x2
+ Bu = 0.
x + G1 x + G2 x x + e1
dt
1 + x2
5
(6.61)
We observe in simulations that the outputs of this system exhibit larger high-order harmonic distortions
than (6.60).
full model x
0.15
full model x
reduced model x
0.1
reduced model x
0.05
0
0.05
0.1
0.5
1
time
1.5
104
40
60
80
100
0
reduced model x
reduced model x
1.5
1
0.5
0
0.5
1
0.5
1.5
2
time
2.5
3.5
full model x
4
6
harmonics
10
full model x1
2.5
reduced model x1
20
full model x1
full model x1
10
reduced model x
20
30
40
50
60
70
80
0
4
6
harmonics
10
Figure 6.4: Fig. 6.4(a), Fig. 6.4(b) show the time-domain waveforms of x1 , x2 and Fourier
coefficients of x1 , respectively, when u(t) = cos(4t). Fig. 6.4(c), Fig. 6.4(d) show the
time-domain waveforms of x1 , x2 and Fourier coefficients of x1 , respectively, when u(t) =
10 cos(2t).
105
The same model reduction procedure is repeated. Results are plotted in Fig. 6.5 for inputs
u(t) = 8 cos(3t) and u(t) = 10 cos(3t). In this case, strong nonlinearities are present in the
waveforms, but the QLMOR model still produces quite accurate results. An interesting point
to note is that at time t 0.5, the reduced model does not quite capture a fast transient
well. This is only to be expected, however, since MOR based on moment matching tends
to reduce size by eliminating fast components in the system, while maintaining fidelity for
slower components: observe that once the fast transient dies out, the QLMOR model is far
more accurate. This inaccuracy at higher frequencies is also apparent from the harmonics
shown in Fig. 6.5(d).
6.5.2
We consider the nonlinear transmission line circuit shown in Fig. 6.6 [2]. All resistors and
capacitors are set to 1 and the diode I-V characteristic is iD = e40vD 1. The input is set to
the current source i = u(t); the output is the voltage at node 1.
The modified nodal equations for this circuit are
v1 = 2v1 + v2 + 2 e40v1 e40(v1 v2 ) + u(t),
vN = vN + vN 1 1 + e40(vN1 vN ) .
2 i N 1,
(6.62)
(6.63)
full model x
full model x2
reduced model x
1.5
reduced model x
1
0.5
0
0.5
0.5
1.5
time
2.5
106
40
60
80
100
0
full model x
10
reduced model x
2
1
reduced model x
8
6
4
2
0
2
2.5
full model x
1.5
time
4
6
harmonics
10
12
reduced model x1
20
0.5
full model x
0
10
full model x
20
reduced model x1
30
40
50
60
70
80
0
4
6
harmonics
10
Figure 6.5: Fig. 6.5(a), Fig. 6.5(b) show the time-domain waveforms of x1 , x2 and Fourier
coefficients of x1 , respectively, when u(t) = 8 cos(3t). Fig. 6.5(c), Fig. 6.5(d) show the
time-domain waveforms of x1 , x2 and Fourier coefficients of x1 , respectively, when u(t) =
10 cos(3t).
107
N-1
......
(6.64)
Therefore, (6.63) and (6.64) together are already in QLDAE form with size 2N. This size
is far less than the system size that would result from a second-order Carleman bilinearization
[107] (i.e., N + N 2 ). This example illustrates that if we use branch voltages as state
variables, a much smaller QLDAE system can result. Intuitively speaking, this is reasonable
since the device models are usually defined in terms of branch voltages and branch
currents. This suggests that sparse tableau like circuit equation formulations [128], where
the state variables are branch voltages and currents, can be useful for quadratic-linearization.
We apply QLMOR on this nonlinear transmission line circuit with N = 10 and compare
against weakly-nonlinear Taylor-expansion-based MOR methods, setting the reduced size to
7.6 To drive the circuit sufficiently hard to produce large harmonic distortions, we apply
u(t) = 1.8 sin(2 0.1t) and plot waveforms in Fig. 6.7. As expected, the Taylor-based
weakly nonlinear reduced model matches the Taylor approximation of the original system
well; however, the original Taylor approximation itself does of course not match the full
model, because of its validity is limited to small inputs. QLMOR, on the other hand,
produces a reduced model which matches the original system well, the relative error being
about 1%.
To illustrate that QLMOR scales to larger system sizes, we apply QLMOR to nonlinear
transmission line circuits of 50,100,150,200 stages. To make a fair comparison, we tune
QLMOR-reduced model sizes such that in several simulations, the ratio of the maximum
6
The size q = 7 is obtained by trial-and-error. With q larger than 7, QLMOR gives almost the same
result as that of q = 7.
108
0.15
full model
reduced quadraticlinear model
reduced 3rdorder polynomial model
3rdorder polynomial model
0.1
0.05
0
0.05
0.1
0.15
0.2
0
10
15
time
20
25
30
Figure 6.7: Time-domain waveforms (x1 ) of full model, quadratic-linear reduced model, 3rdorder polynomial reduced model and 3rd-order polynomial model. (u(t) = 1.8 sin(2 0.1t))
109
0.1
full model x
0.05
0
0.05
0.1
0.15
0
10
15
time
20
25
30
absolute error of the time-domain waveform to its amplitude is within 1%. As a result,
we produce reduced models of size 11, 15, 18, 20 corresponding to 50,100,150,200 stages,
respectively. Again, we excite large harmonic distortions in the simulation, and a typical
transient waveform is shown in Fig. 6.8 (for the circuit with 50 stages).
As comparison, Taylor-based weakly nonlinear MOR did not succeed in producing a
reduced model on account of shortage of memory for computing V T G3 (V V V ). For
3
example, when N = 50, q = 11, V R5011 , we have V V V R(5011) . Hence, suppose
the size of a double-precision number is 8 bytes, V V V requires a minimum memory of
(50 11)3 8 Bytes, which is 1.2396 GB.
0
full model x1
reduced quadraticlinear model x1
20
40
60
80
100
0
4
6
harmonics
10
Figure 6.8: Fig. 6.8(a), Fig. 6.8(b) show the time-domain waveforms and Fourier coefficients
of x1 , respectively, when u(t) = 1.8 cos(2 0.1t). The size of the full model is 50; the size
of the reduced model is 11.
110
Chapter 7
NTIM: Phase/Timing Macromodeling
In this chapter, we present a special macromodeling technique that captures phase responses of nonlinear systems. While the phase response is yet to be defined, it can be
intuitively thought of as timing behaviors such as delay or phase modulation, and therefore
is an extremely important behavior in many practical systems, especially in circuits.
We formally define the phase response and discuss its applications. We derive a phase
macromodel via a nonlinear perturbation analysis, and present numerical methods to compute the phase macromodel. We show its applications in modeling circuits and neurons, and
compare the results with the full model in terms of accuracy and computational cost.
7.1
In many systems, due to design intents or physical constraints, the waveforms of state
variables are almost the same, and therefore can be characterized by the phase information.
For example: in phase modulation circuits, we aim to design a system that transforms
the input signal u(t) and a carrier signal x(t) to a output signal y(t) = x(t + u(t)). In
digital circuits, outputs of waveforms typically saturate at the power supply voltages, and
the waveforms all look like a ramp waveform, either from 0 to 1 or from 1 to 0.
However, the exact shape of waveforms are not the same under different input signals, and
these shapes may be characterized by the point-to-point delay information, or equivalently
the phase information.
In fact, using phase information alone, we may represent a very rich set of waveforms.
Example 7.1.1 Suppose x(t) = sin(2t), and
y(t) = x(t + z(t)) = sin(2(t + z(t)))
(7.1)
By choosing z(t) carefully, y(t) can represent a very rich set of waveforms. As shown in
Fig. 7.1,
111
1. z1 (t) = 1/4, y1 (t) = sin(2(t 1/4)) is x(t) delayed by 1/4;
2. z2 (t) = t, y2 (t) = sin(4t) doubles the frequency of x(t);
1
3. z3 (t) = t + 2
sin1 (2s(t 0.5) 1) (s(t) is a step function) , y3 (t) = 2s(t 0.5) 1
is a square waveform;
Generally, given x(t) and y(t) with maxt x(t) maxt y(t) and mint x(t) mint y(t), we
can always find z(t) such that y(t) = x(t + z(t)). In the case where x(t) = sin t, we have
z(t) = sin1 y(t) t. More generally, if x(t) is a vector, representing a trajectory in the state
space, then y(t) = x(t + z(t)) can be used to represent any other waveforms that have the
same trajectory in the state space.
While the phase information alone is powerful in representing waveforms, it is not capable
of representing any waveform whose trajectory deviates from the original trajectory. In that
case, if the original waveform is xS (t), and the new waveform is xp (t), then
xp (t) = xS (t + z(t)) + y(t),
(7.2)
where y(t) is not identically equal to 0. However, we may choose z(t) such that xS (t + z(t))
approximates xp (t) in some sense. In the setting of differential equations, we will define
such a z(t) to be the phase response. Formally, we have the following definition.
Definition 7.1.1 (Phase response) Consider a nonlinear dynamical system described by
a set of differential algebraic equations
d
q(x(t)) + f (x(t)) + b(t) = 0.
dt
(7.3)
Suppose the (unperturbed) solution of (7.3) to (unperturbed) input bS (t) is xS (t), and the
(perturbed) solution of (7.3) to (perturbed) input bp (t) is xp (t). There exists z(t) : R+ R,
such that xS (t + z(t)) best approximates xp (t) in some sense1 .
We call z the phase variable, and z(t) the phase response.
Notice that xS (t) may be any waveform in this definition. Indeed, we may derive bS (t)
by plug xS (t) in (7.3), i.e.,
d
bS (t) =
q(xS (t)) + f (xS (t)) .
(7.4)
dt
1
112
z (t)
z (t)
0
0.5
1
2
0
0.1
0.2
0.3
0.4
0.5
t
0.6
0.7
0.8
0.9
0
0
0.1
0.2
0.3
0.4
0.5
t
y (t)
1
0.5
0.5
0.5
0.5
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
z (t)
0.5
0.5
0.5
0.3
0.4
0.5
t
y (t)
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0.5
0.5
0.5
0.1
0.2
0.3
0.4
0.6
0.7
0.8
0.9
0.5
t
y (t)
0.6
0.7
0.8
0.9
0.5
0.6
0.7
0.8
0.9
1
0
0.2
0.9
z (t)
0.1
0.8
0.5
1.5
0
0.7
1
0
0.6
y (t)
0.5
0.6
0.7
0.8
0.9
1
0
0.1
0.2
0.3
0.4
0.5
0
0.5
1
1.5
0
0.1
0.2
0.3
0.4
0.5
t
0.6
0.7
0.8
0.9
0.6
0.7
0.8
0.9
y (t)
5
1
0.5
0
0.5
1
0
0.1
0.2
0.3
0.4
0.5
113
This definition can also be understood from a geometrical view point. As shown in
Fig. 7.22 , the solid (blue) orbit is the trajectory to input bS (t), and the dotted (red) orbit is
the trajectory to input bp (t). At time t, on the orbit of xS (t), there exists a point given by
xS (t + z(t)) which best approximates xp (t) on the perturbed orbit. The remaining difference
y(t) = xp (t) xS (t + z(t)) is called the amplitude deviation [29]. Therefore, we can
view the response to input bp (t) as combined effects of the phase response and the amplitude
deviation, i.e.,
xp (t) = xS (t + z(t)) + y(t).
(7.5)
Remark This definition is dependent on the trajectory of xS (t), and therefore bS (t). Therefore, it may be called a xS -PR (phase response).
With this definition of phase response, it is worth noting that although the phase response
is an important characteristic of oscillators as pioneered by many researchers [29, 129], it is
2
We show two closed orbits in Fig. 7.2 for simplicity, but there is no restriction of the orbit being closed
in our definition. The orbit can be arbitrary.
114
not a special characteristic only for oscillators, but generally for all kinds of dynamical
systems. It has broader applications in many fields.
Phase response is of great importance in circuits in particular. For examples, in oscillators, phase noise or timing jitter is the phase response in a stochastic setting where
noise/uncertainties are present; in digital gates and RC interconnects, timing information
such as delay and slope is essentially derived from the phase response; in a pulse generator,
the pulse width can also be derived from the phase response.
Scientists in biology and neuroscience also place a great emphasis on quantitative understanding of the phase response. For examples, normal human circadian rhythm disorders,
which can be caused by insomnia, fatigue and jet lag, are phase responses of underlying
physiological processes under different perturbations; synchronization of a neural network is
a phenomenon where the phase responses of all neurons settle down to the same value; some
neuron firing/bursting behaviors, which are exhibited when currents are injected, can also
be understood as the phase response of the neuron.
7.2
Phase Macromodels
With the above definition of the phase response, the phase macromodel is a model that
describes the dynamics of the phase response. In another word, the phase macromodel is
composed of a differential equation in terms of z(t), and therefore is a 1-dimensional reduced
order model.
To use the phase macromodel, we simulation the differential equation for z(t), instead of
(7.3), and use xS (t + z(t)) to approximate the original solution xp (t). With the assumption
that y(t) is small, xS (t + z(t)) will give a very good approximation to xp (t).
In many systems, the validity of this assumption depends on both the system properties
as well as input perturbations. Fortunately, we see a lot of applications that this assumption
is satisfied. For examples: in oscillators that are asymptotically orbital stable, the periodic
orbit is an attractor which by definition attracts all trajectories nearby; in digital applications, although circuits never reside in a periodic steady state, their trajectories follow
almost the same path since inputs are basically signals switching between 0 and 1 with
variations of the slope, to the first-order approximation.
Among all possible perturbations, there is also a special one that is very useful and
normally makes the above assumption satisfied. This is the phase perturbation, i.e., the
perturbed input bp (t) = bS (t + p(t)) is only a phase-modulated version of bS (t) by some
signal p(t). For example, it can be used to represent an FM/PM signal in RF applications,
or a delayed ramp signal with any slope in digital applications. In particular, when p(t)
is a constant, the perturbed solution xp (t) lies exactly on the same orbit of xS (t) if initial
condition is set properly.
115
7.3
(7.7)
=
=
n
X
i=1
n
X
i=1
n
X
ui (t)
i=1
(7.8)
(7.9)
Z t
i t T
i (ts) T
ui (t) e vi (0)w0 +
e
vi (s)bw (s)ds ,
0
where w0 is the initial condition for (7.8), (t, s) = U(t)D(t s)V T (s) is the state transition matrix of (7.8), i s are the Floquet exponents, D(t) = diag(e1 t , , en t ), U(t) =
[u1 (t), , un (t)], V (t) = [v1 (t), , vn (t)], and V T (t)U(t) = In . More theories and proofs
about LPTV systems can be found in [131].
S (t)
We now introduce a lemma showing that dxdt
, the time derivative of the periodic solution
116
xS (t) of (7.6), satisfies an LPTV system.
Lemma 7.3.1 The time derivative of the periodic solution xS (t) of (1.2), i.e.,
satisfies
dbS (t)
d
w(t) = G(t)w(t) +
,
dt
dt
and can be written as
N
X
dxS
(t) = U(t)c(t) =
ui (t)ci (t),
dt
i=1
d
(xS (t))
dt
(7.10)
(7.11)
where
Z r
dxS
d
i r T
i (rs) T
e vi (0)
ci (t) =
lim
(0) +
bS (s) ds .
e
vi (s)
r=t+kT,k
dt
ds
0
(7.12)
(7.13)
(7.14)
7.4
With the important assumption that the trajectory of the perturbed system stays close
to the trajectory of xS (t), the key idea is to show that under the perturbed input bp (t), the
perturbed response xp (t) can be decomposed into the phase response z(t) and the amplitude
deviation y(t) in a reasonable way, i.e.,
xp (t) = xS (t + z(t)) + y(t),
(7.15)
and that by defining the right differential equation for z(t), y(t) is minimized in some sense.
117
To show this, we start by defining the phase equation, i.e., the differential equation for
the phase response z(t). We then show that the input bp (t) can be decomposed into bz (t)
and by (t) such that when only bz (t) is applied to (7.6), the perturbed response is exactly
xS (t + z(t)). We then derive the first-order approximation of y(t) by linearizing original
differential equations around the phase-shifted solution xS (t + z(t)), and show that certain
decomposition of y(t) is minimized.
Definition 7.4.1 (Phase Equation) We define the phase equation to be
cT (t + z)c(t + z)
dz
= cT (t + z)V T (t + z) [bp (t) bS (t + z)] ,
dt
(7.16)
where c(t) = [c1 (t), , cN (t)] is defined in (7.12) and V (t) is defined in (7.9).
With the definition of z(t), we present a theorem showing that part of the input bp (t)
contributes only to the phase response.
Theorem 7.4.1 Given any perturbed input ~bp (t), define
cT (t + z)
V T (t + z)
cT (t + z)c(t + z)
[bp (t) bS (t + z)] U(t + z)c(t + z),
bz (t) =bS (t + z) +
(7.17)
(7.18)
(7.19)
f (xS ( )) + bS ( ) + U( )c( )
(7.20)
118
Substituting (7.17) in (7.20), we have
T
dz
c ( )
T
U( )c( )
= T
V ( )(bp (t) bS ( )) U( )c( ).
dt
c ( )c( )
(7.21)
Since U( )c( ) is dxdtS ( ) whose elements are not simultaneously zero, we obtain (7.16). Since
all above derivations are reversible, we have also shown that xS (t + z(t)) solves (7.18).
It remains to show that by decomposing the perturbed response xp (t) into phase response
and amplitude deviation according to (7.16), the amplitude deviation is minimized in some
sense. This is proven by the following theorem. The proof also gives another derivation of
the phase equation (7.16).
Theorem 7.4.2 Suppose the perturbed response is xp (t) = xS (t + z(t)) + y(t), then to the
first-order approximation, y(t) is
Z
y(t) = U( )
e(ts) r(s)ds,
(7.22)
0
dz
dt
(7.23)
= a(t), then
bz (t) = bS (t + z) + a(t)
dxS
(t + z).
dt
(7.24)
(7.25)
Therefore, we have
d
y(t) = G(t + z)y(t) + [bp (t) bz (t)] .
dt
(7.26)
119
With the assumption that
bp (t) bz (t), and therefore
dz
dt
d
y(t) = G(t)
y (t) + bw (t).
dt
Therefore, applying (7.9), we have
y(t) = y(t) = U(t)
(7.27)
(7.28)
(7.29)
The minimization of ||r(t)||2 boils down to the problem of minimizing ||Ax b||2 where
A = c(t + z), x = a(t) and b = V T (t + z) [bp (t) bS (t + z)]. The solution to this problem is
simply
cT (t + z)
a(t) = T
V T (t + z)(bp (t) bS (t + z)),
(7.30)
c (t + z)c(t + z)
which again leads to (7.16).
7.4.1
The model shown in Section 7.4 can be interpreted as a special reduced order model in
the projection framework mentioned in Chapter 4, where state variables are forced to lie
in an Nr -dimensional manifold defined by x = v(xr ), and the residual is forced to lie in a
Nr -dimensional manifold defined by rz = w(rx ).
To derive the same phase macromodel under the projection framework, we first enforce
that the perturbed response xp (t) only induces a phase shift, i.e., the state variable lie on the
same orbit as that of the unperturbed solution. This is equivalent to saying that we project
the state space onto the one-dimensional manifold defined by the unperturbed solution xS (t)
which has a natural parameterization by time t. Therefore, the projection function is defined
by
x(t) = v(z(t)) = xS (t + z).
(7.31)
Substituting (7.31) in (7.6), we have
d
xS (t + z) =f (xS (t + z)) + bp (t).
dt
(7.32)
120
Therefore,
d
dxS
d
dxS
dz
[xS (t + z)] =
( )
=
( ) 1 +
dt
d
dt
d
dt
dxS
dz
dxS
( ) +
( ) .
=
d
d
dt
(7.33)
(7.34)
(7.35)
(7.36)
or equivalently,
dz
= bp (t) bS ( ).
(7.37)
dt
Therefore, after the projection onto the periodic orbit, we obtain N equations (7.37) and
1 state variable z. To derive a reduced order model, we further project the residual, i.e.,
(7.37), to a one-dimensional manifold. In our model, we have chosen the projection to be
U( )c( )
w(rx ) = cT ( )V T ( )rx ,
(7.38)
which is shown to minimize the amplitude deviation in the sense of Theorem 7.4.2.
7.4.2
The PPV phase macromodel [29] is specifically devised to characterize the phase response
of an oscillatory system. Its limitation is that it is only applicable to oscillators. This
excludes it from wide applications such as traditional timing analysis and excitory firing
neuron network simulations.
Following previous derivations, the PPV macromodel is indeed a special case of our
general-purpose phase macromodel. We can obtain the PPV model for oscillators by choosing
c(t) in (7.16) to be [1, 0, , 0]T .
This choice of c(t) in fact is the definition in (7.12) for autonomous oscillators. To see
that, we note that for autonomous oscillators, the following two conditions are satisfied: (1)
121
one Floquet exponent of (7.10) to be 0 and others negative (e.g., 1 = 0, i < 0, 2 i n),
and (2) bS (t) is constant. Therefore, (7.12) becomes
d~xS
(0) = 1,
t
dt
d~xS
(0) = 0,
ci (t) = lim ei t~viT (0)
t
dt
(7.39)
2 i n.
7.4.3
Previous derivations are based on the assumption that xS (t) is the asymptotic periodic
orbit in the state space. This is not generally the case for practical systems such as digital
circuits. However, if we force the input to be a periodic signal (for example, a square
wave), then the system will converge to a periodic steady state, and transient behaviors (for
example, the rising and falling edge of waveforms) are well-approximated by its periodic
steady state. Hence, we can force the input to be periodic as long as the periodic trajectory
well-approximates the non-periodic ones. As a result, the derived phase model will also
characterize the transient phase responses.
7.5
Algorithm
Based on the above analysis, the algorithm to generate a general-purpose phase macromodel is shown as follows.
Algorithm 9 General-Purpose Phase Macromodeling
1: Given input bS (t), compute the periodic steady state xS (t).
dxS (t)
, the time derivative of xS (t).
2: Compute dt
3: Perform Floquet decomposition of the LPTV system
dw
= G(xS (t))w(t) (compute U(t), V (t), and D(t)).
dt
dxS (t)
4: Compute c(t) = V T (t) dt
.
5: Precompute cT (t)c(t) and cT (t)V T (t) in (7.16) and store them in the model.
7.6
As shown in Algorithm 9, to compute the phase macromodel, we need to compute the peS (t)
, and perform a full Floquet decomposition
riodic steady state xS (t), its time derivative dxdt
of an LPTV system.
122
PSS analysis has been widely adopted in commercial tools and usually uses shooting
method [132], FDTD method [132] or harmonic balance [132]. The time derivative of the
PSS solution can be computed either in the time domain directly by finite difference approximation or in the frequency domain and be transformed back to time domain.
The main numerical challenge lies in the Floquet decomposition, which is numerically
much more ill-conditioned. Implementations according to textbook definitions using Monodromy matrix fail without question. In this section, we first demonstrate that methods
based on explicitly computing the Monodromy matrix are numerically extremely bad. We
then present3 two much better-behaved methods to solve this challenge, one based on harmonic balance, and the other based on finite-difference time-domain (FDTD) method.
7.6.1
(7.40)
Floquet decomposition [131] refers to computing matrices U(t), V (t) and D = et (where
is a diagonal matrix of all Floquet exponents) such that the state transition matrix of
(7.40) is
(t, s) = U(t)D(t s)V T (s).
(7.41)
The definition for Floquet exponents are derived from the eigenvalues of the monodromy
matrix defined by
B = (T, 0) = P eT P 1,
(7.42)
and U(t) is then defined as (t, 0)P et , and V (t) = U 1 (t).
This definition naturally gives a straightforward implementation of Floquet decomposition:
Algorithm 10 Floquet Decomposition using Monodromy Matrix
1: Starting with initial condition X(0) = I, integrate (7.40) to compute (t, 0) for t [0, T ].
2: Compute the monodromy matrix B = (T, 0).
3: Eigen-decompose the monodromy matrix B = P eT P 1 .
4: Compute U(t) = (t, 0)P et and V (t) = U 1 (t).
This algorithm will work in perfect precision. However, it almost never works in practice.
For real systems, most Floquet exponents are normally negative. This means that using any
initial condition P , the solution to (7.40) (which is (t, 0)P ) tends to 0 exponentially. This
3
While these methods seem to have been developed before, we present our derivation and also explain
the details in the derivation.
123
makes the computed monodromy matrix B extremely numerically singular, and the eigendecomposition of B provides useless results. Also, computing V (t) using the inverse of U(t)
is not a good choice since the smoothness of V (t) can easily be destroyed.
7.6.2
As the first and the most important step to tackle the problems in the Monodromy matrix
method, we derive equations whose unknowns are U(t) and V (t), respectively. Because U(t)
and V (t) are periodic, the resulting differential equations have solutions that are periodic,
instead of decaying to 0. Therefore, this avoids numerical errors in the Monodromy matrix
computation. Moreover, having equations for U(t) and V (t) separately avoids computing
V (t) = U 1 (t).
To derive these equations, we start with a fundamental matrix
X(t) = (t, 0) = U(t)D(t)V T (0),
(7.43)
d
U(t)D(t)V T (0) + G(t)U(t)D(t)V T (0) = 0
dt
d
[U(t)D(t)] + G(t)U(t)D(t) = 0
dt
d
d t
t
+ G(t)U(t)et = 0.
U(t) e + U(t)
e
dt
dt
(7.44)
Therefore,
d
U(t) + U(t) + G(t)U(t) = 0.
(7.45)
dt
Similarly, to derive the equations for V (t), we apply the same technique to the adjoint
system of (7.40)
dy(t)
GT (t)y(t) = 0,
(7.46)
dt
which has a fundamental matrix
Y (t) = V (t)D(t)U T (0).
(7.47)
GT (t)V (t)e(t) = 0.
V (t)e
+ V (t)()e
dt
(7.48)
124
Therefore,
d
V (t) V (t) GT (t)V (t) = 0.
(7.49)
dt
Therefore, we have shown that U(t) and V (t) satisfy differential equations (7.45) and
(7.49), respectively. To solve for U(t) and V (t), we note that they are by definition periodic
with period T , and therefore we need to force the their periodicity during the computation.
Two numerical methods that force the periodicity are HB and FDTD, and we now show how
to apply these methods to solve for U(t) and V (t).
7.6.3
Notations
We follow the notations in [133]:
periodic matrix or vector A(t) which has a
PM For any
jit
truncated Fourier series A(t) = i=M Ai e
( = 2
), we define the block vector VA(t) of
T
Ai s and the block-Toeplitz matrix TA(t) of Ai s to be
VA(t)
..
.
A2
A1
A0
A1
A2
..
.
, TA(t) =
..
..
..
..
.
.
.
.
A0 A1 A2 A3
A1 A0 A1 A2
A2 A1 A0 A1
A3 A2 A1 A0
A4 A3 A2 A1
..
..
..
..
.
.
.
.
..
.
A4
A3
A2
A1
A0
..
.
(7.50)
(7.51)
TX(t)
= TX(t) TX(t) ,
(7.52)
where = 0 In ,
0 = jdiag( , 2, 1, 0, 1, 2, ).
(7.53)
125
HB for Floquet decomposition
Theorem 7.6.3 The eigenvalues and eigenvectors of the HB Jacobian matrix for (7.40)
(JHB = + TG(t) ) are given by the diagonal of 0 In I2M +1 and TU (t) , respectively.
The eigenvalues and eigenvectors of the HB Jacobian matrix for (7.46) (JHB = TGT (t) )
are given by the diagonal of 0 In + I2M +1 and TV (t) , respectively.
Proof Applying lemmas 7.6.1 and 7.6.2 to to (7.45), we have
TU (t) TU (t) + TU (t) T + TG(t) TU (t) = 0
( + TG(t) )TU (t) + TU (t) (T ) = 0,
(7.54)
(7.55)
126
7.6.4
Notations
Given any T -periodic matrix or vector A(t), sampled at N time points t0 , , tN 1 , we
denote block vector VA(t) and the block diagonal matrix DA(t) to be
A(t0 )
A(t0 )
A(t1 )
A(t1 )
VA(t) =
,
D
=
(7.56)
.
..
A(t)
..
.
.
A(tN 1 )
A(tN 1 )
Using these notations, two lemmas follow [133]:
Lemma 7.6.4 If Z(t) = X(t)Y (t), and X(t), Y (t) are T -periodic, then,
VZ(t) = DX(t) VY (t) ,
(7.57)
(7.58)
H = diag(
1
1
, ,
).
h0
hN 1
For example, 0,T D for backward Euler and central difference schemes are
1
1
0
1
1
1 1
1 0
1
0,BE =
, 0,CE =
.
.
.
.
.
.
.
.
.
.
. .. ..
1 1
1 0
1
(7.59)
(7.60)
127
Proof Applying lemmas 7.6.4 and 7.6.5 to to (7.45), we have
(T D + D(G(t) )VU (t) + VU (t) = 0,
(7.61)
(7.62)
7.6.5
Forcing Bi-orthogonality
Theoretically, matrices U(t) and V (t) from Floquet decomposition satisfies bi-orthonormality
V (t)U(t) = I which is an important result that is used in our derivation of the phase model.
However, bi-orthonormality may be destroyed in numerical methods, especially in our
methods where U(t) and V (t) are solved separately in two boundary-value problems of two
LPTV systems. To preserve bi-orthonormality, we need to scale one of them to restore
bi-orthonormality or nearly bi-orthonormality.
The workaround we employed in our implementation is: given U(t) = [u1 (t), , un (t)]
(t)
vn (t)
and V (t) = [v1 (t), , vn (t)], we scale V (t) to obtain V (t) = [ vTv1(t)~
, , vT (t)u
]. This
u1
n (t)
n
1
T
simple procedure forces vi (t)ui (t) = 1. One may further scale V (t) to make it bi-orthogonal
to U(t), but in our experiments, we observe good bi-orthogonality between U(t) and V (t).
T
128
7.7
(7.63)
(7.64)
subject to
dz
= v(t + z)(bp (t) bS (t + z))
dt
v(t)
(7.65)
7.8
129
7.8.1
(7.66)
The (nominal) periodic input bS (t) is set to bS (t) = cos(2t), and the periodic steady state
xS (t) is solved and shown in Fig. 7.3.
periodic steady state
0.1
state space
x1
0.06
x2
0.04
0.05
0.02
0
0.02
0.05
0.04
0.1
0
0.2
0.4
0.6
time
0.8
0.06
0.1
0.05
0
x
0.05
0.1
0.15
130
eigenvalues of HB matrices
300
100
original system
adjoint system
200
original system
adjoint system
50
100
0
100
50
200
300
30
20
10
10
20
30
100
50
50
100
150
200
250
131
x 10
U1, 1(t)
U1, 2(t)
0.005
7.8409
(t)
1, 1
100
63.7682
0.01
7.8409
200
63.7682
0.015
7.8409
0
3
x 10
0.5
time
U2, 1(t)
7.8409
0.015
7.8409
0.01
0.5
time
U2, 2(t)
0.5
time
300
1 63.76820
0.5
time
V (t)
T
1
63.7682
0
0.5
time
0.9
0
1
0
0.5
time
0.6
0.8
(b) V (t).
12
100
0.5
time
x 10
200
63.7682
0.5
time
V (t)
300
(a) U (t).
v (t) u (t)
2, 2
63.7682
1.1
2, 1
0.005
7.8409
(t)
1, 2
v1 (t) u2 (t)
4
2
0.5
time
T
v2
12
x 10
(t) u1 (t)
0
1
0
0.5
time
0.5
time
T
v
2
(t) u (t)
2
1.1
0.9
0
0.5
time
6
0
0.2
0.4
time
(d)
cT (t)V T (t)B
cT (t)c(t)
in (7.16).
132
variables of the phase model are confined on the orbit of xS (t), but they try to follow the
correct trajectory by choosing almost the closest point on the orbit of xS (t). The simulation
results for another initial condition z = 0.2 is also plotted in the x t space in Fig. 7.6(d),
and this is equivalent to applying an impulse function at time t = 0. Similar behaviors of
the phase model are observed, and it finally converges to the correct asymptotic response.
Indeed, if a constant phase-shifted input bS (t + p) is applied, no matter what the initial
condition is, our phase model will reproduce the correct phase-shifted asymptotic output.
This can be seen by the fact that z(t) = p is a fixed point of (7.16).
transient simulation of the phase model
0.1
x1 (full model)
x2 (full model)
x1 (phase model)
x2 (phase model)
0.1
0.05
0.2
0.3
0.4
0.05
0.5
0.6
0.1
0.7
0
(a) z(t).
(b) x(t).
0.06
0.04
0.04
0.02
0.02
x
0
0.02
0
0.02
0.04
0.04
0.06
0.1
0.1
0
0.1
x
time
0.1
time
Figure 7.6: Transient simulation of the phase model when b(t) = cos(2(t + 0.4)). In
Fig. 7.6(c) and Fig. 7.6(d), red(circled): phase model, Blue(solid): full model.
We then make the phase perturbation time-varying we apply an PM signal bp (t) =
cos(2(t + 0.1 sin(0.2t))), and the simulation results for 100 cycles are shown in Fig. 7.7. It
is seen in Fig. 7.7(b) that the response of the full model almost lies on the periodic orbit of
xS (t), and therefore, the phase model works perfectly.
We then apply a frequency perturbed input bp (t) = bS (1.2t), and the simulation results
133
transient simulation of the phase model
0.1
0.1
0.05
0.05
0
0
0.05
0.1
0.2
0.05
100
0
0.1
80
85
90
(a) z(t).
95
100
50
0.2
time
Figure 7.7: Transient simulation when bp (t) = cos(2(t + 0.1 sin(0.2t))). In Fig. 7.7(b),
red(circled): phase model, Blue(solid): full model.
are shown in Fig. 7.8. It is seen that the periodic orbit has a large deviation from that
of xS (t), and the phase model is doing its best to approximate the right trajectory using
points on xS (t). Most importantly, although the resulting time-domain waveforms do not
match exactly, the timing information is captured the frequency of the output waveform
is 1.2 which is the same as that of bp (t). Also note that the frequency perturbation can be
interpreted as a phase perturbation p(t) = 0.2t which can grow unboundedly as time evolves.
This shows that the perturbation can be very large as long as the underlying assumption
(that the trajectory does not change much) is satisfied.
Then an amplitude perturbed signal bp (t) = 1.2bS (t) is applied, and the simulation
results are shown in Fig. 7.9. Similar to previous results, the periodic orbit deviates from
that of xS (t), and the phase model produces a reasonably well-approximated waveforms.
Note that in many applications such as digital circuits and neuron models, voltage/potential
waveforms reach a saturated value in the nominal periodic solution xS (t). This fact makes
the periodic orbit insensitive to the amplitude perturbation, and therefore the assumption
that the trajectory stays close to xS (t) is satisfied. Therefore, the phase model generates
good results in these cases, as we will show in next two examples.
We have not provided results of LTI and LPTV reduced models since they simply produce
meaningless results. The assumption that the additive perturbation input must be small
is generally not satisfied. Specifically, in the case of the phase and frequency perturbations,
suppose the magnitude of bS (T ) is A, then the corresponding additive perturbation b(t) =
bS (t + p(t)) bS (t) can have a magnitude of 2A, and this usually breaks the small signal
assumption.
The speedups are not easy to measure considering various factors including memory
allocation and MATLAB optimization. For the examples we show in this paper, the measured
134
0.05
0.8
0.6
0.4
0.2
0.05
0.1
0.2
0
0.1
(a) z(t).
time
Figure 7.8: Transient simulation of the phase model when bp (t) = cos(2(1.2t)).
Fig. 7.8(b), red(circled): phase model, Blue(solid): full model.
In
0.06
0.05
0.04
0.04
0.02
0
x2
0.03
0.02
0.02
0.04
0.01
0.06
0.1
0
0.01
0.02
0
x1
1
(a) z(t).
0.1
0
2
time
Figure 7.9: Transient simulation of the phase model when bp (t) = 1.2 cos(2t). In Fig. 7.9(b),
red(circled): phase model, Blue(solid): full model.
135
transient runtime speedup is generally about 10 to 20. However, since the size of the
phase model is 1, the speedups are expected to be much larger for larger systems, similar to
PPV models [29].
7.8.2
The firing neuron model we consider is known as Morris-Lecar model [134]. The differential equations for a single neuron are
1
d
(gL (V VL ) gCa M (V VCa ) gK N(V VK ) + I)
V =
dt
CM
d
1
N =(N N) ,
dt
N
(7.67)
60
55
50
45
10
40
15
35
30
0
20
0
50
100
50
T
100
150
time
150
(b)
c (t)V (t)B
cT (t)c(t)
in (7.16).
136
doubled, the periodic orbit only deviates a little bit, and the phase model results almost
match those of the full model.
transient simulation of the phase model
9
30
V (full model)
N (full model)
V (PSS)
N (PSS)
V (phase model)
N (phase model)
20
7
10
10
4
3
20
30
40
0
0
100
200
300
400
500
(a) z(t).
100
200
300
400
500
7.8.3
We apply the phase macromodeling technique to a nine-stage inverter chain circuit. The
nominal periodic input is set to a square wave shown in Fig. 7.12(a), and the time-varying
T
T (t)B
function c c(t)V
is plotted in Fig. 7.12(b). This circuit is a typical digital circuit in the
T (t)c(t)
sense that each node voltage settles to a saturated value in a short time, and this makes
the trajectory in the state space less sensitive to input perturbations. However, the timing
properties are sensitive to input perturbations.
We have tried several kinds of perturbations and similar results to those of previous
examples are observed. Here, we show a special perturbation that is of interest in timing
analysis the perturbation is on the slope of the input ramp. In timing analysis and
library generation, people normally compute a look-up table storing a map from input slope
to output delay and slope. We show that our phase model reproduces accurate timing
properties.
In Fig. 7.13, we plot the transient simulation results of the phase model and the full model
when the slope of the rising edge is 0.6 (instead of 0.5 for the nominal input). It is observed
that delays of waveforms change due to perturbation, and the phase model matches that of
the full model.
137
5000
4000
3000
0.5
2000
1000
0
1000
0.5
2000
3000
0
1
0
10
10
15
time
15
(b)
cT (t)V T (t)B
cT (t)c(t)
in (7.16).
x1 (full model)
x2 (full model)
x1 (unperturbed PSS)
x2 (unperturbed PSS)
x1 (phase model)
x2 (phase model)
0.8
0.1
0.6
0.05
0.4
0.2
0
0.05
0.2
0.1
0.4
0.15
0.6
0.2
0.25
0
0.8
5
10
(a) z(t).
15
10
(b) ~x(t).
15
138
Chapter 8
DAE2FSM: Finite State Machine
Abstraction from Differential
Equations
In this chapter, we study a different problem in which a finite state machine, instead
of a smaller set of differential equations, is abstracted from original differential equations.
The problem is motivated by the fact that in modern digital designs, circuit behaviors may
not be accurately modeled by the ideal pre-specified finite state machines, due to parameter
variations, crosstalk and coupling, environment changes, etc.. As a result, a ground-up
abstraction of the finite state machine from differential equations, deemed as the golden
model, is needed. Other than the circuit applications, the model is also found useful in other
engineering problems, such as biological systems.
We start with a detailed comparison of DAE/ODE models and FSM models that shed
light on the validity of such abstractions and the approach we may follow. Then we review several existing approaches in the community of machine learning1 , and discuss how
we adapt these algorithms in learning differential equations. The key challenge is to develop
corresponding concepts/subroutines in those of machine learning, such as state merging and
equivalence checking. We present experimental results of several systems, and show applications of the finite state machines learned in further system-level simulation and analysis.
8.1
Problem Formulation
The problem of FSM abstraction is one variation of the problem (2.1) in Section 2.1,
where M is the set of ODEs/DAEs, and Mr is the set of finite state machines. The fundamental questions are similar, i.e.,
1
The term machine learning here refers to the actual finite state machine learning, rather than statistical
learning such as learning Bayesian networks.
139
1. What are F (, ), G(, ) and H(, ) that lead to a reasonable FSM model?
2. Can we, or how do we efficiently solve the optimization problem?
To answer these questions, we compare DAE models and FSM models to identify challenges and clarify ambiguities.
8.1.1
DAEs and FSMs are defined in a similar way they both model the dynamics of states
under input excitations and produce outputs as a function of time. The major difference is
that the DAE model is a continuous model while the FSM model is a discrete model.
Take Mealy FSMs as an example. A Mealy machine is defined by a 6-tuple (Q, Q0 , , , , G),
where Q = {q1 , , q|Q| } is a finite set of states, Q0 Q is the initial state, = {1 , , || }
is a finite set of input alphabet, = {1 , , || } is a finite set of output alphabet,
: Q Q is the transition function that maps a pair of a state and an input symbol to
the next state, and G : Q is the output function that maps a pair of a state and an
input symbol to an output symbol. The evolution of the state in an FSM is determined by
the transition function, and is a discrete-time discrete-value trace si , i = 0, 1, 2, , where
si Q. Similarly, the input-output pair of an FSM is a pair of discrete-time discrete-value
traces (ui , yi ), i = 0, 1, 2, , where ui and yi .
In contrast, in differential equations, the evolution of state variables x RN is determined
by the underlying differential equation, or equivalently, a state transition equation
x(t) = (x0 , u(t)),
(8.1)
which maps a pair of current state x0 and input signals u(t) to states in the future x(t). The
evolution x(t), t R is therefore a continuous-time continuous-value trajectory. Similarly,
the input-output pair of the DAEs is a pair of continuous-time continuous-value trajectories
(u(t), y(t)), t R.
Comparison of different components in DAE models and FSM models are summarized
in Table 8.1.
Table 8.1: Comparison of differential equations
Model
Differential Equations
Time
continuous, t R
State Space
x RN (infinite)
Input
u R (infinite)
Output
y R (infinite)
Update Equation x(t) = (x0 , u), : RN R RN
Output Equation
y = h(x, u), h : RN R R
140
8.1.2
The difficulty of formulating the FSM abstraction problem lies in the inherent differences
between DAE models and FSM models: DAE models are continuous and the state space is
infinite, while FSM models are discrete and the state space is finite. Then, does it make
sense for an FSM model to be equivalent to a DAE model? If not, what does it mean for an
FSM model to approximate a DAE model?
For the first question, the answer is at least partially yes. In digital circuit design,
we essentially build up circuits that realize a specific FSM (e.g., a latch). With the valid
consensus that circuits are well modeled by differential equations, it makes sense to say that
these DAEs are equivalent to the corresponding FSMs.
Example 8.1.1 (An RC circuit) Consider an RC circuit whose state variable x is governed by the differential equation
d
x = x + u,
(8.2)
dt
where u is the binary input being 0 or 1. The output y = x is sampled at time t = n where
n Z, and is deemed as 0 if x(n) > 0.5, and 1 otherwise.
Let the initial condition be x0 . If input is 0, the output at the next sample time is
x(1) = x0 e1 .
(8.3)
(8.4)
From (8.3) and (8.4), we see that x(t) cannot overshoot, and therefore 0 x(t) 1.
In order for the circuit to emulate a latch circuit, i.e., the next output always equals to
the current input, we compute the worst-case scenarios. Assuming the current input is 0,
then the next output is
x(1) = x0 e1 e1 0.38 < 0.5.
(8.5)
Therefore, the next output is always 0 if the current input is 0. Similarly, we can show that
the next output is always 1 if the current input is 1.
In summary, (8.2) realizes the FSM shown in Fig. 8.1, which is the FSM of an ideal latch.
From Example 8.1.1, we have the following notion of equivalence between a DAE model
and an FSM model:
Definition 8.1.1 (Equivalence between a DAE model and an FSM model) Suppose
we have a mapping : {ui } 7 u(t) from discrete-time discrete-value input sequences (of the
FSM model) to continuous-time continuous-value input waveforms (of the DAE model), and
141
0/0
0
1/1
1/0
0/1
142
Definition 8.1.2 (U-Equivalence between a DAE model and an FSM model) Suppose
we have a mapping : {ui } 7 u(t) from discrete-time discrete-value input sequences (of the
FSM model) to continuous-time continuous-value input waveforms (of the DAE model), and
a mapping : y(t) 7 {
yi } from continuous-time continuous-value output waveforms (of the
DAE model) to discrete-time discrete-value output sequences (of the FSM model.
If for all input sequences {ui } U which is a subset of input sequences, the discretized
output of the DAE model matches the output of the FSM model, yi = yi , i = 0, 1, 2, , then
the DAE model and the FSM model are called to be U-equivalent with each other.
For example, U can be the set of input sequences of fixed length L. This is indeed
reasonable in many practical applications because system state evolves from a steady state
for some time and gets reset once in a while. For example, in a communication system or a
high-speed link, the circuits are reset after each packet is transmitted/received. Therefore,
the models will be good enough if they matches the responses to input sequences of length
equal to the packet size.
This relaxed definition therefore specifies an FSM approximation to the original DAE
model. By specifying different Us in different applications, we obtain different levels of
approximations.
Another notion of equivalence may be defined by discretizing the continuous state space
of the DAE model and adding transitions among different regions. Denote X(si ) to be the
region in RN corresponding to the state si . The transition from one state si to another state sj
is added if there exists a point xi X(si ) and an input ui such that (xi , (ui )) X(sj ).
For example, in Fig. 8.2, the state 0 corresponds to the half space x < 0.5 and the state 1
corresponds to the half space x 0.5. and the transitions are added by the rule above.
While this conversion from a DAE model to an FSM model seems to be straightforward
for the above simple example, it is not generally the case. Consider a similar example as
Example 8.1.1, where the dynamics of x is defined by
d
x = 0.5x + u.
dt
(8.6)
In this case, if we still use the same partitioning, then we may have an FSM with more than
4 transitions, leading to a non-deterministic FSM. Also, in general, it may not be possible to
examine every point in a region, and some transitions may be lost. While a non-deterministic
FSM may be useful, it may be useless if there are too many transitions that little information
of the DAE behaviors is modeled.
143
8.1.3
Based on the above definition, we can formulate a simplistic FSM abstraction problem
in the form of (2.1):
minimize
FSMs
(8.7)
where we set F (, ) to be constant 1, H(, ) to be empty, and G(, ) is the difference between
the discretized outputs of the DAE model and the outputs of the FSM model.
In another word, in problem (8.7), we are given a set of input-output trajectories of a
DAE model, and we try to find an FSM model that matches these input-output trajectories.
It appears that there have been efficient algorithms developed in other areas such as software verification and machine learning that solves similar problems. We will study some of
them in the following sections, and develop sufficient adaptations and modifications of these
algorithms to solve our problem.
The optimization problem (8.7), however, does not have a unique solution. To ensure a
unique solution, we may formulate more sophisticated problems by adding/choosing other
F, G, H functions. For example, we may minimize the number of states of the FSM model,
or enforcing deterministicity of the FSM model, by modifying F and H respectively.
8.2
8.2.1
144
observation table is further split into the top part and the bottom part, where row labels
of the bottom part are obtained by concatenating row labels of the top part with all input
alphabets. We will regard each row in the top part of the table as one state in the state
machine, and therefore for each state in the top table, there must exist a row in the table
that corresponds to its next state.
With the observation table, the key idea is then to differentiate different states by experiments, i.e., we regard two rows as the same state if entries of the two rows are the same.
For example, row 1 has the same entries as row 110 in Table 8.2, and therefore represent
the same state. However, to conjecture a DFA from an observation table, we need to place
two constraints on the table:
1. The table is closed, i.e., every row in the bottom part of the table has a corresponding
row in the top part with identical results in all columns. This guarantees that none of
the rows in the bottom part corresponds to a new state.
2. The table is consistent, i.e., every pair of rows in the top part of the table with identical
results in all columns also has identical results when any alphabet symbol is added.
Mathematically, s1 , s2 S, if row(s1 ) = row(s2 ), then a , row(s1 a) = row(s2 a).
This guarantees that any pair of states that we think are the same are still the same
if the experiment string is extended by one input alphabet symbol.
The procedures to make the observation table closed and consistent are simple, and can be
found in [135].
Given a closed and consistent observation table (S, E, T ), we can then conjecture a DFA
(Q, , , Q0 , F ) defined by
Q = {row(s) : s S of top table},
Q0 = row(),
(row(s), a) = row(s a), a ,
F = {row(s) : T (s, ) = 1}.
(8.8)
Using this DFA, we can feed it to an equivalence checker to check if the DFA is equivalent to the original one. If equivalence checker returns yes, we terminate the algorithm.
Otherwise, it returns an counter-example. In that case, we add the counter-example into the
observation table, and repeat the procedure of conjecturing a new DFA until convergence.
Summarizing the above procedure, the Angluins algorithm is sketched in Algorithm 13.
8.2.2
As defined in Section 8.1.1, the Mealy machine is similar to a DFA. It does not have
accept states in DFA, but has an output function. This difference requires us to re-define
145
Algorithm 13 Angluins DFA learning algorithm
1: Construct the initial observation table (S, E, T ) with S = {} and E = {}.
2: repeat
3:
Make the observation table (S, E, T ) closed and consistent.
4:
Conjecture a DFA according to (8.8).
5:
Check the equivalence between the conjectured DFA and the original DFA.
6:
if not equivalent then
7:
Add the counter-example in the observation table.
8:
end if
9: until Equivalence checker returns yes.
the meaning of the observation table.
In Angluins algorithm, each cell in the table is filled in by yes/no answers to whether the
DFA accepts an input sequence. This can be viewed as the output of the DFA. Similarly,
for Mealy machines, we can fill in each cell in the table by the output of the Mealy machine.
It turns out that using this new definition of the observation table, the concepts of closedness and consistence of the observation table still holds [136], and we can conjecture a Mealy
machine (Q, Q0 , , , , G) from this new observation table by
Q = {row(s) : s S of top table},
Q0 = row(),
(row(s), a) = row(s a), a ,
G(row(s), a) = T (s, a), a .
(8.9)
Therefore, in order to learn a Mealy machine from another Mealy machine using Angluins
idea, we just need to use the new observation table, and replace step 4 in Algorithm 13 by
Conjecture a Mealy machine according to (8.9).
8.3
146
obtain the output. Therefore, we need a simulator that can produce I/O traces of a
machine/system.
2. The equivalence checker. The algorithm decides whether to terminate by making
equivalence queries (i.e., is the new machine M equivalent to the original machine?).
8.3.1
8.3.2
To implement the equivalence checker, we use the U-equivalence defined in Section 8.1.2.
Therefore, we call the SPICE simulator to perform a set of simulations of circuit differential
equations using u(t) derived from {ui } U. If the sampled outputs {y(ih)} lie in the corresponding quantization regions defined for {yi } (the output of the FSM to input sequence
{ui }), then we return yes, i.e., the FSM is equivalent to the differential equations. Otherwise, the input sequence is reported as an counter-example, and is added into the observation
table to refine the FSM.
147
The training set used in the equivalence checking is important. One common implementation is to use a Monte-Carlo sampling of input signals if a priori statistics of inputs
are known. Another implementation useful in digital applications is to try out all digital sequences of a fixed length N. Although this leads to 2N possible sequences, it is expected that
for practical circuits with small number of I/O ports, short sequences are enough to characterize behaviors of the circuit. Furthermore, we can add/remove some special sequences
into the training set to ensure some desirable properties of the system. For example, we may
say that we want to capture results by applying five 1s at the input, and therefore we add
the input sequence of five 1s into the training set. For another example, we may have a
constraint on the input signal that it is not possible to have inputs of six or more consecutive
1s (due to properties of the circuit generating the input signal), and then the training set
can be pruned to speedup the equivalence checking procedure.
8.3.3
The Mealy machine extracted from differential equations has a finite set of states. Since
the FSM is sometimes extracted in order to understand system behaviors, it poses a question
of what these states actually mean in the original system.
According to Angluins algorithm, a state is defined to be the state that the system
reaches by applying a corresponding input sequence. Therefore, strictly speaking, each state
in the FSM corresponds to a single point in the continuous state space.
On the other hand, similar to U-equivalence of two FSMs, we may define the U-equivalence
of two states. Namely, suppose x1 , x2 are two states in the full continuous state space, if for
all u(t) U, the sampled and quantized versions of y1 = h(x1 , u) and y2 = h(x2 , u) match,
then x1 and x2 are U-equivalent. Using this definition, we see that there is a region around
each point so that all point in that region are U-equivalent.
For U being all binary inputs of length k, this definition partitions the continuous state
k
space into 2k2 regions. This is a huge number. However, in real circuits, many of these
regions are either never/rarely visited or extremely small, giving the premise that our algorithm can work well.
Solving for the boundaries of such regions, however, may be complicated and expensive.
It is therefore preferable to obtain an approximation. To do this, we use the discrete points
obtained above to construct a Voronoi tessellation/decomposition/diagram [137] which partitions a space with n points into convex polygons such that each polygon contains exactly
one generating point and every point in a given polygon is closer to the generating point
than to any other. These polygons, known as Voronoi polygons, are then associated with
states in the FSM. As an example, a 3D Voronoi diagram with 4 points is shown in Fig. 8.7
in Section 6.5.
Due to properties of the Voronoi tessellation, it can be viewed as a natural partitioning
148
of the continuous state space given that the n points are obtained by the definition of each
conjectured state in the algorithm. With this partitioning, it is then possible to apply
other techniques to further exploit low-level dynamics of the circuit. For example, we can
construct simple nonlinear/linear systems for each partition, and this is very similar to ideas
in hybrid system abstraction of continuous dynamical systems [138], and trajectory-based
MOR techniques [2].
The key difference between this approach and many previous methods (described in
Section 8.2) on partitioning the continuous state space is that previous methods often suffer
from curse of dimensionality due to the discretization of all state variables, while our approach
avoids this problem by discretizing the state space directly. As long as the original system
can be represented by a small FSM, our approach will be much more efficient. The number
of partitions we obtain is the same as the size of the FSM, and we just need to store one
point for each state to store the partitioning.
8.4
8.4.1
We consider a simplified latch circuit shown in Fig. 8.3. The latch is composed of a
sampler (consisting of a multiplexer) and a regenerator (consisting of two cross-coupled
inverters). Each output of the multiplexer and inverters is followed by an RC circuit that is
not shown in Fig. 8.3. The latch samples the DAT A input when CLK is 1. When CLK
is 0, two inverters are cross-coupled, and re-generate the output (QBB).
Q
QB
INV
CLK
MUX
DATA
QBB
INV
149
The differential equations of this latch circuit are
RC
dx1
+ (x1 fM U X (u1 , u2 , x3 )) =0,
dt
dx2
RC
+ (x2 fIN V (x1 )) =0,
dt
dx3
+ (x3 fIN V (x2 )) =0,
RC
dt
(8.10)
where state variables x1 , x2 , x3 are node voltages vQ , vQB , vQBB , respectively; inputs u1 , u2
are CLK and DAT A, respectively; fIN V and fM U X are I/O functions defined by
fIN V (x) = tanh(10x),
A+B BA
+
tanh (10C) ,
fM U X (C, A, B) =
2
2
(8.11)
F
C
INV
MUX
B A
Note that Mealy machines shown in Fig. 8.6 are indeed Moore machines, simply due to the fact that the
150
Figure 8.5: Response of the latch when the input sequence is 10011000001001100000 and
Vsw = 1.6V .
are not arbitrary. They are values of rows in the observation table, i.e., the quantized output
with respect to corresponding input sequences. As an example, the observation table of the
latch with input swing 1.6 is shown in Table 8.2. We see that the state 0,0,0,0 corresponds
to row , the state 0,0,1,0 corresponds to row 1, etc.. From Table 8.2, we can derive the
FSM in Fig. 8.6(b) using (8.9).
Based on the Mealy machines in Fig. 8.6, we can claim that when input swing is 2V , the
latch behaves like an ideal latch; when input swing is 1.6V , the input signal must persist
for two consecutive cycles to avoid errors; when input swing is 1.3V , the input signal must
persist for three consecutive cycles to avoid errors. These properties can either be directly
observed in simple FSMs, or be checked by model checking tools. To validate these claims, we
perform many further SPICE simulations, and the same qualitative behaviors are observed.
Using these Mealy machines, we then construct Voronoi diagrams using discrete points
corresponding to each state. For example, the Voronoi diagram of the Mealy machine in
Fig. 8.6(b) is shown in Fig. 8.7. The four points correspond to four states in Fig. 8.6(b), and
the (blue) trajectories in Fig. 8.7 are simulation traces that leads to the discrete points of
different states, and show typical transition dynamics of the latch. This partitioning of the
state space also provides direct intuition on how circuit behave: the two states at two corners
represent solid 0 and solid 1 states; the two states in the middle capture dynamics in
the metastable region, and can be viewed as weak 0 and weak 1 states.
output y is just the output voltage which corresponds to a state variable x3 . This is very common in circuit
simulations, and a simple post-processing can be performed to obtain a Moore machine.
151
0/0
1/0
0/1
0,0
1/1
0/0
1,1
0,0,0,0
1/1
1/0
0/0
0,0,1,0
1/0
1,1,1,0
0/1
(a) Vsw = 2V .
1/1
0/1
1,1,1,1
0/0
1,1,1,1,1,1
1/0
0,0,0,0,0,0
0,0,0,1,0,0
1/1
1/0
1/1
0/0
0,0,1,1,0,0
0/0
1/0
0/1
1,1,1,1,0,0
0/1
0/1
0 0 0 0
1
0 0 1 0
11
1 1 1 0
110
0 0 1 0
111
1 1 1 1
1110 1 1 1 0
11100 0 0 1 0
0
0 0 0 0
10
0 0 0 0
1100 0 0 0 0
1101 1 1 1 0
1111 1 1 1 1
11101 1 1 1 1
111000 0 0 0 0
111001 1 1 1 0
1,1,1,1,1,0
152
153
8.4.2
An Integrator Circuit
The second example is an integrator circuit shown in Fig. 8.8, where R = 1, C = 1, and
the opamp is approximated by a controlled voltage source which clips at 1V and has gain
10 and resistance 0.1. The input is a binary sequence where 1 and 0 are at voltage level
+Vin and Vin , respectively, and each bit lasts for 1 second.
Vin
Vout
154
components of the DAG (Directed Acyclic Graph) of the FSM, and this typically filters out
a lot of transient states. Therefore, to discard useless transient states, we apply algorithms
for finding strongly connected components, and fortunately, linear-time (with respect to the
number of states) are available, such as [139].
0/1
1/0
1/1
0/1
1,1,1,1
1/1
1,1,1,0
1/1
0/0
1,1,0,0
0/1
0,0,0,0
1/1
1,1,1,0,0,0,0,0,1,1,0,1,1,1,1
1/1
0/1
1,1,1,1,0,0,1,0,1,1,0,1,1,1,1
1,1,0,0,0,0,0,0,1,1,0,0,1,1,0
1/1
0/1
0/1
1/1
1,1,1,1,1,1,1,1,1,1,0,1,1,1,1
0/0
1,1,1,0,0,0,1,0,1,1,0,0,1,1,1
1/1
0,0,0,0,0,0,0,0,1,1,0,0,1,1,0
1/1
0/1
1/1
1/1
1,1,1,1,1,0,1,1,1,1,1,1,1,1,1
0/1
1,1,0,0,0,0,1,0,1,1,0,1,1,1,1
0/1
1,1,1,0,0,0,1,1,1,1,0,1,1,1,1
1/0
0/1
1,1,1,1,0,0,1,1,1,1,1,1,1,1,1
0/0
0,0,0,0,0,0,0,0,0,1,0,0,1,0,0
0/0
0/1
1,1,1,1,0,0,1,1,1,1,0,1,1,1,1
0/1
1/1
0/1
1,1,1,0,0,0,1,0,1,1,0,1,1,1,1
1/0
1/1
0/1
1/0
1,1,0,0,0,0,0,0,1,1,0,0,1,1,1
0/0
1/1
1/1
0,0,0,0,0,0,0,0,0,1,0,0,1,1,0
1/0
1/0
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
0/0
155
Chapter 9
Conclusion and Future Work
Model order reduction for nonlinear dynamical systems has been a difficult problem for
many years. This dissertation provides a self-contained view of the MOR problem, and
presents several algorithms, ManiMOR, QLMOR, NTIM, DAE2FSM, to solve the problem
from different perspectives.
Starting with the linear projection framework used by most previous techniques, we study
and analyze its drawbacks in applications to reduce nonlinear systems. With that, we present
a more general nonlinear projection framework which we believe is the natural and correct
framework for nonlinear model order reduction. Indeed, the techniques presented in previous
chapters may all be interpreted in this framework.
In its most straightforward form, we have presented ManiMOR which explicitly constructs and parameterizes a nonlinear manifold, and projects the full model on the manifold
to obtain a reduced model. The projected model, by construction, matches the DC and
AC responses of the full model, and can be proven to be much more efficient than linear
projection-based methods.
By exploiting internal structures of the full model, we have proposed QLMOR which
creates a QLDAE reduced model. The key contribution of this work is the development of a
symbolic transformation procedure that creates an equivalent quadratic-linear model to the
full model. By doing so, it has avoided all problems caused by Taylor series truncation, such
as local convergence, high computational cost, passivity/stability loss, etc.. The other key
contribution is the development of a reduction algorithm of QLDAEs which guarantees to
match certain moments of the Volterra kernel for QLDAEs. Therefore, this reduced model
is especially of great use for capturing harmonic distortions in nonlinear systems.
Using a very specific nonlinear manifold/projection, we have presented NTIM, a phase
macromodel that is proven to capture phase responses in nonlinear systems. While we develop the model using Floquet theory, the phase macromodel naturally fits into the nonlinear
projection framework. Due to the importance of timing/phase responses in practical systems,
this model is shown to be useful in many scenarios.
Besides reduced order models in the form of differential equations, we also try to see if
156
other types of reduced models are appropriate. In particular, we presented DAE2FSM that
abstracts a finite state machine from a set of differential equations. This model is especially
useful if the full model behaves as a finite state machine in that case, the algorithm
guarantees to find the corresponding finite state machine. FSM models can be further fed into
model checking tools to verify higher-level system properties, and have various applications
such as mixed-signal system simulation and biological system abstraction.
Rather than completely solving the nonlinear MOR problem, our exploration leads to
many more research directions to pursue. In the following, we comment on further improvements that may be done for the methods we present in this dissertation. We then sketch a
few potential theoretical and practical improvements/extensions/considerations.
9.1
9.1.1
1. The Krylov distribution defined by (4.22) is generally not involutive. However, there
are systems of which (4.22) is involutive. It is interesting to see what systems have such
properties and what such properties correspond to in real circuit implementations.
2. By construction, ManiMOR ensures that the DC attractors are included in the manifold. This satisfies the necessary condition that the reduced order model can reproduce
attractors of the full model. However, this is not sufficient for the stability to be preserved. It is useful to show under what conditions can the stability be preserved.
3. In the derivation of ManiMOR, we have assumed that f () in the differential equations
is invertible. However, this does not always hold. For example, in bi-stable systems
such as SRAMs, one DC input corresponds to two stable DC solutions. One possible
solution is to find a manifold that covers all DC solutions. Another possible solution is
to modify the circuit to a decoupled (open-loop) one such that f () becomes invertible,
and then reduce the open-loop model.
4. MIMO systems may be handled by generalizing the manifold definition. For example,
the DC manifold, which is a one-dimensional manifold for SIMO models, may be
defined as a manifold of larger size. For example, Fig. 9.1 shows a 2-input NAND
circuit and the 2-dimensional DC manifold. Following this idea, we may also extend
the definition of AC and ManiMOR manifolds.
5. For oscillators, DC equilibria are not the attractors. Instead, the periodic solution,
or the limit cycle in the state space is the attractor. Therefore, the definition of the
manifold should be changed for oscillators. It seems that this has several connections
to the NTIM reduced model that deserve further investigations.
157
Vdd
in1
in2
out
1.5
in1
1
0.5
0
0
0.5
in2
0.5
1.5
in
1
1.5
in
(b) DC surface
9.1.2
QLMOR
1. Polynomial systems are easier to handle not only in model order reduction, but also
in simulation and verification. For example, SOStools [140] may be used to construct
Lyapunov functions for polynomial systems; the simulation of quadratic systems may
be much faster than more complicated nonlinear systems. Therefore, the quadraticlinearization procedure serves as a key transformation to enable the application of tools
that handle polynomial systems to non-polynomial systems.
2. Chemical rate equations contain only polynomial functions in the ODEs. Therefore,
there is a one-to-one correspondence between a set of polynomial ODEs and a set of
chemical reactions. Therefore, we may synthesize chemical reactions from ODEs even
if they are non-polynomial.
158
3. From a different perspective, the quadratic-linearization procedure shows that any
nonlinear system may be written as a QLDAE. Therefore, the QLDAE form is an
acceptable form for nonlinear system identification, and it represents a nonlinear system
by matrices. We may modify current system identification methods to identify QLDAE
models.
9.1.3
NTIM
1. NTIM in Chapter 7 generates phase models for only single-input systems. It is not
applicable to multiple-input systems because these system will settle in a quasi-periodic
steady state. From a projection view point, the manifold is no longer the limit cycle,
but a torus in the state space. It appears that multiple phase variables and multiple
phase equations need to be defined for systems in quasi-periodic steady states.
2. The phase model may be used in timing analysis of a large digital circuit we may
replace each digital gate by its phase model, therefore reducing the number of equations
describing each gate to 1. This greatly speeds up the simulation, while preserves timing
behaviors.
9.1.4
DAE2FSM
159
9.2
1. While we have mainly considered the model reduction of ODE/DAE models, the model
order reduction of other models in Table 1.1 deserve further investigation. For example,
we may reduce an FSM model by merging/removing states. The question is again what
are the important states, and how to identify these states. A straightforward idea
is to preserve the attractors of the FSM. Similarly, we can explore the model order
reduction of other models.
2. Chemical reactions may also be viewed as a model, although it not listed in Table 1.1.
While it describes chemical kinetics, it is not strictly a mathematical model because it
may be interpreted as many different models, deterministic or stochastic, finite-state
or infinite-state. The reduction of the chemical reactions is possible by ignoring the
reactions that have a high reaction rate, we may obtain a model that captures the slow
kinetic dynamics. However, this method is not optimal generally in other senses
(such as attractor-preserving, moment-matching, etc.). It is interesting to explore the
application of other MOR techniques or develop new MOR techniques for the chemical
reactions.
3. We may also look at the problem of abstracting one model from another. For example,
from an ODE/DAE model to a hybrid model. Such abstractions can have a lot of
applications, e.g., to simplify the implementations/synthesis of a controller.
4. One important aspect of the projection framework is that it is computationally efficient
to obtain a projection that leads to a reasonable reduced order model. However,
recently there have been many progresses in improving the efficiency of optimization
algorithms. The direct application of the state-of-the-art optimization routines may
be possible.
5. We have mentioned in Section 2.1.3 the problem of formal verification of analog/mixedsignal systems. However, such algorithms and tools are extremely immature, and
deserves further investigation. While this is not about model order reduction, it closely
relates to MOR in that such tools may be used to check the correctness of reduced
models.
6. Practical systems may be redundant. For example, the LTI system matrices may be
rank-deficient. It looks like MOR method may remove the redundancy in such systems.
It is interesting to see if MOR methods can actually preserve certain redundancies.
7. In several structure-preserving MOR techniques, certain structures (of LTI systems)
are preserved so that the reduced model corresponds to an RLC circuit implementation.
A similar problem is how to obtain a nonlinear model that corresponds to a transistorlevel implementation.
160
8. There have been some literature discussing how to do nonlinear balanced truncation to
obtain a reduced model. However, the theory of nonlinear balancing has not reached
a consensus, and the existing methods are computationally extremely inefficient.
9. The theory for the balanced realization of LTV systems seems to be well-established.
However, the circuit community have not been aware of these techniques. We may try
these techniques and see how they work for circuit examples.
10. We have shown that some nonlinear systems, under appropriate coordinate transformations, may be converted to equivalent linear systems. However, the question is why
would we want to implement an LTI system using a nonlinear system? One possibility is that the LTI implementation may require unrealistic resistance and capacitance
values that cannot be mapped to an integrated circuit. In that case, it would be great
if we can map the LTI system to a nonlinear system that can be implemented on chip.
11. The sparsity is typically destroyed by MOR methods. However, the discussion in
Chapter 3 reveals that sometimes we have extra degree of freedom in choosing a projected model. These degree of freedom may be used to achieve the sparsity of the
reduced order model.
161
Bibliography
[1] R.T. Chang et al. Near Speed-of-Light Signaling Over On-chip Electrical Interconnects.
Solid-State Circuits, IEEE Journal of, 38:834838, May 2003.
[2] Michal Rewienski and Jacob White. A Trajectory Piecewise-Linear Approach to Model
Order Reduction and Fast Simulation of Nonlinear Circuits and Micromachined Devices. IEEE Transactions on Computer-Aided Design, 22(2), February 2003.
[3] Y. Wang and H.G. Dohlman. Pheromone Signaling Mechanisms in Yeast: A Prototypical Sex Machine. Science, 306:15081509, Nov 2004.
[4] K.C. Sou, A. Megretski, and L. Daniel. A quasi-convex optimization approach to
parameterized model order reduction. Computer-Aided Design of Integrated Circuits
and Systems, IEEE Transactions on, 27(3):456469, 2008.
[5] L. Daniel, O.C. Siong, L.S. Chay, K.H. Lee, and J. White. A multiparameter momentmatching model-reduction approach for generating geometrically parameterized interconnect performance models. Computer-Aided Design of Integrated Circuits and
Systems, IEEE Transactions on, 23(5):678693, 2004.
[6] B.N. Bond and L. Daniel. A piecewise-linear moment-matching approach to parameterized model-order reduction for highly nonlinear systems. Computer-Aided Design
of Integrated Circuits and Systems, IEEE Transactions on, 26(12):21162129, 2007.
[7] P. Li, F. Liu, X. Li, L.T. Pileggi, and S.R. Nassif. Modeling interconnect variability
using efficient parametric model order reduction. In Proceedings of the conference on
Design, Automation and Test in Europe-Volume 2, pages 958963. IEEE Computer
Society, 2005.
[8] X. Li, P. Li, and L.T. Pileggi. Parameterized interconnect order reduction with explicitand-implicit multi-parameter moment matching for inter/intra-die variations. In iccad,
pages 806812. IEEE, 2005.
[9] T. Bui-Thanh, K. Willcox, and O. Ghattas. Model reduction for large-scale systems
with high-dimensional parametric input space. SIAM Journal on Scientific Computing,
30(6):32703288, 2008.
162
[10] Y. Liu, L.T. Pileggi, and A.J. Strojwas. Model order-reduction of RC (L) interconnect
including variational analysis. In Proceedings of the 36th annual ACM/IEEE Design
Automation Conference, pages 201206. ACM, 1999.
[11] A. Odabasioglu, M. Celik, and L.T. Pileggi. PRIMA: Passive reduced-order interconnect macromodeling algorithm. In Proceedings of the 1997 IEEE/ACM international
conference on Computer-aided design, pages 5865. IEEE Computer Society, 1997.
[12] J.M. Wang, C.C. Chu, Q. Yu, and E.S. Kuh. On projection-based algorithms for
model-order reduction of interconnects. Circuits and Systems I: Fundamental Theory
and Applications, IEEE Transactions on, 49(11):15631585, 2002.
[13] I.M. Elfadel and D.D. Ling. A block rational Arnoldi algorithm for multipoint passive model-order reduction of multiport RLC networks. In Proceedings of the 1997
IEEE/ACM international conference on Computer-aided design, pages 6671. IEEE
Computer Society, 1997.
[14] Y. Shin and T. Sakurai. Power distribution analysis of VLSI interconnects using model
order reduction. Computer-Aided Design of Integrated Circuits and Systems, IEEE
Transactions on, 21(6):739745, 2002.
[15] A.C. Cangellaris, M. Celik, S. Pasha, and L. Zhao. Electromagnetic model order reduction for system-level modeling. Microwave Theory and Techniques, IEEE Transactions
on, 47(6):840850, 1999.
[16] H. Wu and A.C. Cangellaris. Model-order reduction of finite-element approximations of
passive electromagnetic devices including lumped electrical-circuit models. Microwave
Theory and Techniques, IEEE Transactions on, 52(9):23052313, 2004.
[17] T. Wittig, I. Munteanu, R. Schuhmann, and T. Weiland. Two-step Lanczos algorithm
for model order reduction. Magnetics, IEEE Transactions on, 38(2):673676, 2002.
[18] R.D. Slone, R. Lee, and J.F. Lee. Multipoint Galerkin asymptotic waveform evaluation for model order reduction of frequency domain FEM electromagnetic radiation
problems. Antennas and Propagation, IEEE Transactions on, 49(10):15041513, 2001.
[19] M.E. Kowalski and J.M. Jin. Model-order reduction of nonlinear models of electromagnetic phased-array hyperthermia. Biomedical Engineering, IEEE Transactions on,
50(11):12431254, 2003.
[20] H. Yu, Y. Shi, and L. He. Fast analysis of structured power grid by triangularization
based structure preserving model order reduction. In Proceedings of the 43rd annual
Design Automation Conference, pages 205210. ACM, 2006.
163
[21] H. Su, E. Acar, and S.R. Nassif. Power grid reduction based on algebraic multigrid
principles. 2003.
[22] T.V. Nguyen and J.M. Wang. Extended Krylov subspace method for reduced order
analysis of linear circuits with multiple sources. In dac, pages 247252. ACM, 2000.
[23] J. Wood, D.E. Root, and N.B. Tufillaro. A behavioral modeling approach to nonlinear
model-order reduction for RF/microwave ICs and systems. Microwave Theory and
Techniques, IEEE Transactions on, 52(9):22742284, 2004.
[24] P. Li and L.T. Pileggi. NORM: compact model order reduction of weakly nonlinear
systems. 2003.
[25] J.R. Phillips. Projection-based approaches for model reduction of weakly nonlinear,
time-varying systems. Computer-Aided Design of Integrated Circuits and Systems,
IEEE Transactions on, 22(2):171187, 2003.
[26] X. Li, P. Li, Y. Xu, and L.T. Pileggi. Analog and RF circuit macromodels for systemlevel analysis. 2003.
[27] A. Zhu and T.J. Brazil. Behavioral modeling of RF power amplifiers based on pruned
Volterra series. Microwave and Wireless Components Letters, IEEE, 14(12):563565,
2004.
[28] J. Roychowdhury. Reduced-order modeling of time-varying systems. Circuits and
Systems II: Analog and Digital Signal Processing, IEEE Transactions on, 46(10):1273
1288, 1999.
[29] A. Demir, A. Mehrotra, and J. Roychowdhury. Phase noise in oscillators: A unifying
theory and numerical methods for characterization. Circuits and Systems I: Fundamental Theory and Applications, IEEE Transactions on, 47(5):655674, 2000.
[30] X. Lai and J. Roychowdhury. TP-PPV: piecewise nonlinear, time-shifted oscillator
macromodel extraction for fast, accurate PLL simulation. In Proceedings of the 2006
IEEE/ACM international conference on Computer-aided design, pages 269274. ACM,
2006.
[31] X. Lai and J. Roychowdhury. Macromodelling oscillators using Krylov-subspace methods. In Proceedings of the 2006 Asia and South Pacific Design Automation Conference,
pages 527532. IEEE Press, 2006.
[32] B. Gu, K.K. Gullapalli, S. Hamm, B. Mulvaney, X. Lai, and J. Roychowdhury. Implementing nonlinear oscillator macromodels using Verilog-AMS for accurate prediction
of injection locking behaviors of oscillators. In Behavioral Modeling and Simulation
164
Workshop, 2005. BMAS 2005. Proceedings of the 2005 IEEE International, pages 43
47. IEEE, 2005.
[33] J.S. Han, E.B. Rudnyi, and J.G. Korvink. Efficient optimization of transient dynamic
problems in MEMS devices using model order reduction. Journal of Micromechanics
and Microengineering, 15:822, 2005.
[34] J. Lienemann, E.B. Rudnyi, and J.G. Korvink. MST MEMS model order reduction:
Requirements and benchmarks. Linear algebra and its applications, 415(2-3):469498,
2006.
[35] E. Rudnyi and J. Korvink. Model order reduction for large scale engineering models
developed in ANSYS. In Applied Parallel Computing, pages 349356. Springer, 2006.
[36] J. Chen and S.M. Kang. Model-order reduction of nonlinear mems devices through
arclength-based Karhunen-Loeve decomposition. In Circuits and Systems, 2001. ISCAS 2001. The 2001 IEEE International Symposium on, volume 3, pages 457460.
IEEE, 2001.
[37] E.B. Rudnyi and J.G. Korvink. Review: Automatic Model Reduction for Transient
Simulation of MEMS-based Devices. Sensors Update, 11(1):333, 2002.
[38] YC Liang, WZ Lin, HP Lee, SP Lim, KH Lee, and H. Sun. Proper orthogonal decomposition and its applications-part II: Model reduction for MEMS dynamical analysis.
Journal of Sound and Vibration, 256(3):515532, 2002.
[39] P. Benner, V. Mehrmann, and D.C. Sorensen. Dimension reduction of large-scale
systems. 2005.
[40] K. Willcox and J. Peraire. Balanced model reduction via the proper orthogonal decomposition. AIAA journal, 40(11):23232330, 2002.
[41] P.V. Kokotovic et al. Singular perturbations and order reduction in control theoryan
overview. Automatica, 12(2):123132, 1976.
[42] M. Bordemann and J. Hoppe. The Dynamics of Relativistic Membranes I: Reduction
to 2-dimensional Fluid Dynamics. Arxiv preprint hep-th/9307036, 1993.
[43] C.W. Rowley, T. Colonius, and R.M. Murray. Model reduction for compressible flows
using POD and Galerkin projection. Physica D: Nonlinear Phenomena, 189(1-2):115
129, 2004.
[44] B.F. Farrell and P.J. Ioannou. Accurate low-dimensional approximation of the linear
dynamics of fluid flow. Journal of the atmospheric sciences, 58(18):27712789, 2001.
165
[45] X. Ma and G.E. Karniadakis. A low-dimensional model for simulating threedimensional cylinder flow. Journal of Fluid Mechanics, 458(-1):181190, 2002.
[46] B.R. Noack, M. Schlegel, M. Morzy
nski, and G. Tadmor. System reduction strategy
for Galerkin models of fluid flows. International Journal for Numerical Methods in
Fluids, 63(2):231248, 2010.
[47] I.K. Fodor. A survey of dimension reduction techniques. Center for Applied Scientific
Computing, Lawrence Livermore National Laboratory, 2002.
[48] K.R. Schneider and T. Wilhelm. Model reduction by extended quasi-steady-state
approximation. Journal of mathematical biology, 40(5):443450, 2000.
[49] C.W. Rowley, I.G. Kevrekidis, J.E. Marsden, and K. Lust. Reduction and reconstruction for self-similar dynamical systems. Nonlinearity, 16:1257, 2003.
[50] D. Lebiedz, V. Reinhardt, and J. Kammerer. Novel trajectory based concepts for model
and complexity reduction in (bio) chemical kinetics. Model Reduction and CoarseGraining Approaches for Multiscale Phenomena, pages 343364.
[51] S. Dan, M.F. Madsen, H. Schmidt, and G. Cedersund. Reduction of a biochemical
model with preservation of its basic dynamic properties. FEBS Journal, 273(21):4862
4877, 2006.
[52] T.B. Kepler, LF Abbott, and E. Marder. Reduction of conductance-based neuron
models. Biological Cybernetics, 66(5):381387, 1992.
[53] D. Golomb, J. Guckenheimer, and S. Gueron. Reduction of a channel-based model for
a stomatogastric ganglion LP neuron. Biological cybernetics, 69(2):129137, 1993.
[54] W.M. Kistler, W. Gerstner, and J.L. Hemmen. Reduction of the Hodgkin-Huxley
equations to a single-variable threshold model. Neural Computation, 9(5):10151045,
1997.
[55] E.M. Izhikevich. Simple model of spiking neurons. Neural Networks, IEEE Transactions on, 14(6):15691572, 2003.
[56] D. Hansel, G. Mato, and C. Meunier. Phase dynamics for weakly coupled HodgkinHuxley neurons. EPL (Europhysics Letters), 23:367, 1993.
[57] Y. Maeda, K. Pakdaman, T. Nomura, S. Doi, and S. Sato. Reduction of a model for
an Onchidium pacemaker neuron. Biological cybernetics, 78(4):265276, 1998.
[58] B. Cartling. Response characteristics of a low-dimensional model neuron. Neural
computation, 8(8):16431652, 1996.
166
[59] K. Doya and A.I. Selverston. Dimension reduction of biological neuron models by
artificial neural networks. Neural Computation, 6(4):696717, 1994.
[60] S. Peles, B. Munsky, and M. Khammash. Reduction and solution of the chemical
master equation using time scale separation and finite state projection. The Journal
of chemical physics, 125:204104, 2006.
[61] L. Petzold and W. Zhu. Model reduction for chemical kinetics: An optimization
approach. AIChE journal, 45(4):869886, 1999.
[62] M.S. Okino and M.L. Mavrovouniotis. Simplification of mathematical models of chemical reaction systems. Chemical reviews, 98(2):391408, 1998.
[63] D. Lebiedz. Computing minimal entropy production trajectories: An approach to
model reduction in chemical kinetics. The Journal of chemical physics, 120:6890, 2004.
[64] S. Lall, J.E. Marsden, and S. Glavaski. A subspace approach to balanced truncation
for model reduction of nonlinear control systems. International Journal of Robust and
Nonlinear Control, 12(6):519535, 2002.
[65] U. Maas and S.B. Pope. Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space. Combustion and Flame, 88(3-4):239264, 1992.
[66] Z. Ren, S.B. Pope, A. Vladimirsky, and J.M. Guckenheimer. The invariant constrained
equilibrium edge preimage curve method for the dimension reduction of chemical kinetics. The Journal of chemical physics, 124:114111, 2006.
[67] E.A. Mastny, E.L. Haseltine, and J.B. Rawlings. Two classes of quasi-steady-state
model reductions for stochastic kinetics. The Journal of chemical physics, 127:094106,
2007.
[68] G. Kerschen, J.C. Golinval, A.F. Vakakis, and L.A. Bergman. The method of proper
orthogonal decomposition for dynamical characterization and order reduction of mechanical systems: an overview. Nonlinear Dynamics, 41(1):147169, 2005.
[69] S. Lall, P. Krysl, and J.E. Marsden. Structure-preserving model reduction for mechanical systems. Physica D: Nonlinear Phenomena, 184(1-4):304318, 2003.
[70] P. Koutsovasilis and M. Beitelschmidt. Comparison of model reduction techniques for
large mechanical systems. Multibody System Dynamics, 20(2):111128, 2008.
[71] V. Duindam and S. Stramigioli. Energy-based model-reduction of nonholonomic mechanical systems. In Robotics and Automation, 2004. Proceedings. ICRA04. 2004
IEEE International Conference on, volume 5, pages 45844589. IEEE, 2004.
167
[72] T. Kailath. Linear systems, volume 1. Prentice-Hall Englewood Cliffs, NJ, 1980.
[73] S. Sastry. Nonlinear systems: analysis, stability, and control. Springer Verlag, 1999.
[74] A. Odabasioglu, M. Celik, and L. T. Pileggi. Prima: passive reduced-order interconnect
macromodeling algorithm. 17(8):645654, 1998.
[75] C. Gu and J. Roychowdhury. ManiMOR: Model Reduction via Projection onto Nonlinear Manifolds, with Applications to Analog Circuits and Biochemical Systems. In
Proceedings of the IEEE International Conference on Computer Aided Design, 10-13
Nov. 2008.
[76] C. Gu. QLMOR: a New Projection-Based Approach for Nonlinear Model Order Reduction. In Proceedings of the IEEE International Conference on Computer Aided
Design, 2-5 Nov. 2009.
[77] C. Gu and J. Roychowdhury. Generalized Nonlinear Timing/Phase Macromodeling:
Theory, Numerical Methods and Applications. In Proceedings of the IEEE International Conference on Computer Aided Design, 7-11 Nov. 2010.
[78] C. Gu and J. Roychowdhury. FSM Model Abstraction for Analog/Mixed-Signal Circuits by Learning from I/O Trajectories. Proc. IEEE Asia and South Pacific Design
Automation Conference, 2011.
[79] S. Gupta, B. H. Krogh, and R. A. Rutenbar. Towards formal verification of analog
designs. In Proc. ICCAD-2004 Computer Aided Design IEEE/ACM Int. Conf, pages
210217, 2004.
[80] S. Blouin. Finite-state machine abstractions of continuous systems.
Queens University, 2003.
PhD thesis,
[81] P. Ye, E. Entcheva, R. Grosu, and S.A. Smolka. Efficient modeling of excitable cells
using hybrid automata. In Proc. of CMSB, pages 216227. Citeseer, 2005.
[82] E. Asarin, T. Dang, and A. Girard. Hybridization methods for the analysis of nonlinear
systems. Acta Informatica, 43(7):451476, 2007.
[83] R. Alur, T.A. Henzinger, G. Lafferriere, and G.J. Pappas. Discrete abstractions of
hybrid systems. Proceedings of the IEEE, 88(7):971984, 2000.
[84] J. Raisch and S.D. OYoung. Discrete approximation and supervisory control of continuous systems. Automatic Control, IEEE Transactions on, 43(4):569573, 1998.
[85] M.T. Heath. Scientific computing. McGraw-Hill, 1997.
168
[86] R.K. Brayton. Logic minimization algorithms for VLSI synthesis. Springer, 1984.
[87] S. Devadas, A. Ghosh, and K. Keutzer. Logic synthesis. McGraw-Hill, Inc. New York,
NY, USA, 1994.
[88] R. Harjani, R.A. Rutenbar, and L.R. Carley. OASYS: A framework for analog circuit
synthesis. Computer-Aided Design of Integrated Circuits and Systems, IEEE Transactions on, 8(12):12471266, 1989.
[89] E.J. Grimme. Krylov Projection Methods for Model Reduction. PhD thesis, University
of Illinois, EE Dept, Urbana-Champaign, 1997.
[90] B.N. Bond. Stability-Preserving Model Reduction for Linear and Nonlinear Systems
Arising in Analog Circuit Applications. PhD thesis, Massachusetts Institute of Technology, 2010.
[91] B.N. Bond and L. Daniel. Guaranteed stable projection-based model reduction for
indefinite and unstable linear systems. In Proceedings of the 2008 IEEE/ACM International Conference on Computer-Aided Design, pages 728735. IEEE Press, 2008.
[92] S.K. Tiwary, A. Gupta, J.R. Phillips, C. Pinello, and R. Zlatanovici. fSpice: a boolean
satisfiability based approach for formally verifying analog circuits. In Proc. Workshop
on Formal Verification of Analog Circuits (FAC), Princeton, 2008.
[93] S.K. Tiwary, A. Gupta, J.R. Phillips, C. Pinello, and R. Zlatanovici. First steps
towards SAT-based formal analog verification. In Proceedings of the 2009 International
Conference on Computer-Aided Design, pages 18. ACM, 2009.
[94] A. Singh and P. Li. On behavioral model equivalence checking for large analog/mixed
signal systems. In Computer-Aided Design (ICCAD), 2010 IEEE/ACM International
Conference on, pages 5561. IEEE.
[95] David C. Walter. Verification of Analog and Mixed-Signal Circuits Using Symbolic
Methods. PhD thesis, University of Utah, ECE Dept, 2007.
[96] Scott Little. Efficient Modeling and Verification of Analog/Mixed-Signal Circuits Using
Labeled Hybrid Petri Nets. PhD thesis, University of Utah, ECE Dept, 2008.
[97] L. Ljung and E.J. Ljung. System identification: theory for the user, volume 280.
Prentice-Hall Upper Saddle River, NJ, 1987.
[98] Rolf Isermann and Marco Mnchhof. Identification of Dynamic Systems: An Introduction with Applications. Springer, 2011.
[99] G.A. Baker and P.R. Graves-Morris. Pade approximants. Cambridge Univ Pr, 1996.
169
[100] L.T. Pillage and R.A. Rohrer. Asymptotic waveform evaluation for timing analysis.
IEEE Transactions on Computer-Aided Design, 9:352366, April 1990.
[101] P. Feldmann and R.W. Freund. Efficient linear circuit analysis by Pade approximation
via the Lanczos process. Computer-Aided Design of Integrated Circuits and Systems,
IEEE Transactions on, 14(5):639649, 1995.
[102] G.H. Golub and C.F. Van Loan. Matrix computations. Johns Hopkins Univ Pr, 1996.
[103] Y. Saad. Iterative methods for sparse linear systems. Society for Industrial Mathematics, 2003.
[104] Bruce Moore. Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduction. IEEE Trans. Automatic Control, 26:1732, February
1981.
[105] W. Rugh. Nonlinear System Theory - The Volterra-Wiener Approach. Johns Hopkins
Univ Press, 1981.
[106] P. Li and L.T. Pileggi. Compact reduced-order modeling of weakly nonlinear analog
and rf circuits. Computer-Aided Design of Integrated Circuits and Systems, IEEE
Transactions on, 24(2):184203, 2005.
[107] Joel R. Phillips. Projection-Based Approaches for Model Reduction of Weakly Nonlinear Time-Varying Systems. IEEE Transactions on Computer-Aided Design, 22(2):171
187, 2003.
[108] L. Pillage. Electronic circuit and system simulation methods. McGraw-Hill Professional, 1998.
[109] J. Roychowdhury. Numerical Simulation and Modelling of Electronic and Biochemical
R in Electronic Design Automation, 3(2-3):97303,
Systems. Foundations and Trends
2008.
[110] D. Garfinkel, L. Garfinkel, M. Pring, S.B. Green, and B. Chance. Computer applications to biochemical kinetics. Annual Review of Biochemistry, 39(1):473498, 1970.
[111] P.A. Absil, R. Mahony, and R. Sepulchre. Optimization algorithms on matrix manifolds. Princeton Univ Pr, 2008.
[112] J.M. Lee. Introduction to smooth manifolds. Springer Verlag, 2003.
[113] M. Spivak and M. Spivak. A comprehensive introduction to differential geometry.
1979.
170
[114] Dmitry Vasilyev, Michal Rewienski, and Jacob White. A TBR-based Trajectory
Piecewise-Linear Algorithm for Generating Accurate Low-order Models for Nonlinear
Analog Circuits and MEMS. Proceedings of the IEEE Design Automation Conference,
2003.
[115] J. Savoj and B. Razavi. A 10-Gb/s CMOS clock and data recovery circuit with a
half-rate linear phase detector. IEEE Journal of Solid-State Circuits, 36:761768, May
2001.
[116] Ning Dong and Jaijeet Roychowdhury. Piecewise Polynomial Nonlinear Model Reduction. Proceedings of the IEEE Design Automation Conference, 2003.
[117] Ning Dong and Jaijeet Roychowdhury. Automated Nonlinear Macromodelling of Output Buffers for High-Speed Digital Applications. Proceedings of the IEEE Design
Automation Conference, 2005.
[118] BSIM3. http://www-device.eecs.berkeley.edu/bsim3.
[119] B. Kofahl and E. Klipp. Modelling the dynamics of the yeast pheromone pathway.
Yeast, 21:831850, 2004.
[120] T.M. Thomson and D. Endy. Rapid Characterization of Cellular Pathways Using
Time-Varying Signals. In Sixth International Conference on Systems Biology, Oct
2005.
[121] J. Roychowdhury.
Reduced-order modelling of time-varying systems.
Trans. Ckts. Syst. II: Sig. Proc., 46(10), November 1999.
IEEE
[122] T. Sakurai and A.R. Newton. Alpha-power law MOSFET model and its applications
to CMOS inverter delay and other formulas. Solid-State Circuits, IEEE Journal of,
25(2):584594, 1990.
[123] G. Gildenblat, X. Li, W. Wu, H. Wang, A. Jha, R. van Langevelde, G.D.J. Smit, A.J.
Scholten, and D.B.M. Klaassen. PSP: An advanced surface-potential-based MOSFET
model for circuit simulation. Electron Devices, IEEE Transactions on, 53(9):1979
1993, 2006.
[124] J. Shawe-Taylor and N. Cristianini. Kernel methods for pattern analysis. Cambridge
Univ Pr, 2004.
[125] H. Shichman and D.A. Hodges. Modeling and simulation of insulated-gate field-effect
transistor switching circuits. Solid-State Circuits, IEEE Journal of, 3(3):285289, 1968.
[126] D.L. Boley.
Krylov space methods on state-space control models.
Trans. Ckts. Syst. II: Sig. Proc., 13(6):733758, 1994.
IEEE
171
[127] J. Phillips. Model Reduction of Time-Varying Linear Systems Using Approximate Multipoint Krylov-Subspace Projectors. In Proceedings of the IEEE International Conference on Computer Aided Design, November 1998.
[128] G. Hachtel, R. Brayton, and F. Gustavson. The sparse tableau approach to network
analysis and design. Circuit Theory, IEEE Transactions on, 18(1):101113, 1971.
[129] E.M. Izhikevich. Dynamical systems in neuroscience: The geometry of excitability and
bursting. The MIT press, 2007.
[130] A. Demir. Floquet theory and non-linear perturbation analysis for oscillators with
differential-algebraic equations. International journal of circuit theory and applications,
28(2):163185, 2000.
[131] E.A. Coddington and N. Levinson. Theory of ordinary differential equations. Tata
McGraw-Hill, 1972.
[132] K.S. Kundert, J.K. White, and A. Sangiovanni-Vincentelli. Steady-state methods for
simulating analog and microwave circuits. Kluwer Academic Publishers, 1990.
[133] A. Demir and J. Roychowdhury. A reliable and efficient procedure for oscillator PPV
computation, with phase noise macromodeling applications. Computer-Aided Design
of Integrated Circuits and Systems, IEEE Transactions on, 22(2):188197, 2003.
[134] K. Tsumoto, H. Kitajima, T. Yoshinaga, K. Aihara, and H. Kawakami. Bifurcations
in Morris-Lecar neuron model. Neurocomputing, 69(4-6):293316, 2006.
[135] D. Angluin. Learning regular sets from queries and counterexamples* 1. Information
and computation, 75(2):87106, 1987.
[136] K. Li, R. Groz, and M. Shahbaz. Integration testing of components guided by incremental state machine learning. In Testing: Academic and Industrial Conference-Practice
And Research Techniques, 2006. TAIC PART 2006. Proceedings, pages 5970. IEEE,
2006.
[137] M. De Berg, O. Cheong, M. Van Kreveld, and M. Overmars. Computational geometry:
algorithms and applications. Springer-Verlag New York Inc, 2008.
[138] G.J. Pappas. Hybrid systems: Computation and abstraction. PhD thesis, UNIVERSITY of CALIFORNIA, 1998.
[139] B. Aspvall. A linear-time algorithm for testing the truth of certain quantified boolean
formulas. Inform. Process. Lett., 8:121123, 1979.
172
[140] S. Prajna, A. Papachristodoulou, P. Seiler, and P. A. Parrilo. SOSTOOLS: Sum of
squares optimization toolbox for MATLAB, 2004.
[141] B. law Jakubczyk. Introduction to Geometric Nonlinear Control; Controllability and
Lie Bracket.
[142] Michal Rewienski. A Trajectory Piecewise-Linear Approach to Model Order Reduction
of Nonlinear Dynamical Systems. PhD thesis, MIT, 2003.
173
Appendix A
Introduction of Differential Geometry
The following introduction is based on notes in [141].
A.1
(A.1)
(A.2)
Then starting from a point p, the solution of (A.2), is denoted by a map (t, p) 7 t (p). The
family t of local maps are called the flow of the vector field f . The following properties
hold for t (p):
t1 t2 = t1 +t2 , t = (t )1 , 0 = id.
(A.3)
If t, p, t (p) are well-defined, the f is called complete, and its flow forms a oneparameter group of diffeomorphisms of X. On the other hand, given any t (p), we can
define a unique vector field f (p) by
(A.4)
f (p) = t (p).
t t=0
We will denote the flow of f by tf or exp(tf ).
174
A.2
The Lie bracket of two vector fields is another vector field which, roughly speaking,
measures the non-commutativeness of the flows of both vector fields. Non-commutativeness
means the dependence of the result of applying these flows.
There are several definitions (interpretations) of the Lie bracket:
Given a coordinate system, the Lie bracket of two vector fields f and g, denoted by [f, g]
is
g
f
[f, g](x) =
(x)f (x)
(x)g(x).
(A.5)
x
x
We will call this the Jacobian definition of Lie bracket.
Proposition A.2.1 The Lie bracket of vector fields f and g is equal identically to zero if
and only if their flows commute, i.e.,
[f, g] 0
tf sg (p) = sg tf (p),
s, t RN , p X.
(A.6)
Two vector fields having the property in (A.6) are called commuting.
g
f
Proposition A.2.2 Consider a curve (t) = t
t
tg tf (p). Then we have that its
first derivative at zero vanishes, (0) = 0, and the second derivative is given by the Lie
bracket (0) = 2[f, g](p).
A.3
.
x
We also write
(A.8)
Note that the coordinate change p = (q) transforms the differential equations p = f (p)
to q = Ad f (q). We also have Ad (1 f1 + 2 f2 ) = 1 Ad (f1 ) + 2 Ad (f2 ), and Ad =
Ad Ad (f ).
Proposition A.3.1 Consider the vector field Ad (f ). The local flow of this vector field is
t = t 1 .
(A.9)
175
Geometric definition of Lie bracket: we define the Lie bracket of f and g as the
derivative with respect t, at t = 0, of the vector field g transformed by the flow of f , i.e., we
define
f
[f, g](p) = Dt
(tf (p))g(tf (p)) = (Ad f g)(p).
(A.10)
t
t
t
Proposition A.3.2 If is a diffeomorphism of X, then
[Ad f, Ad g] = Ad [f, g].
(A.11)
In the following, we will also use adf g for [f, g], therefore adf is a linear operator in the
space of vector fields V (X). We also denote
ad0f g = g,
(A.12)
Proposition A.3.3 If the vector fields f and g are real analytic, then we have the following
expansion formula for the vector field g transformed by the flow of the vector field f :
(Ad f g)(p) =
t
X
(t)i
i=0
i!
(adif g)(p),
(A.13)
A.4
A smooth vector field f on X defines a linear operator Lf on the space of smooth functions
C (X) in the following way
n
X
f
fi (p)
(p).
(Lf )(p) = (t (p)) =
t t=0
xi
i=1
(A.14)
n
X
i=1
ai (x)
,
xi
(A.15)
and it defines a unique vector field defined by f = (a1 (x), , an (x))T . Therefore, there is a
one-to-one correspondence between the vector field f and the differential operator Lf .
Define the commutator of Lf and Lg by
[Lf , Lg ] = Lf Lg Lg Lf .
(A.16)
176
We have
[Lf , Lg ] = L[f,g]
(A.17)
The Lie bracket also satisfies two important properties, antisymmetry, and Jacobi
identity:
[f, g] = [g, f ]
[f, [g, h]] + [g, [h, f ]] + [h, [f, g]] = 0
(A.18)
The linear space V (X) of smooth vector fields f and g on X, with the Lie bracket as
product, is called the Lie algebra of vector fields on X.
A.5
F = {fu }uU
(A.19)
We shall see these points locally equivalent at the points p and p, respectively, if there
p 7 p which transforms the vector fields fu
is a local analytic diffeomorphism : X X,
into fu locally, i.e., Ad fu = fu , u U locally around p.
Denote by L and L the Lie algebras of vector fields generated by the families of F and
F.
A family of vector fields is called transitive at a point if its Lie algebra is of full rank
at this point, i.e., the vector fields in this Lie algebra span the whole tangent space at this
point.
We will use the notation
f[u1 , ,un ] = [fu1 , [fu2 , [ , [fuk1 , fuk ], , ]]
(A.20)
(A.21)
The distribution is called of class C if, locally around each point in X, there is a family
of vector fields {f } (called local generators of ) which spans , i.e., (p) = span f (p).
is called locally generated if the above family of vector fields is finite.
Finally, the distribution is called of dimension k, if dim(p) = k for all points in X.
177
Definition A.5.2 A vector field f belongs to a distribution , i.e., f , if f (p) (p)
for all p in X.
A distribution is called involutive if for any vector fields f, g , the Lie bracket is
also in , i.e., [f, g] .
Theorem A.5.1 (Frobenius theorem (local)) If is an involutive distribution of class
C and of dimension k on X, then locally around any point in X, there exists a smooth
change of coordinates which transforms the distribution to the following constant distribution span{e1 , , ek }, where ei = (0, , 0, 1, 0, , 0)T with 1 at the i-th place.
Proof The proof consists of two steps.
Firstly, we show the distribution is locally generated by k pairwise commuting vector
fields. We fix a point p in X and let F = (f1 , , fk ) RN k be any vector fields which
generate the distribution in a neighborhood of p. Multiplying F from the right by an
invertible k k matrix of smooth functions does not change the distribution spanned by the
columns of F (it only changes the generators). By a possible permutation of variables we can
make the upper k k matrix of F to be identity matrix. Therefore, the first k coefficients
of the Lie bracket [fi , fj ] are 0. On the other end, from involutivity, it follows that the Lie
bracket is a linear combination of columns of F , and therefore, the Lie bracket has to be 0.
This shows that the new vector fields commute.
Secondly, assume that the vector fields f1 , , fk generate the distribution locally
around p, and they commute. We can then choose fk+1 , , fN so that f1 , , fN are
linearly independent at p. Define a map by
(t1 , , tn ) exp(t1 f1 ) exp(t2 f2 ) exp(tN fN )(p).
(A.22)
A.6
178
Definition A.6.2 We call a subset S X an immersed sub-manifold of X of dimension
k if
[
S=
Si , where S1 S2 S,
(A.24)
i=1
Proposition A.6.1 If two vector fields f, g are tangent to an sub-manifold S, then also
their Lie bracket [f, g] is tangent to this sub-manifold.
Definition A.6.3 A foliation {S }A of X of dimension k is a partition
[
X=
S
(A.25)
of X into disjoint connected sub-manifolds S , called leaves, which has the following property: for any x X there exists a neighborhood U of x and a diffeomorphism : U V
RN onto an open subset V such that
\
N
((U
S )cc ) = {x = (x1 , , xN ) V |xk+1 = ck+1
(A.26)
, , xN = c },
where Pcc denotes a connected component of the set P .
A.7
(A.27)
179
Definition A.7.2 Let be the smallest distribution on X which contains the vector fields
in the family F (i.e., fu (p) (p) for all u U) and is invariant under any flow tu , u U
, ie,
Dtu (p)(p) (tu (p))
(A.28)
for all p X, u U.
Equivalently, we have
g = Adtu g ,
u U, t R.
(A.29)
Theorem A.7.1 (Orbit Theorem) Each orbit S = Orb(p) of a family of vector fields
F = {fu }uU is an immersed sub-manifold. Moreover, the tangent space to this sub-manifold
is given by the distribution , i.e.,
Tp S = (p),
p X.
(A.30)
Corollary A.7.2 If the vector fields fu are analytic, then the tangent space to the orbit can
be computed as
Tp X = L(p) = {g(p)|g Lie{fu }uU },
(A.31)
where Lie{fu }uU denotes smallest family of vector fields which contains the family F and is
closed under taking linear combinations and Lie bracket. In the smooth case, the following
inclusion holds
L(p) (p)
(A.32)
A.8
(A.33)
180
Appendix B
Benchmark Examples
This chapter lists a few benchmark examples we have tried. Table B.1 summarizes the
MATLAB scripts/codes for each benchmark example. The scripts are to be available on
http://www.eecs.berkeley.edu/~gcj.
Benchmark
NTL1
NTL2
INVC
NEURON ML
NEURON FN
LATCH1
B.1
The nonlinear transmission line circuit, shown in Fig. B.1 [2], is an example that has been
tested in lots of MOR literature. All resistors and capacitors are set to 1 and the diode I-V
characteristic is iD = e40vD 1. The input is the current source i = u(t); the output is the
voltage at node 1.
The modified nodal equations for this circuit are
v1 = 2v1 + v2 + 2 e40v1 e40(v1 v2 ) + u(t),
vN = vN + vN 1 1 + e40(vN1 vN ) .
2 i N 1,
(B.1)
181
N-1
......
B.2
The circuit was tested in papers on the TPWL method [2] [142]. The circuit diagram is
shown in Fig. B.2, where the resistances and capacitances are all set to 1, and the nonlinear
resistor at each stage has the I V characteristic in (v) = sgn(v)v 2. The input is the current
source i = u(t); the output is the voltage at node 1.
The differential equations are
v1 = 2v1 + v2 sgn(v1 )v12 + u(t),
vi = 2vi + vi1 + vi+1 sgn(vi )vi2 , 2 i N 1,
2
vN = vN + vN 1 sgn(vN )vN
.
1
N-1
(B.2)
......
182
B.3
The inverter chain circuit is shown in Fig. B.3. The resistances and capacitances are all
set to 1. The inverter Vout Vin characteristic is
Vout = finv (Vin ) = Vdd tanh(AVin ),
(B.3)
where Vdd is the supply voltage and A is a parameter. We choose Vdd = 1 and A = 5. The
input is the voltage source u(t) and the output is the node voltage at node N.
The differential equations are
x1 + x1 u(t) = 0
xi + xi finv (xi1 ) = 0,
1
(B.4)
2 i N.
3
......
B.4
The Morris-Lecar model [134] is one neuron model that describes potential dynamics of
firing neurons. The input is the injection current I.
The differential equations are
1
d
V =
(gL (V VL ) gCa M (V VCa ) gK N(V VK ) + I)
dt
CM
d
1
N =(N N) ,
dt
N
(B.5)
183
B.5
The FitzHugh-Nagumo model is a PDE neuron model. The input is the injection current
i0 (t) which is included in the boundary condition.
Let x [0, L], t 0, the PDEs are
2
1
v(x, t) = 2 v(x, t) + (f (v(x, t)) w(x, t) + c)
t
x
(B.6)
(B.7)
with
boundary conditions
v(0, t) = i0 (t),
x
v(L, t) = 0
x
(B.8)
w(x, 0) = 0,
(B.9)
= 0.015,
b = 0.5,
= 2,
c = 0.05.
(B.10)
wi (t) = w(ih),
i = 1, , N
(B.11)
vi (t) =
+ (f (vi (t)) wi (t) + c)) , i = 2, , N 1,
2
t
h
1 v2 v1
1
v1 (t) =
+ i0 (t) + (f (v1 (t)) w1 (t) + c)) ,
t
h
h
vN vN 1
1
1
(B.12)
184
B.6
The circuit diagram of the latch is shown in Fig. B.4. The latch is composed of a sampler
(consisting of a multiplexer) and a regenerator (consisting of two cross-coupled inverters).
Each output of the multiplexer and inverters is followed by an RC circuit that is not shown
in Fig. B.4. The latch samples the DAT A input when CLK is 1. When CLK is 0, two
inverters are cross-coupled, and re-generate the output (QBB).
Q
QB
INV
CLK
MUX
DATA
QBB
INV
dx1
+ (x1 fM U X (u1 , u2 , x3 )) =0,
dt
dx2
+ (x2 fIN V (x1 )) =0,
RC
dt
dx3
RC
+ (x3 fIN V (x2 )) =0,
dt
(B.13)
(B.14)
185
F
C
INV
MUX
B A