Theoretical and practical aspects of linear and nonlinear model order reduction techniques

Download Presentation

Theoretical and practical aspects of linear and nonlinear model order reduction techniques

Loading in 2 Seconds...

- 69 Views
- Uploaded on
- Presentation posted in: General

Theoretical and practical aspects of linear and nonlinear model order reduction techniques

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Theoretical and practical aspects of linear and nonlinear model order reduction techniques

Dmitry Vasilyev

Thesis supervisor: Jacob K White

December 19, 2007

- Motivation
- Overview of existing methods:
- Linear MOR
- Nonlinear MOR

- TBR-based trajectory piecewise-linear method
- Modified AISIAD linear reduction method
- Graph-based model reduction for RC circuits
- Case study: microfluidic channel
- Conclusions

Nonlinear

Linear

Q: How to reduce the cost of simulating the big system?

System to simulate

Device 2

Device 1

A: Reduce the complexity of each sub-system, i.e.

approximate input-output behavior of the system by a system of lower complexity.

Device 3

…

…

Device 10

The goal of MOR in a nutshell

Modern processor or

system-on-a chip

> Millions of transistors

> Kilometers of interconnects> Linear and nonlinear devices

Figures thanks to Mike Chou, Michał Rewienski

inputs

outputs

Many (> 104)

internal states

inputs

outputs

few (<100)

internal states

- Reduction should be automatic
- Must preserve input-output properties

- Original complex model:

- Model can represent:
- Finite-difference spatial discretization of PDEs
- Circuits with linear capacitors and inductors

Need accurate input-output behavior

- Original complex model:

- Requirements for reduced model
- Want q << n (cost of simulation is q3)
- f r should be fast to compute
- Want yr(t) to be close to y(t)

- Reduced model:

- state

A – stable, n xn (large)

- vector of inputs

- vector of outputs

- Model can represent:
- Spatial discretization of linear PDEs
- Circuits with linear elements

Transfer function:

Laplace transform of the output

Laplace transform of the input

Matrix-valued rational function of s

- Motivation
- Overview of existing methods:
- Linear MOR
- Nonlinear MOR

- TBR-based trajectory piecewise-linear method
- Modified AISIAD linear reduction method
- Graph-based model reduction for RC circuits
- Case study: microfluidic channel
- Conclusions

n – large

(thousands)!

q – small (tens)

Need the reduction to be automatic and preserve

input-output properties (transfer function)

- Wide-band applications: model should have small worst-case error

maximal difference

over all frequencies

ω

- Narrow-band approximation: need good approximation only near particular frequency:

Elmore delay:

preserved if the first derivative at zero frequency is matched.

ω

frequency response

Linear MOR

Projection-based

Non projection-based

Most widely used.

Will be the central topic of this work.

- Pick projection matrices V and U:
such that VTU=I

x

Uzx

x

n

q

U

z

Ax

VTAUz

- Important: reduced system depends only on column spans of V and U
- D = Dr, preserves response at infinite frequency

Linear MOR

Projection-based

Non projection-based

V, U =

eig{PQ, QP}

Balancing-based (TBR)

V, U =

K((si I-A)-1,B), K((si I-A)-T,CT)

Krylov-subspace

methods

Proper Orthogonal

Decomposition methods

V=U = {x(t1)… x(tq )}

Linear MOR

Projection-based

Non projection-based

Will describe

next.

Balancing-based (TBR)

V, U =

K((si I-A)-1,B), K((si I-A)-T,CT)

Krylov-subspace

methods

Proper Orthogonal

Decomposition methods

V=U = {x(t1)… x(tq )}

u

y

LTI SYSTEM

t

t

input

output

P (controllability)

Which states are easier to reach?

Q (observability)

Which states produce more output?

X (state)

- Reduced model retains most controllable and most observable states
- Such states must be both very controllable and very observable

Compute controllability and observability gramians P and Q : (~n3)AP + PAT + BBT =0 ATQ + QA + CTC = 0

Reduced model keeps the dominant eigenspaces of PQ : (~n3)

PQui= λiuivTiPQ = λivTi

Reduced system: (VTAU, VTB, CU, D)

Very expensive. P and Q are dense even for sparse models

- Guaranteed stability
- In practice provides more reduction than Krylov
- H-infinity error bound => ideal for wide-band approximations

Hankel

singular values

Linear MOR

Projection-based

Non projection-based

