Towards an improved Eigensystem Realization Algorithm for low-error guarantees
Mohammad N. Murshed, Moajjem Hossain Chowdhury, Md. Nazmul Islam Shuzan, M. Monir Uddin
TTowards an improved Eigensystem RealizationAlgorithm for low-error guarantees
Mohammad N. Murshed ∗ , Moajjem Hossain Chowdhury † , Md. Nazmul Islam Shuzan ‡ , andM. Monir Uddin § Department of Mathematics and PhysicsNorth South UniversityDhaka, Bangladesh
June 4, 2020
Abstract
Eigensystem Realization Algorithm (ERA) is a tool that can produce a reduced order model (ROM)from just input-output data of a given system. ERA creates the ROM while keeping the number ofinternal states to a minimum level. This was first implemented by Juang and Pappa (1984) to analyzethe vibration of aerospace structures from impulse response. We reviewed ERA and tested it on singleinput single output (SISO) system as well as on multiple input single output (MISO) system. ERAprediction agreed with the actual data. Unlike other model reduction techniques (Balanced truncation,balanced proper orthogonal decomposition), ERA works just as fine without the need of the adjointsystem, that makes ERA a promising, completely data-driven, thrifty model reduction method. Inthis work, we propose a modified Eigensystem Realization Algorithm that relies upon an optimallychosen time resolution for the output used and also checks for good performance through frequencyanalysis. Four examples are discussed: the first two confirm the model generating ability and the lasttwo illustrate its capability to produce a low-dimensional model (for a large scale system) that is muchmore accurate than the one produced by the traditional ERA.
Introduction
We consider the discrete linear system, x i +1 = Ax i + Bu i (1) y i = Cx i + Du i (2)where x ∈ R n contains the internal states, u ∈ R p is the input, y ∈ R q refers to the output, and i is thetime index. A set of inputs can cause a system to react in a particular manner via the internal states ∗ Lecturer † Research Assistant ‡ Research Assistant § Associate Professor a r X i v : . [ phy s i c s . d a t a - a n ] M a y o result in the output. A ∈ R n × n , B ∈ R n × p , C ∈ R q × n and D ∈ R n × p are called the realizations.The number of equations in the system appears to be large for most of the common, complex prob-lems around us. An example is the Navier Stokes equation where the system is large owing to thehigh number of spatial nodes in the domain. It is very difficult to extract knowledge from such largesystems. Our focus in this work is on Eigensystem Realization Algorithm that makes a model fromjust impulse response data. We begin by reviewing some of the existing powerful model reductiontechniques.Model Order Reduction (MOR) is a way of reducing the complexity of models by means of pro-jection. MOR aids in creating a low-dimensional version of the large scale system and enables a goodenough understanding of the phenomenon in terms of fewer dominant states.This notion has beenused to analyze random generated linear systems, flow past a flat plate, combustion and many otherproblems of interest.Proper Orthogonal Decomposition (POD) is a statistical method to derive a low rank version of aset of data, [1]. The idea can be tailored to obtain a reduced order model, but research has been inprogress to find out better projections than the orthogonal ones. An interesting study of POD in thefield of turbulence is available in [2]. Balanced truncation (Moore, 1981) and balanced proper orthog-onal decomposition (BPOD) are two powerful model reduction methods successfully implemented onCFD problems and randomly generated systems [3]. Rowley et. al. used the idea of model reductionand extended it to solve non-linear complex compressible flows [4, 5]. The available techniques rely onthe direct system and a transformed system also known as the adjoint system.Eigensystem Realization Algorithm, [6, 7], is the only model reduction tool that is based just onthe direct system, hence, making it applicable on experimental data. It leverages just the input andoutput measurements to create a reduced order model for a given problem. The connection betweenERA and BPOD is shown in [8]. ERA has been tested on many occasions. It works well to makelow-order models for unstable flows, [9], which are then used to design controllers. Tangential inter-polation based eigensystem realization algorithm (TERA), [10], built on ideas of ERA to handle thehuge amount of input-output data for multi-input multi-output (MIMO) system. TERA is appliedon mass spring damper system and a cooling model for steel. A modified ERA, [8], is also developedand compared with the performance of balanced POD on the flow past a flat plate at a low Reynoldsnumber.Since measured data may contain noise, the way the data gets separated into the signal and thenoise has also been a subject under study. A noteworthy work can be found in [11] that discusses theeffect of noise, if any, on the modal parameters for a system.This paper is about the development of an improved version of Eigensystem Realization Algorithmthat monitors and empirically determines the rank and the time resolution of the output measure-ment to produce a reduced order model that is much more reasonable than the one from conventionalEigensystem Realization Algorithm. Background
In this section, we review the Eigensystem Realization Algorithm, its derivation and how it relatesto Dynamic Mode Decomposition. These will be used to elaborate on the modified Eigensystem2nput/Output MeasurementHankel (H) and time shifted Hankel Matrix (H’)SVD of H A r , B r , C r Figure 1: Eigensystem Realization AlgorithmRealization Algorithm in the next section.
Eigensystem Realization Algorithm
Eigensystem Realization Algorithm is a system identification method that was first used to createmodels for vibration in aerospace structures. This tool borrows from the idea of Ho’s algorithm, [12],to find the realizations while keeping the number of internal states to a minimum that is to say thatkeeping the dimension of matrix A as low as possible. Completely data-driven, ERA uses only impulseresponse of the system i.e. just the inputs and the outputs, [13].The discrete, linear time-invariant system in Eq.(1) and Eq.(2) can be excited by a pulse definedas u = (cid:40) k = 00 k > . For x = 0, the system reduces to x = Ax + B (1) to give x = B . We can iterate through the systemto get, x = Ax + Bu = ABx = Ax + Bu = A Bx = Ax + Bu = A B and so on, while the outputs appear to be y = Cx = 0 y = Cx = CBy = Cx = CABy = Cx = CA B y k = CA k − B which are also known as Markov parameters. Note thatthe dimension of Markov parameter is q × p . These are then used to construct the Hankel matrix andthe time shifted Hankel matrix, H = y y y y ... y m − s − y y y y ... y m − s ... ... ... ... . . . ... y s − ... ... ... ... y m − H (cid:48) = y y y y y m − s − y m − s y y y y y m − s y m − s +1 ... ... ... ... . . . ... y s ... ... ... y m − y m − . After performing the singular value decomposition of H = U Σ V ∗ , the truncated version of U , V , Σare computed as: U r = U (1 : r, :) V r = V (1 : r, :)Σ r = Σ(1 : r, r ) , where r is the rank of H . The reduced order model is then given by, A r = Σ − / U ∗ r H (cid:48) V r Σ − / B r = the first p columns of Σ / r V ∗ r C r = the first q rows of U r Σ / r . There are a few important points about the notation used. The output measurement is a functionof time and variable s controls the way we stack the time shifted output measurement. The Hankelmatrices above refer to a single input and single output system (SISO). The dimensions would changeas we deal with a different system e.g. MISO or MIMO. Note on Hankel Singular Values.
In general, eigenvalues give hint on the system stability. But,Hankel singular values identify the highly energetic states that contribute the most to characterize thesystem. That means the states with low energy can be truncated to obtain an approximate model.The Hankel singular values are computed from the SVD of the product of the controllability and theobservability Gramian.
Connection to Dynamic Mode Decomposition
Dynamic Mode Decomposition (DMD) is a strategy of creating a model from time series data, [14].It can be viewed as a special case of Koopman operator which is a way of representing a non-lineardynamical system as a infinite-dimensional linear system. DMD has numerous applications in variousdisciplines like fluid dynamics, neuroscience, epidemiology and many others. There are a lot of variantsof DMD that are to be utilized based on the type of data. For instance, time delay coordinate basedDMD is an option when the data is highly oscillatory [15].DMD aims to map the current states to the future states as, X i +1 = AX i . This equation from DMD is essentially Eq.(1) with no control. This implies that DMD is related toERA from a dynamical systems point of view. 4 odified Eigensystem Realization Algorithm
We propose an improved version of Eigensystem Realization Algorithm that searches for the optimalnumber of temporal nodes, N t , used to express the output from the actual system and also the optimalrank, r , used in ERA to attain a reduced order model that can predict the output with high accuracy.The idea behind this modified ERA is to run the conventional ERA multiple times so to identifythe best possible number of temporal nodes and the rank that keeps the error, (cid:15) = || y actual − y ERA || (3)as low as possible. N t controls the time resolution: certain N t values are optimal while others canyield large error. The steps in the middle are the same as the ones in the traditional ERA. At thevery end of this modified version, the A r , B r , and C r are computed based on the optimal rank. Theroutine is provided in Algorithm 1. The specialty of this updated version of ERA is that it uses theappropriate time resolution and identifies the ’best’ possible rank to keep the error to a minimum. Algorithm 1 Modified Eigensystem Realization AlgorithmPre-run ERA to identify N t and r that keep the error, (3), as low as possible Utilize the output measurements, y actual , based on the optimal time resolution, T /N t Construct the Hankel matrix ( H ) and the time-shifted Hankel matrix ( H (cid:48) )Compute the SVD of the Hankel matrix, H = U Σ V ∗ Find the truncated U, Σ, and V using the rank ( r ) from the pre-runCalculate the reduced system matrices just as in traditional ERAGenerate the output from ERA, y ERA , via the reduced order modelWe also recommend that a frequency analysis be performed after this modified ERA scheme isenacted. Frequency analysis is often helpful for engineering purposes. A way to do it is by using tf estimate on MATlab. This function takes in the input and output to generate an approximatetransfer function for a certain range of frequencies. Note that bodeplot on MATlab results in themagnitude and phase of the the system, but visual comparison is well done via tf estimate . Numerical Results
We have tested ERA on four different problems. The first two aim to stress on the model identificationfunction of ERA and the last two prove the ability of ERA to work as a model reduction tool. Theresults are generated on a personal computer (HP Pavilion 14) with CORE i5 8th Gen processor1.6-3.4 GHz and RAM of 8 GB via MATLAB version 2019b.
Example 1. Pitch Model (SISO)
The 3 D motion of an aircraft is governed by the pitch, plunge and surge models, [16]. Many statevariables come into play. Velocity, density, temperature and pressure are a few of them. Computingall these state variables in a grid is not easy since there may be spatial nodes as many as 10 . Inaeronautics, engineers care a lot about what is called the lift coefficient per unit span, C L = 2 LρU ∞ c where L is the lift force on the wing, ρ the air density and U ∞ free-stream velocity,and c the chord. Theangle of attack, α , is the angle between the airfoil chord and the flow direction. It can be thought of as5 a) Angular Acceleration Variation (b) Lift coefficient data and ERA prediction for pitchmodel Figure 2: Pitch model Input and Responsesome angular displacement which automatically makes ¨ α the angular acceleration. The pitch motionof an aircraft is the one that is observed with the nose moving up or down. We use the linearizedpitch model from [16] that reads: ddt x α = A 0 B ˙ α x α ˙ α + ¨ αC L = (cid:2) C C α C ˙ α (cid:3) x α ˙ α C ¨ α ¨ α where x is a vector containing the states. A non-dimensionalized version of time is used, τ = t U ∞ c .The pitfall in applying ERA on a state space model is that we are bound to use an impulseresponse. Thus, it is worth correctly identifying the right input for the pitch motion of an aircraft.We examine the behavior of various state variables over time, available in [16]. The angle of attackand angular velocity vary in the form of a ramp and step function, respectively, whereas the angularacceleration, Figure 2(a), follows a dirac delta function which is the requirement of ERA. Hence,angular acceleration would be a suitable input. Note that we consider the input from τ = 5 upto τ = 7 .
5. In this example, we use exponential functions to approximate the lift coefficient behavior, C l = 1105 ( 1( σ (cid:112) (2 π )) e − . τ − /σ ) + 2 . . τ − . ) , with σ set at √ . τ = 5. Example 2. Second Order State Space Model (MISO)
We also test ERA on a multi input single output system and compare the reduced order model to theoriginal model by feeding arbitrary inputs (step, ramp, unit and random). The MISO, in consideration,6 a) Input: Unit Step function (b) Input: Ramp function(c) Input: Random numbers
Figure 3: Arbitrary inputs used to check performance7 a) Input:Dirac Delta (b) Input: Unit Step function(c) Input: Ramp function (d) Input: Random numbers
Figure 4: Response comparison between ERA and actual system for various inputs8s defined as per the following matrices, A = (cid:20) − . − . . (cid:21) ; B = (cid:20) −
10 2 (cid:21) C = (cid:2) . . (cid:3) ; D = (cid:2) (cid:3) . The reduced order model extracted by ERA came out to be,ˆ A = (cid:20) . . − .
004 0 . (cid:21) ; ˆ B = (cid:20) − . − . − . . (cid:21) ˆ C = (cid:2) − . . (cid:3) ; ˆ D = (cid:2) . . (cid:3) . After extracting the Markov parameters and creating a model, we tested the model with differentinputs. The outputs, Figure 4, from these different inputs, Figure 3, are compared with outputsfrom the original model. The error between these outputs were minimal. Even when the input wasa randomly generated signal, the model was able to capture the characteristics of the signal. TheERA output resembles the actual output despite the difference in the Markov parameters of the twosystems. We can see that ERA found the parameters from the data.
Example 3. Heat Diffusion Equation (Sparse Model)
In this example, the heat diffusion equation for a rod of unit length (1D) is defined as, ∂T ( x, t ) ∂t = α ∂ T∂x + u ( x, t ) T (0 , t ) = T (1 , t ) = 0 T ( x,
0) = 0 , where x refers to space and t is the time. The equations above consist of the partial differentialequation showing the evolution of the temperature, T ( x, t ), the boundary conditions, and the initialcondition. The input is controlled by u ( x, t ).We utilized the system matrices available in [17] (the state matrix shown in Figure 5(a)) to com-pute the output for this problem. The output is then fed into ERA to identify a reduced order model.The crux is to observe the behavior of the norm of the error with the rank set in ERA as in Figure5(d). This allows us to find the optimal number of temporal nodes to be used that maintains a lowrank for the system. We have identified that N t = 300 works fine for r = 6. Thus, the state matrixin the original system is reduced from 200 ×
200 to 6 × lsim function, also agrees well with the output from the actual large system, Figure 5(b). Theerror over time is also displayed in Figure 5(c). The frequency analysis done by tf estimate shows that r greater 1 yields a model as good as the actual system, Figure 5(e). Example 4. Atmospheric Storm Track (Dense Model)
The atmospheric storm track is a model from oceanography used to analyze the velocity of the airflowin the zonal (latitude wise) and meridional (longitude wise) setting. We can imagine this of a flow ina channel, the physical domain of which is defined as,0 < x < π a) The sparsity pattern in state matrix, A (b) Comparison of the output from the actual systemand ERA(c) Uncertainty with time (d) Error behavior with the number of temporalnodes(e) Frequency analysis for different rank used in theERA Figure 5: Heat-Diffusion model reduction10 . π < y < . π < z < . In the z axis, z = 0 is the ground level and z = 1 is the tropopause.The mean velocity is set to vary with the altitude, U ( z ) = 0 . z, and time is non-dimensionalized as T = LU where L = 1000 km and U = 30 m/s . The system isthought to have a uniform flow, but can be disturbed by a linear damping at the entrance and theexit of the track. The details of the dynamics can be found in [17]. The governing equation of thestates is, dψdt = Aψ, (4)where ψ is the velocity variable.The output from the actual system is put into ERA. Figure 6(c) shows that r ≈
55 for N t = 500 and N t = 750 whereas use of 300 temporal nodes allows for r = 30. Setting the rank to 30 and numberof temporal nodes to 300 results in ERA output that resembles the original output, 6(a). The normof the difference between the approximate output and the actual output is also plotted in 6(b). Weobserve that the intractable original system of dimension 598 ×
598 gets reduced to a 30 ×
30 systemby ERA.The transfer function estimate for several different rank values are plotted in Figure 6(d). The ERAmodel with r = 1 is far away from the actual system, r = 5 and r = 10 show improvement and r = 15enables reduced order modeling that is as efficient as the actual system. Thus, proper selection of thetime resolution and rank along with a frequency analysis in the Eigensystem Realization Algorithmshows promise in building reduced order models with low error. Conclusion and Future Work
In this work, we delineated the steps in our proposed modified Eigensystem Realization Algorithm andimplemented this method on four test problems. Modified ERA identifies the model in the first twoexamples and performs as a tool to reduce the order of the model in the third and fourth example. Bymodel identification, we mean finding the state, input and output matrix and model order reductionrefers to the minimization of the size of the state matrix, also known as the system matrix. The outputpredicted by ERA agree well with that from the original system. The first example is concerned aboutthe pitching motion of an aircraft, and the second one a second order state space model. The thirdexample is the heat-diffusion equation and the last one a model for the airflow velocity when a stormor a cyclone surges. Indeed, the third and the fourth numerical tests demonstrate that the rank shouldbe carefully set at or above 5 to minimize error in the output predicted by ERA and also to get atransfer function estimate that is much close to that of the original system.We plan to work on a survey of all the model order reduction techniques and apply them on anarray of synthetic and practical data and finally weigh the pros and cons of each technique. At thesame time, our work would also be to establish any connection between model order reduction methodand DMD, [18]. 11 a) Comparison of the output from the actual systemand ERA (b) Absolute value of the difference in the output from theactual system and ERA(c) Error behavior with the number of temporalnodes (d) Frequency analysis for different rank used in theERA
Figure 6: Atmospheric Storm Track model reduction12 cknowledgement
This project is partially funded by Office of Research, North South University. (Grant number:
CTRG-19/SEPS/06 ) References [1] A. Chatterjee, “An introduction to the proper orthogonal decomposition,”
Current science , pp.808–817, 2000.[2] L. Sirovich, “Turbulence and the dynamics of coherent structures. i. coherent structures,”
Quar-terly of applied mathematics , vol. 45, no. 3, pp. 561–571, 1987.[3] K. Willcox and J. Peraire, “Balanced model reduction via the proper orthogonal decomposition,”
AIAA journal , vol. 40, no. 11, pp. 2323–2330, 2002.[4] C. W. Rowley, T. Colonius, and R. M. Murray, “Model reduction for compressible flows usingpod and galerkin projection,”
Physica D: Nonlinear Phenomena , vol. 189, no. 1-2, pp. 115–129,2004.[5] C. W. Rowley, “Model reduction for fluids, using balanced proper orthogonal decomposition,”
International Journal of Bifurcation and Chaos , vol. 15, no. 03, pp. 997–1013, 2005.[6] R. Pappa and J.-N. Juang, “Galileo spacecraft modal identification using an eigensystem real-ization algorithm,” in , 1984, p.1070.[7] J.-N. Juang and R. S. Pappa, “An eigensystem realization algorithm for modal parameter iden-tification and model reduction,”
Journal of guidance, control, and dynamics , vol. 8, no. 5, pp.620–627, 1985.[8] Z. Ma, S. Ahuja, and C. W. Rowley, “Reduced-order models for control of fluids using theeigensystem realization algorithm,”
Theoretical and Computational Fluid Dynamics , vol. 25, no.1-4, pp. 233–247, 2011.[9] T. L. Flinois and A. S. Morgans, “Feedback control of unstable flows: a direct modelling approachusing the eigensystem realisation algorithm,”
Journal of Fluid Mechanics , vol. 793, pp. 41–78,2016.[10] B. Kramer and S. Gugercin, “Tangential interpolation-based eigensystem realization algorithmfor mimo systems,”
Mathematical and Computer Modelling of Dynamical Systems , vol. 22, no. 4,pp. 282–306, 2016.[11] P. Li, S. Hu, and H. Li, “Noise issues of modal identification using eigensystem realization algo-rithm,”
Procedia engineering , vol. 14, pp. 1681–1689, 2011.[12] H. p. Zeiger and A. McEwen, “Approximate linear realizations of given dimension via ho’s algo-rithm,”
IEEE Transactions on Automatic Control , vol. 19, no. 2, pp. 153–153, 1974.[13] J. N. Kutz,
Data-driven modeling & scientific computation: methods for complex systems & bigdata . Oxford University Press, 2013.[14] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,”
Journal offluid mechanics , vol. 656, pp. 5–28, 2010. 1315] M. N. Murshed and M. Monir Uddin, “Time delay coordinate based dynamic mode decompositionof a compressible signal,” in , 2019, pp. 1–5.[16] S. L. Brunton, S. T. Dawson, and C. W. Rowley, “State-space model identification and feedbackcontrol of unsteady aerodynamic forces,”
Journal of Fluids and Structures , vol. 50, pp. 253–270,2014.[17] Y. Chahlaoui and P. Van Dooren, “A collection of benchmark examples for model reduction oflinear time invariant dynamical systems.” 2002.[18] J. H. Tu and C. W. Rowley, “An improved algorithm for balanced pod through an analytictreatment of impulse response tails,”
Journal of Computational Physics , vol. 231, no. 16, pp.5317–5333, 2012.
Appendix
Derivation of the reduced system matrices from ERA
Let’s consider a possible scenario where the Hankel matrix and the time shifted Hankel matrix aredefined as, H = y y y y y y y y y = CB CAB CA BCAB CA B CA BCA B CA B CA B = ¯ O ¯ CH (cid:48) = y y y y y y y y y = CAB CA B CA BCA B CA B CA BCA B CA B CA B = ¯ OA ¯ C. It is important to note that the Hankel matrices can also be written in terms of the observability andcontrollability. Controllability refers to how the inputs can excite the states and observability meanshow the states can affect the outputs. The singular value decomposition of the Hankel matrix reads, H = U Σ V T H = U T V T where Σ = T . T will then be used to define the observability and controllability as, H = ¯ O ¯ C = U T V T → ¯ O = U T, ¯ C = T V T . Finally, we construct the reduced system matrices H (cid:48) = ¯ OA ¯ CH (cid:48) = U T AT V T → ˆ A = T − U T H (cid:48) V T − . escriptor system for the SLICOT based problems The last two examples in this paper used systems from [17] which is a collection of benchmark problemsthat have real-life applications. The collection essentially gives the matrices (
A, B, C, D, E ) for differentdense, sparse and second order state space models. They are based on the following descriptor system, E x i +1 = A x i + B u i (5) y i = C x i + D u i , (6)where E is invertible. Additional information like Hankel singular values, frequency and fre-quency responsefre-quency response