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

1. Theoretical and practical aspects of linear and nonlinear model order reduction techniques Dmitry Vasilyev Thesis supervisor: Jacob K White December 19, 2007

2. Outline • 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

3. Main motivation for MOR: system-level simulation 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

4. Main motivation for MOR: system-level simulation Modern processor or system-on-a chip > Millions of transistors > Kilometers of interconnects> Linear and nonlinear devices Figures thanks to Mike Chou, Michał Rewienski

5. Model reduction problem inputs outputs Many (> 104) internal states inputs outputs few (<100) internal states • Reduction should be automatic • Must preserve input-output properties

6. Differential equation model • Original complex model: • Model can represent: • Finite-difference spatial discretization of PDEs • Circuits with linear capacitors and inductors Need accurate input-output behavior

7. Nonlinear model reduction problem • 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:

8. Linear 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

9. Transfer function of LTI system Transfer function: Laplace transform of the output Laplace transform of the input Matrix-valued rational function of s

10. Outline • 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

11. Linear MOR problem n – large (thousands)! q – small (tens) Need the reduction to be automatic and preserve input-output properties (transfer function)

12. Approximation error • Wide-band applications: model should have small worst-case error maximal difference over all frequencies ω

13. Approximation error • Narrow-band approximation: need good approximation only near particular frequency: Elmore delay: preserved if the first derivative at zero frequency is matched. ω frequency response

14. Linear MOR methods roadmap Linear MOR Projection-based Non projection-based Most widely used. Will be the central topic of this work.

15. Projection-based linear MOR • Pick projection matrices V and U: such that VTU=I x Uzx x n q U z Ax VTAUz

16. Projection-based linear MOR • Important: reduced system depends only on column spans of V and U • D = Dr, preserves response at infinite frequency

17. Linear MOR methods roadmap 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 )}

18. Linear MOR methods roadmap 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 )}

19. TBR idea 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

20. Balanced truncation reduction (TBR) 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

21. TBR benefits • Guaranteed stability • In practice provides more reduction than Krylov • H-infinity error bound => ideal for wide-band approximations Hankel singular values

22. Linear MOR methods roadmap 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]

23. Outline • 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

24. Nonlinear MOR framework • 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

25. Nonlinear MOR framework Problem: evaluation of V Tf(Uz)is still expensive Two solutions: Use Taylor series off Use TPWL approximation

26. Taylor series for nonlinear MOR 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.

27. Nonlinear MOR framework Problem: evaluation of V Tf(Uz)is still expensive Two solutions: Use Taylor series off Use TPWL approximation Will be discussed next

28. Trajectory piecewise linear (TPWL) approximation of f( ) [Rewieński, 2001] Training trajectory x0 x2 x1 … wi(x)is zero outside circle xs Simulating trajectory

29. VT Projection and TPWL approximation yields efficient f r( ) q x1 Air Ai q U = Air q n n Evaluating fTPWLr( ) requires only O(sq2) operations

30. TPWL approximation of f( ). Extraction algorithm • 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

31. Outline • 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

32. Krylov-subspace methods Balanced-truncation method The matter of this contribution What are projection options for TPWL? Can we use it? Used in the original work [Rewienski ‘02]

33. Example problem RLC line Linearized system has non-symmetric, indefinite Jacobian

34. Numerical results – nonlinear RLC transmission line System response for input current i(t) = (sin(2π/10)+1)/2 • Input: training input testing input Voltage at node 1 [V] Time [s]

35. Numerical results –RLC transmission line 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

36. Micromachined switch example Finite-difference model of order 880 non-symmetric indefinite Jacobian Model description [Hung ‘97]

37. TPWL-TBR results– MEMS switch example Errors in transient Unstable! Odd order models unstable! Even order models beat Krylov ||yr – y||2 Why??? Order of reduced system

38. Explanation of even-odd effect – Problem statement 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?

39. Perturbation behavior of TBR basis is similar to symmetric eigenvalue problem Eigenvectors ofM0 : Eigenvectors ofM0 + Δ: Mixing of eigenvectors (assuming small perturbations): ciklarge when λi0 ≈ λk0

40. Hankel singular values, MEMS beam example This is the key to the problem. Singular values are arranged in pairs! # of the Hankel singular value

41. Explaining even-odd behavior 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 Helps to ensure that linearizations are stable

42. Summary • 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.

43. Outline • 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

44. Balanced truncation reduction (TBR) 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

45. Most reduction algorithms effectively separately approximate dominant eigenspaces of Pand Q: • 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

46. RC line (symmetric circuit) V(t) – input i(t) - output • Symmetric Jacobian, B=CT, P=Qall controllable states are observable and vice versa

47. RLC line (nonsymmetric circuit) 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

48. Lightly damped RLC circuit 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

49. Xi= (PQ)Ui Ui+1= qr(Xi) “iterate” AISIAD model reduction algorithm Idea of AISIAD approximation: Approximate eigenvectors using power iterations: Uiconverges to dominant eigenvectors ofPQ Need to find the product (PQ)Ui How?

50. Approximation of the product Ui+1=qr(PQUi), AISIAD algorithm Vi≈ qr(QUi) Ui+1≈ qr(PVi) Approximate using solution of Sylvester equation Approximate using solution of Sylvester equation