Twice better error

bound than TBR

[Glover ’84]

Hankel-optimal MOR

Singular perturbation

approximation

Match at zero frequency

instead of infinity [Liu ‘89]

Transfer function

fitting methods

Promising topic

of ongoing research [Sou ‘05]

- Motivation
- Overview of existing methods:
- Linear MOR
- Nonlinear MOR

- TBR-based trajectory piecewise-linear method
- Modified AISIAD linear reduction method
- Graph-based model reduction for RC circuits
- Case study: microfluidic channel
- Conclusions

- Consider original (large) system:

- Projection of the nonlinear operator f(x):
substitute x ≈ Uz and project residual onto VT

Problem: evaluation of V Tf(Uz)is still expensive

Problem: evaluation of V Tf(Uz)is still expensive

Two solutions:

Use Taylor series off

Use TPWL approximation

Problem: evaluation of V Tf(Uz)is still expensive

Two solutions:

Use Taylor series off

Use TPWL approximation

- Accurate only near expansion point or weakly nonlinear systems
- Storing of dense tensors is expensive; limits the series to orders no more than 3.

Problem: evaluation of V Tf(Uz)is still expensive

Two solutions:

Use Taylor series off

Use TPWL approximation

Will be discussed next

Training trajectory

x0

x2

x1

…

wi(x)is zero outside circle

xs

Simulating trajectory

VT

q x1

Air

Ai

q

U

=

Air

q

n

n

Evaluating fTPWLr( ) requires only O(sq2) operations

- Compute A1
- Obtain V1 and U1using linear reduction for A1
- Simulate training input, collect and reduce linearizations Air = W1TAiV1f r (xi)=W1Tf(xi)

Initial system position

x0

x2

x1

…

xs

Training trajectory

Non-reduced state space

- Motivation
- Overview of existing methods:
- Linear MOR
- Nonlinear MOR

- TBR-based trajectory piecewise-linear method
- Modified AISIAD linear reduction method
- Graph-based model reduction for RC circuits
- Case study: microfluidic channel
- Conclusions

Krylov-subspace methods

Balanced-truncation method

What are projection options for TPWL?

Can we use it?

Used in the original work

[Rewienski ‘02]

RLC line

Linearized system has

non-symmetric, indefinite Jacobian

System response for input current i(t) = (sin(2π/10)+1)/2

- Input:

training

input

testing

input

Voltage at node 1 [V]

Time [s]

TBR-based TPWL beat

Krylov-based

4-th order TBR TPWL reaches the limit of TPWL representation

Error in transient

||yr – y||2

Order of the reduced model

Finite-difference

model of order 880

non-symmetric

indefinite Jacobian

Model description [Hung ‘97]

Errors in transient

Unstable!

Odd order models unstable!

Even order models beat Krylov

||yr – y||2

Why???

Order of reduced system

Consider two LTI systems:

Perturbed:

Initial:

TBR reduction

TBR reduction

~

Projection basis V

Projection basis V

Define our problem:

How perturbation in the initial system

affects TBR projection matrices?

Eigenvectors ofM0 :

Eigenvectors ofM0 + Δ:

Mixing of eigenvectors (assuming small perturbations):

ciklarge when λi0 ≈ λk0

This is the key to the problem.

Singular values are arranged in pairs!

# of the Hankel singular value

The closer Hankel singular

values lie to each other, the

more corresponding eigenvectors

of V tend to intermix!

- Analysis implies simple recipe for using TBR
- Pick reduced order to ensure that
- Remaining Hankel singular values are small enough
- The last kept and the first removed Hankel singular values are well separated

- Pick reduced order to ensure that

Helps to ensure that linearizations are stable

- We used TBR-based linear reduction procedure to generate TPWL reduced models
- Order reduced 5 times while maintaining comparable accuracy with Krylov TPWL method (efficiency improved 125 times!)
- Simple recipe found which helps to ensure stability.

- Motivation
- Overview of existing methods:
- Linear MOR
- Nonlinear MOR

- TBR-based trajectory piecewise-linear method
- Modified AISIAD linear reduction method
- Graph-based model reduction for RC circuits
- Case study: microfluidic channel
- Conclusions

Compute controllability and observability gramians P and Q : (~n3)AP + PAT + BBT =0 ATQ + QA + CTC = 0

Reduced model keeps the dominant eigenspaces of PQ : (~n3)

