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