L23 – Fast Time Domain Simulation via Decomposition Prof. Olivier de Weck Simulation of LTI Systems large order systems medium order systems small order systems number of stat e variables n 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 xt Axt But yt Cxt Dut ? =+ =+ () () () () () () nn nm p npm AB CD ×× ×× ∈∈ ∈∈ nullnull nullnull ? LTI systems are very common in engineering applications and simulating them is of great importance. Mathematical model: Problems of interest are usually complex! 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) HS w[n] S z[n] LTI plant discretized S w(t) w[n] z[n] ode45 method lsim.m method 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(n s 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(n s ) cost). ? Multiple sampling rates trade efficiency with accuracy. ? However, newlsim must first diagonalize the A matrix. Newlsim Flowchart Diagonalization Subsystem Planner DownsamplingDiscretization 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-7 data process mfiles original 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 23! AT d AT AT Ae IAT==++ + +L ? Block diagonal A matrix allows fast discretization. 1 2 3 1 2 00 exp 0 0 AT AT A T AT e AT e e ? ??? ?? ? ??? ?? = ? ? ? ? ?? ??? ? LL OO MOO MO ? 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 D DT plantAA filter Iu[n] y[n] 1. Frequency domain perspective (Input downsampling): LPF Downsampling SS Interpolation IDT plantLu[n] y[n] 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. 10 0 10 1 10 2 10 3 10 -8 10 -6 10 -4 10 -2 10 0 10 2 10 4 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 s 2 ) can be reduced to O(n s ). ? 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 Ver 2.2 modal Magellan RWA 6 inputs Starlight OPD #1 ~ 4 WFT 1 ~ 2 261 sec spent ? Using newlsim to simulate SIM model v2.2. Computation Time 10 0 10 1 10 2 10 3 10 -2 10 -1 10 0 10 1 10 2 10 3 n s 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) 934.4N/A949.31046.4T cpu [s]2000 1.334E-5N/A1.325E-51.346E-5σ 1500 382.7N/A393.3425.3T cpu [s] 5.228E-5N/A5.184E-55.452E-5σ 1000 0.68801.37500.07800.2960T cpu [s] 2.733E-5N/A2.714E-52.869E-5σ 96.0320N/A98.562104.25T cpu [s] 3.061E-5N/A3.002E-53.006E-5σ 15.2030N/A13.85915.3280T cpu [s]500 1.187E-51.187E-51.159E-51.167E-5σ 1.01504.53100.23400.6410T cpu [s]100 1.122E-61.722E-61.1713E-61.714E-6σ 50 newlsimlsimnewlyapfreqresultsn s Computation Cost Budget Main computation cost for lsim.m: Discretization State transitionOther Other 2( ) ss nmpnn×++×× 3 17.5742 s n× Diagonalization State transitionOther Other 2(2 ) s mp nn× ++×× 3 27.8276 s n× Main computation cost for new_lsim.m: ? Although diagonalization still takes O(n s 3 ) operations, it is generally faster than discretization. ? State transition of lsim is O(n s 2 n) but that of newlsim is O(n s n). 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 Ey [] [$ ] [] 22 2 ? Ey Ey Ey [] [ $ ] [] ? 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 Ey [ $ ] [] ? 2 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.