PQui= λiuivTiPQ = λivTi

Reduced system: (VTAU, VTB, CU, D)

Very expensive. P and Q are dense even for sparse models

- Arnoldi [Grimme ‘97]:U = colsp{A-1B, A-2B, …}, V=U, approx. Pdomonly
- Padé via Lanczos [Feldman and Freund ‘95]colsp(U) = {A-1B, A-2B, …}, - approx. Pdomcolsp(V) = {A-TCT, (A-T )2CT, …},- approx. Qdom
- Frequency domain POD [Willcox ‘02], Poor Man’s TBR [Phillips ‘04]

colsp(U) = {(jω1I-A)-1B, (jω2I-A)-1B, …}, - approx.Pdom

colsp(V) = {(jω1I-A)-TCT, (jω2I-A)-TCT, …},- approx.Qdom

However, what matters is the product PQ

V(t) – input

i(t) - output

- Symmetric Jacobian, B=CT,
P=Qall controllable states are observable and vice versa

Vector of states:

- P and Q are no longer equal!
- By keeping only mostly controllable and/or only mostly observable states, we may not find dominant eigenvectors of PQ

y(t) = i1

R = 0.008,

L = 10-5

C = 10-6

N=100

- Exact low-rank approximations of P and Q of order < 50 leads to PQ≈ 0

Xi= (PQ)Ui Ui+1= qr(Xi)

“iterate”

Idea of AISIAD approximation:

Approximate eigenvectors using power iterations:

Uiconverges to dominant eigenvectors ofPQ

Need to find the product (PQ)Ui

How?

Vi≈ qr(QUi)

Ui+1≈ qr(PVi)

Approximate using solution of Sylvester equation

Approximate using solution of Sylvester equation

Right-multiply by Vi

(original AISIAD)

X

H, qxq

X

M, nxq

Right-multiply by Vi

^

X

H, qxq

X

Approximate!

M, nxq

Right-multiply by Vi

^

X

H, qxq

X

Approximate!

M, nxq

We can take advantage of various

methods, which approximate P and Q

-M

X

X

A

=

+

H

qxq

nxq

nxn

Need only column span of X

Schur decomposition of H :

-M

X

X

A

~

~

~

=

+

~

Solve for columns of X

X

- Applicable to any stable A
- Requires solving q times

Schur decomposition of H :

Solution can be accelerated via fast MVP

Original method suggests IRA, needs A>0 [Zhou ‘02]

LR-sqrt

^

^

- Obtain low-rank approximations of Pand Q
- Solve AXi+XiH+ M = 0, =>Xi≈ PVi where H=ViTATVi, M = P(I - ViViT)ATVi + BBTVi
- Perform QR decomposition of Xi =UiR
- Solve ATYi+YiF+ N = 0, =>Yi≈ QUi where F=UiTAUi, N = Q(I - UiUiT)AUi + CTCUi
- Perform QR decomposition of Yi =Vi+1 Rto get new iterate.
- Go to step 2 and iterate.
- Bi-orthogonalize VandUand construct reduced model:

^

^

(VTAU, VTB, CU, D)

H-infinity norm of reduction error

(worst-case discrepancy over all frequencies)

N = 1000,

1 input

2 outputs

- Fast approximation to TBR
- Especially useful if gramians do not share common dominant eigenspace
- Improved accuracy and extended applicability over AISIAD
- Generalized to the systems in descriptor form

- Motivation
- Overview of existing methods:
- Linear MOR
- Nonlinear MOR

- TBR-based trajectory piecewise-linear method
- Modified AISIAD linear reduction method
- Graph-based model reduction for RC circuits
- Case study: microfluidic channel
- Conclusions

Reduction quality

Hankel-optimal

TBR

Krylov-subspace methods

Graph-based reduction:

manipulates RC network by removing nodes and inserting new elements

Cost of reduction

State-space model in the frequency domain:

vs

vm

Vector of node voltages (state):

Jm

vk

symmetric

Jk

External ports

Conductance matrix is analogous, ground node is excluded.

- Consider removing a single internal node (Nth), partition matrices and vectors:

- Substitute vN in the system equations (one step of Gaussian elimination):

- Where

- Problem: last capacitance term is negative! Potentially inserting a negative capacitor???
- The term was ignored in the TICER algorithm [Sheehan ‘99]. Leads to inconsistent diagonal update.

Added

conductance

Capacitance-like

