L23- Fast Time Domain simulation via decomposition Prof olivier de Weck imulation of lti systems LTI systems are very common in engineering applications and simulating them is of great importance small dynamic urge dynamic Mathematical l0000 range model 5000 ext generatio x(1)=Ax(t)+B() 2000 order Telescope Interest 1000 stems y(1)=Cx(1)+D()500 medium order A∈mB∈囗"m司 切100 stems Middeck Starlight Active origi C∈D∈日p Testbed rder MIT stems 3-mass Problems of interest are I single mdl usually complex oscillator 10 normalized dynamic bandwidth min(o,)
L23 – Fast Time Domain Simulation via Decomposition Prof. Olivier de Weck Simulation of LTI Systems large order systems medium order systems small order number of state variables n systems s 10000 5000 2000 1000 500 200 100 50 20 10 5 2 1 normalized dynamic bandwidth max( ) min( ) n n ω ω small dynamic range large dynamic range 0 10 2 10 4 10 6 10 8 10 Next Generation Space Telescope region of interest SIM v2.2 NEXUS Starlight Middeck Active Control Experimen t Origins Testbed (MIT) 3-mass model single DOF oscillator x t Ax t Bu t y t Cx t Du t • = + = + () () () () () () nn nm p n pm A B C D × × × × ∈ ∈ ∈ ∈ • LTI systems are very common in engineering applications and simulating them is of great importance. Mathematical model: Problems of interest are usually complex!
Obiective of Newlsim Total computation time=# of simulation X time per simulation Reduction of simulation time can push the fixed total time contour further to Utopia Are there better designs? · Newlsim is a runtime reduction Fixed total time budet simulation Tia=M 1pu scheme · Newlsim is an Goal is to push the tradeoff attempt to be the curve further towards Utopia equivalence of落 Isim. m for large Are the results current order systems believable? 三| state of the art Number of designs explored N(# of simulations) Simulation schemes Solving ode ivp 1. Standard Ode IVP routines like Runge Kutta LTI plant ode45 continuous Z method time 2. Discretization and propagation of state(lsim m) LTI plant continuous z(t) z[n]: Isim.m time method LTI plant [n] discretized zin
Objective of Newlsim Number of designs explored N (# of simulations) Are there better designs? Are the results believable? Utopia Goal is to push the tradeoff curve further towards Utopia current state of the art Fixed total time budet: Ttot= N*Tcpu • Total computation time = # of simulation X time per simulation. • Reduction of simulation time can push the fixed total time contour further to Utopia. • Newlsim is a runtime reduction simulation scheme. • Newlsim is an attempt to be the equivalence of lsim.m for large order systems. Simulation Schemes LTI plant continuous time Solving ODE IVP 1. Standard ODE IVP routines like Runge Kutta 2. Discretization and propagation of state (lsim.m) w(t) z(t) LTI plant continuous time w(t) z(t) S H w[n] S z[n] LTI plant discretized S w(t) w[n] z[n] ode45 method lsim.m method
sim and newlin State transition is faster than ODE solvers for LTI systems Large time step stability is guaranteed Discretization is slow for large order systems(e.g. SIM v2. 2 has 2184 state variables) Block diagonal a matrix allows fast discretization Memory is not enough for the big state matrix required by Ititr (14Gb for SIM simulation) Subsystem simulations are less memory consuming State transition cost goes with O(n 2), even worse for large-order ystems Lsim. m does not recognize special structure(e.g. sparse matrix) Simulink discrete state space(DSS)solver exploits matrices sparsity (i.e. O(n)cost) Multiple sampling rates trade efficiency with accuracy However newlsim must first diagonalize the a matrix Newlsim flowchart x Original System(ABCD) Diagonalization md n modal system↓ Subsystem seg_plan_xm subsystems PI anner subsystem bandwidth build- multirate sys.m Discretization Downsampling:Original Input DT LTI downsampled input DT subsystems ] System Solver st_sim.m process subsystem responses Interpolation interp. m Superposition mfiles Final response
Lsim and Newlsim • State transition is faster than ODE solvers for LTI systems. • Large time step stability is guaranteed. • Discretization is slow for large order systems (e.g. SIM v2.2 has 2184 state variables). • Memory is not enough for the big state matrix required by ltitr (14Gb for SIM simulation). • State transition cost goes with O(ns 2), even worse for large-order systems. • Lsim.m does not recognize special structure (e.g. sparse matrix). • Block diagonal A matrix allows fast discretization. • Subsystem simulations are less memory consuming. • Simulink discrete state space (DSS) solver exploits matrices sparsity (i.e.O(ns) cost). • Multiple sampling rates trade efficiency with accuracy. • However, newlsim must first diagonalize the A matrix. Newlsim Flowchart Diagonalization Subsystem Planner Discretization Downsampling DT LTI System Solver Interpolation Superposition Original System (ABCD) Final Response md.m seg_plan_x.m build_multirate_sys.m st_sim.m interp.m Original Input 0 20 40 60 80 100 120 140 160 180 200 0 20 40 60 80 100 120 140 160 180 200 nz = 292 modal system subsystems subsystem bandwidth DT subsystems downsampled input 0 50 100 150 200 250 -25 -20 -15 -10 -5 0 5 10 15 20 25 subsystem responses 0 50 100 150 200 250 -0.5 0 0.5 1 1.5 2 2.5 3 x 10 data -7 process mfiles original
Subsystem Planning Determine how to break the total system into subsystems and assign appropriate sampling rates Calculation using SF init(max SM=I dSF down sample factor SM start mode EM end mode TM total mode bz block size SF= DSF /2 or EM-SM+1≥BZ DSF= DSF EM-SM≥B EM=SM+ BZ max-I or TM SM=EM+I EM= EM or TM SM=EM+I DSF= DSF/2 or Fast discretization Discretization is to compute matrix exponential I+at+ 2 Block diagonal A A,T 0 0 matrix allows fast exp 0 A,T discretization Generic algorithm for fast discretization for i= 1 to total number of subsystems index= location of subsystem; DT_ system(index)=C2d(subsystem) end Computation time for matrix exponential of a 1000 by 1000 diagonal matrix tells the story dense matrix: 27 11 seconds(on the test machine) sparse matrix 0.031 seconds(on the test machine
Subsystem Planning DSF down sample factor SM start mode EM end mode TM total mode BZ block size DSF_init (max) SM = 1 EM = ? EM ≥ TM? EM – SM + 1 ≥ BZ_min? DSF = DSF EM – SM ≥ BZ_max? EM = EM or TM SM = EM + 1 DSF = DSF / 2 or 1 EM = SM + BZ_max - 1 or TM SM = EM + 1 DSF = DSF / 2 or 1 Calculation using bisection method Y N N N Y Y Terminate • Determine how to break the total system into subsystems and assign appropriate sampling rates. Fast Discretization • Discretization is to compute matrix exponential. 22 33 2 3! AT d AT AT A e I AT = =+ + + +L • Block diagonal A matrix allows fast discretization. 1 2 3 1 2 0 0 exp 0 0 A T A T A T AT e AT e e ⎛ ⎞ ⎡ ⎤ ⎡ ⎤ ⎜ ⎟ ⎢ ⎥ ⎢ ⎥ = ⎜ ⎟ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎝ ⎠ ⎣ ⎦ ⎣ ⎦ L L O O M OO M O • Generic algorithm for fast discretization for i = 1 to total number of subsystems index = location of subsystem; DT_system(index) = c2d(subsystem); end • Computation time for matrix exponential of a 1000 by 1000 diagonal matrix tells the story: • dense matrix: 27.11 seconds (on the test machine) • sparse matrix 0.031 seconds (on the test machine)
Downsampling Downsamp pling trades efficiency with accuracy. Two downsampling schemes are introduced I Frequency domain perspective(Input downsampling) u[n]-AA filterDH DT plant yIn LPF Downsampling Ss Interpolation 2. Time domain perspective(Output downsampling): L DT plant Lifting SS Interpolation large step transition small step transition tput(state)traj Interpolation Many interpolation schemes are available. Newlsim chooses MATLAB command interp. m Downsampled by 4 expe showing Sinusoidal signals are downsampled nd th Lineal inte pol with 1000 Hz original sampling rate f frequency domain method is fairly accurate
Downsampling u[n] y[n] AA filter D DT plant I 1. Frequency domain perspective (Input downsampling): LPF Downsampling SS Interpolation u[n] y[n] L DT plant I Lifting SS Interpolation 2. Time domain perspective (Output downsampling): • Downsampling trades efficiency with accuracy. Two downsampling schemes are introduced: output (state) trajectory small step transition large step transition Interpolation • Many interpolation schemes are available. Newlsim chooses MATLAB command interp.m. • A small experiment showing accuracy: • Sinusoidal signals are downsampled and then interpolated, with 1000 Hz original sampling rate. • Frequency domain method is fairly accurate. 100 101 102 103 10-8 10 -6 10-4 10 -2 100 10 2 104 f Linear intepolation Cubic spline Frequency domain interpolation Downsampled by 4
Fast state transition State transition of a discrete-time system is matrix-vector product Taking advantage of the block diagonal structure of the A matrix, O(n 2)can be reduced to O(ns) Newlsim recognizes the structure but lsim does not Computation time for a dense system omputation time for a diagonalized system O"." O Isit.m △DISs 8 Computation Example Using newlsim to simulate Sim model v2.2 SIM model 261 sec spent Magellan rWA 6 inputs ver 2.2 modal Starlight OPD #I-4 WET 1-2 SIM ver 2.2 MIMo result time in sec time in sec
Fast State Transition • State transition of a discrete-time system is matrix-vector product. • Taking advantage of the block diagonal structure of the A matrix, O(ns 2) can be reduced to O(ns). • Newlsim recognizes the structure but lsim does not. Computation Example 0 50 100 150 200 250 -8 -6 -4 -2 0 2 4 6 8 x 10-7 time in sec starlight opd 1~6 SIM ver 2.2 MIMO result 0 50 100 150 200 250 -25 -20 -15 -10 -5 0 5 10 15 20 25 time in sec SIM model Magellan RWA 6 inputs Ver 2.2 modal Starlight OPD #1 ~ 4 WFT 1 ~ 2 261 sec spent • Using newlsim to simulate SIM model v2.2
Computation Time Logarithmic scaled plot of computation time. Newlsim is more efficient for large order systems 灵10 ovel K newish Computation Time(1) n results freq newiyap Isim newish 50TcmS]02960 0.0780 1.3750 0.6880 1714E61.1713E-61.722E-61.122E6 100Tcs]0640023404531010150 a167E51159E517-51187E5 50Tcs]115.32801389 N/A 15.2030 a3006E3005NA306E5 1000811042598562 N/A 96.0320 5452E-55.184E5 N/A 5.228E-5 100Tcms14253 393.3 N/A 382.7 1346E-51.325E-5 N/A 1.334E-5 2000s]10464 949.3 N/A 934.4 2.869E-52714E-5 N/A 2.733E-5
Computation Time 100 101 102 10 3 10 -2 10 -1 10 0 10 1 102 10 3 ns lsim.m newlsim.m crossover small systems large systems lsim.m newlsim.m • Logarithmic scaled plot of computation time. Newlsim is more efficient for large order systems. Computation Time (1) 2000 Tcpu [s] 1046.4 949.3 N/A 934.4 σ 1.346E-5 1.325E-5 N/A 1.334E-5 1500 Tcpu [s] 425.3 393.3 N/A 382.7 σ 5.452E-5 5.184E-5 N/A 5.228E-5 1000 Tcpu [s] 0.2960 0.0780 1.3750 0.6880 σ 2.869E-5 2.714E-5 N/A 2.733E-5 T 104.25 98.562 N/A 96.0320 cpu [s] σ 3.006E-5 3.002E-5 N/A 3.061E-5 500 Tcpu [s] 15.3280 13.859 N/A 15.2030 σ 1.167E-5 1.159E-5 1.187E-5 1.187E-5 100 Tcpu [s] 0.6410 0.2340 4.5310 1.0150 σ 1.714E-6 1.1713E-6 1.722E-6 1.122E-6 50 ns results freq newlyap lsim newlsim
Computation Cost Budget Main computation cost for Isim. m Other Discretization State transition Other 175742×n232×(n,+m+p)×n,xn Main computation cost for new Isim.m Other Diagonalization*State transition Other 27.8276×n32×(2+m+p)×n2×n Although diagonalization still takes o(n )operations, it is generally faster than discretization State transition of lsim is O(n, n) but that of newlsim is O(n n) A ccuracy Time response(starlight OPD #1)of a 200-mode subsystem of SIM V2.2 model driven by Magellan RWA Fx data Relative averaged point to oint error is about 0.015% (002%)E[y- ELy Relative rms value error is about 0.017% time response by new Isim 0yi-√A Relative mean value error is about 0y]-E[ (<0.02 Ely -0
Computation Cost Budget Main computation cost for lsim.m: Other Discretization State transition Other 2( ) s s × ++ ×× n mp nn 3 17.5742× ns Other Diagonalization State transition Other 2 (2 ) × ++ ×× mp nn s 3 27.8276× ns Main computation cost for new_lsim.m: • Although diagonalization still takes O(ns 3) operations, it is generally faster than discretization. • State transition of lsim is O(ns 2n) but that of newlsim is O(nsn). Accuracy • Relative averaged point to point error is about 0.015% (< 0.02%) • Relative RMS value error is about 0.017% (< 0.02%) • Relative mean value error is about 0.017% (< 0.02%) Time response (starlight OPD #1) of a 200-mode subsystem of SIM V2.2 model driven by Magellan RWA Fx data Ey Ey E y [] [$ ] [ ] 2 2 2 − Ey Ey E y [] [$] [ ] − 0 1 2 3 4 5 6 7 8 -0.5 0 0.5 1 1.5 2 2.5 x 10-7 time response by direct lsim 0 1 2 3 4 5 6 7 8 -0.5 0 0.5 1 1.5 2 2.5 x 10-7 time response by new lsim Ey y E y [ $ ] [ ] − 2
Summary Newlsim is developed to extend simulation design capability to arge order L systems Newlsim is meant to be complementary to lsim Newlsim assumes diagonalizable a matrix and must diagonalize it before simulation Newlsim divides a large system into smaller subsystems Newlsim assigns multiple sampling rates Newlsim discretizes subsystems using O(n)discretization scheme Newlsim applies lifting to downsampling Newlsim employs o(n)state transition scheme Newlsim is notably more efficient than Isim for large order systems simulations
Summary • Newlsim is developed to extend simulation design capability to large order LTI systems. • Newlsim is meant to be complementary to lsim. • Newlsim assumes diagonalizable A matrix and must diagonalize it before simulation. • Newlsim divides a large system into smaller subsystems. • Newlsim assigns multiple sampling rates. • Newlsim discretizes subsystems using O(n) discretization scheme. • Newlsim applies lifting to downsampling. • Newlsim employs O(n) state transition scheme. • Newlsim is notably more efficient than lsim for large order systems simulations