- Claim: keeping the exact Taylor series is OK:

Gnew

Cnew

- Proof: Define projection:

Congruence transform

Model is always

stable and passive

- When is it safe to eliminate a node?

- Denominator expansion:

(used in TICER)

- Numerator term ~s2 (element-by-element)

(overlooked in TICER)

Using these criteria the reduced order will be chosen on-the-fly

- Given the initial circuit and maximal frequency of interest
- Using lowest-degree ordering (minimize fill-ins)
- Perform the elimination of the “qualified” nodes by inserting new capacitors and resistors:

(for every nodes i and j which were connected via the node N)

- Until no nodes satisfy elimination conditions.

Testing only substitution rules, 1-CDF of the reduction error

tested more than

30,000 circuits

Same elimination rules, same average reduction

different elimination criteria:

Narrower distribution

Better worst-case accuracy

- Improved accuracy and error control over TICER by using correct Taylor series and elimination criteria
- Preserves stability and passivity
- Generalized to parameter-dependent case
- Fastest, though conservative

- Motivation
- Overview of existing methods:
- Linear MOR
- Nonlinear MOR

- TBR-based trajectory piecewise-linear method
- Modified AISIAD linear reduction method
- Graph-based model reduction for RC circuits
- Case study: microfluidic channel
- Conclusions

- Inside the carrying fluid, marker fluid spreads governed by 3D convection-diffusion equation:

u(t) = Cin(t)

(input)

r1

(outputs)

y2(t)

w

- Using mapped-domain finite-difference volume discretization, obtained model has 2842 unknowns (large)

y1(t) =<Cout>(t)

y3(t)

V

- In case of constant mobility and diffusivity the model is linear:

Linear reduction techniques are extremely efficient for such models

[Vasilyev, Rewienski, White ‘06]

Modified AISIAD

runtime:

73s

TBR runtime:

2207s

(Matlab implementation)

For arbitrary nonlinearity in convection and diffusion coefficients and TPWL, this problem is verychallenging!

[Vasilyev, Rewienski, White ‘06]

However, the problem becomes more tractable, if one considers a quadratic problem

This is the case for affine μand D:

μ(C) = μ0C+ μ1

D(C) = D0C+ D1

- Affine mobility and diffusivity leads to quadratic model:

- Use orthogonal projection V = U, V TV = I

Reduced quadratic system

Projected reduced quadratic model of size 60 approximates original system of size 2842 quite well:

Krylov-subspace basis,

Quadratic reduction

- Performed applicability analysis of TBR-based TPWL models based on matrix perturbation theory
- Developed modified AISIAD method which is aimed at approximating TBR for the cases where gramians do not necessarily share common dominant eigenspaces
- Developed graph-based parameterized RC reduction method and improved nominal reduction

Prof. Jacob White – my supervisor,

Profs. Luca Daniel and Alexandre Megretski,

Profs. Karen Willcox, John Kassakian, John Wyatt, Dr. Yehuda Avniel, Dr. Joel Phillips, Dr. Mark Reichelt

My groupmates: Anne, Bo, Brad, Carlos, Dave, Jay, Jung Hoon, Homer, Kin, Laura, Lei, Michał, Shihhsien, Steve, Tarek, Tom, Xin, Zhenhai

My wife, Patrycja

Thank you! Спасибо!Grazie! Dziękuje!

Komapsumnida!Xie Xie! Dua Netjer en ek!

Generalized Lyapunov equations:

Lead to similar

approximate power iterations

Low-rank gramians

(cost varies)

mAISIAD

LR-square root

(inexpensive step)

(more expensive)

For the majority of non-symmetric cases,

mAISIAD works better than low-rank square root

H-infinity norm of reduction error

(worst-case discrepancy over all frequencies)

N = 1000,

1 input

2 outputs

Taken from Oberwolfach benchmark collection, N=1357 7 inputs, 6 outputs

For symmetric systems (A = AT, B = CT) P=Q, therefore

mAISIAD is equivalent to LRSQRT for P,Q of order q

^

^

RC line example

- Cost of the algorithm is directly proportional to the cost of solving a linear system:(where sjj is a complex number)
- Cost does not depend on the number of inputs and outputs

(non-descriptor case)

(descriptor case)

Top 5 eigenvectorsof Q

Top 5 eigenvectors of P

Union of eigenspaces of P and Q

does not necessarily approximate

dominant eigenspace of PQ .