# Modified Empirical Eigenfunctions and Its Applications for Model Reduction of Nonlinear Spatiotemporal Systems

Authors: Mian Jiang; Shuangqi Liu; Jigang Wu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Model reduction can greatly reduce complexity and difficulty of control design for spatiotemporal systems (STS) in engineering applications. Empirical eigenfunctions (EEFs) are widely used for the model reduction of spatiotemporal systems, however, truncation of higher modes may describe the behaviours of nonlinear spatiotemporal systems inaccurately. In this paper, modified EEFs are proposed and applied to model reduction of nonlinear spatiotemporal systems. Modified EEFs are obtained via modifying the weights matrix in the method of snapshots, which can be rewritten as linear combinations of initial EEFs. The coefficient matrix for combinations is computed according to the nonlinear temporal dynamics of STSs. Thus, the effects of higher modes are considered into modified EEFs with less computational requirements. The reduced model can give a more accurate description for behaviours of the system. The performance of the proposed method is further proved theoretically, and a numerical example demonstrates the effectiveness of the proposed method.


---

## Body

## 1. Introduction

With their nonuniformly distributed dynamics in space, many engineering problems belong to a class of nonlinear spatiotemporal systems (STSs) or distributed parameter systems (DPSs). Their infinite-dimensional spatiotemporal coupling and complex nonlinear behaviours make modelling, system analysis, numerical simulations and control design very difficult. Thus, model reduction for nonlinear STSs at a reasonable cost and accuracy has a great significance in practical engineering applications. Empirical eigenfunctions (EEFs) identified by proper orthogonal decomposition (POD) [1] are widely used as global spatial basis functions in advanced methods for model reduction of nonlinear STSs [2–8]. Modes obtained from POD can be calculated quite easily as the solutions of an eigenvalue problem involving second-order correlation tensors. When establishing a low-dimensional dynamic model, few POD modes are used to capture the dominant dynamics of nonlinear STSs according to the energetic optimality. However, the traditional POD is a linear dimension reduction (i.e., linear projection and linear reconstruction) method, and it produces a linear approximation to the measured data with nonlinear spatiotemporal structure, which may not guarantee the assumption that minor components do not contain important information. As a result, truncation of higher modes with small “energies” may have great influences on accuracy of the reduced model [9]. Without considering the effects of higher modes, the reduced dynamic model will give an inaccurate description for the spatiotemporal dynamic behaviours of STSs.This situation has been addressed in some literatures. With more computational power for training, the nonlinear dimension reduction method named nonlinear POD was developed for nonlinear problems to retain more information using fewer components. The nonlinear POD has a more powerful capability of dimension reduction than the traditional POD for nonlinear systems. This method has been widely used to deal with the nonlinear dimension reduction problems in the field of machine learning, image processing and pattern recognition. Examples include principal curves, multilayer autoassociative neural networks, kernel function approach, and radial basis function(RBF) networks. ISOMAP-based spatiotemporal modelling [10] and LLE- based nonlinear spatiotemporal modelling [11] are also proposed to improve the performance of POD-based reduced mode in recent years. However, these methods mainly focus on the reduction and analysis of multivariate time series or the order reduction of DPSs with no analytical models [12].To deal with the inaccuracy of traditional POD method, some literatures related to fluid dynamics are also given. In 1990s, Aubry [13] and Armbruster [14] respectively pointed out that the reduced models using the first several EEFs can have difficulties reproducing behaviours dominated by irregular transitions between different dynamical states. In 2003, Noack [15] had found that the low-dimensional Galerkin model by traditional POD method cannot approximate the transient behaviours of flow around a circular cylinder well. He introduced a new “shift mode” to improve the performance of the low-dimensional model, which incorporates the steady solution in the POD framework. In 2010, Sengupta [16] built the relationship of instability modes with POD modes in the study of global spatiotemporal nonlinear instabilities for flow past a cylinder. In recent years, the effects of higher modes are highlighted [17] and discussed [18] in the field of fluid dynamic systems. Some approaches such as spectral viscosity method [19] and variational multiscale model [20] based on the traditional POD are developed to deal with the problems. However, some additional empirical parameters are introduced into these methods for model correction, and their optimizations depend on specific problems. Another important method to improve the accuracy and dynamical system properties of POD-based reduced order models is sensitivity-based enhancement of modes [21–23], which uses sensitivity analysis (SA) to include the flow and shape parameters influence during the basis selection process to develop more robust reduced order models for varying parameters. It has shown promising results for parametric variation due to richness in the eigenfunctions and presents modified eigenfunctions though in a different setting.To introduce the effects of higher modes into the POD-based reduced model, nonlinear closure modelling [24, 25], nonlinear Galerkin method, and improved EEFs [26–28] are proposed to improve the model reduction performance. The higher modes contain relatively low energies but play a vital role in the overall dynamics of complex flows. Nonlinear closure modelling [24, 25] enhances the performance of POD-based reduced order model at low computational cost. Though it still remains a challenging task for 3D turbulent flows, closure modelling already shows promising results for Burgers’ equation and the Navier-Stokes equations. Nonlinear Galerkin method considers the effects of higher modes via approximate inertial manifolds (AIM) [26] and is used to approximate the nonlinear DPSs for fluids systems. The space spanned by EEFs is split into two subspaces and the interactions between lower and higher modes are built via AIM to compensate the modelling accuracy. However, the calculations for AIM in this method are theoretically complex and computational power-consuming. The improved EEFs [28] use the transformation for initial EEFs from traditional POD to derive a new set of basis functions, while the transformation matrix is obtained according to balanced truncation method for an approximated linear temporal system of a STS. This method introduces the dynamical information of higher modes into new basis functions to derive the reduced model. Because the transformation matrix is obtained from the linear approximation for nonlinear dynamics of the STS, this method can be improved by utilizing the global nonlinear dynamics directly to calculate the matrix.In this paper, modified EEFs are developed and applied to model reduction of nonlinear STSs. The EEFs are modified by introducing an extra weights matrix into the method of snapshots, which is commonly used to calculate the EEFs in traditional POD. This procedure transforms the derivation of modified EEFs to linear combinations of initial EEFs, where the coefficients matrix is computed directly according to the nonlinear temporal dynamics of STSs. Thus, the effects of higher modes are considered into modified EEFs and give a better accuracy of the reduced-order model. The effectiveness of the proposed method is proved theoretically. This method requires less computational consumption and a numerical example shows that it has a better performance than the improved EEFs based modelling [28].


## 2. EEFs Based Model Reduction for STSs

Assuming that a kind of nonlinear STSs can be governed by a partial differential equation (PDE) with the following state description:(1)∂X∂t=AX+BU+FX,∂X∂z,⋯,U,∂U∂z,⋯subject to a number of boundary and initial conditions. In (1), X=X(z,t) denotes the vector of state variable, where t∈[0,∞) is the time variable, z∈Ω is the spatial coordinate, and only one spatial-dimension is considered here. U=U(z,t) denotes the vector of manipulated spatio-temporal inputs, where U=∑i=1mhi(z)ui(t) and ui(t) is the ith temporal signal with certain spatial distribution hi(z). A and B are two linear operators that involve linear spatial derivatives on the state variable and spatiotemporal input. ∂X/∂z and ∂U/∂z denote the partial derivatives of X and U, respectively. F is a nonlinear function containing spatial derivatives for X and U. A scalar product is defined on spatial domain Ω introduced in the phase space of (1), which is given by (f(z),g(z))Ω=∫Ωf(z)g(z)dz.The spatiotemporal variableX(z,t) is assumed to be expanded onto the set of EEFs φ1(z),φ2(z),⋯,φN(z) (the computational details are given in Appendix) with corresponding temporal coefficients(2)Xz,t=∑i=1Nxitφiz.Suppose that the spatial-temporal outputY(z,t) is measured at M spatial locations. Substituting the expansion (2) into (1) and taking the inner product with φi(z) (orthogonality), the integration of equation will give a finite-dimensional nonlinear ordinary differential equation (ODE) system in a general form as follows:(3)x˙t=Axt+But+fxt,utyt=Sxtwhere x(t)=[x1(t),x2(t),⋯,xN(t)]T with xi(t) is the temporal coefficient and y(t)=[Y(z1,t),Y(z2,t),⋯,Y(zM,t)]T with Y(zi,t) denotes the measured output on the location zi.The matrices and the nonlinear terms in (3) are given as follows:(4)A=AijN×N,Aij=Aφjz,φizΩ,B=BijN×N,Bij=Bφjz,φizΩ,S=SijM×N,Sij=φjzi,where z1,z2,⋯,zM are the spatial locations.(5)fxt,ut=f1xt,ut,f2xt,ut,⋯,fNxt,utTis the vector of nonlinear terms with(6)fixt,ut=FX,∂X∂z,⋯,U,∂U∂z,⋯,φizΩ.The EEFs-based modeling could be one of the most commonly used DPS modeling methods. It has been applied to system analysis, model reduction, simulations for many complex processes. However, truncation of higher modes with critical dynamics information will greatly influence the model reduction performance in many situations.


## 3. Modified EEFs and Its Applications for Model Reduction

### 3.1. The Modified Principle

It is assumed that{Yl(z)} denote the ensemble data on space location z in the method of snapshots, where l is the time sampling point. EEFs is usually calculated in the situation that NTim, the number of grid points for each snapshots, is much larger than N, the total number of the ensemble {Yl(z)}. It is computationally advantageous to transform the calculation of EEFs as a N×N eigenvalue problem. Thus, the nth EEFs of the eigenvalue problem (Appendix, (A.6)) can be reconstructed from the coefficients cjn and eigenvalue λn via(7)φnz=1λnNc1n⋯cNnY1z⋮YNz.The maximum number of the EEFs in the eigenvalue problem isN, which can be given as follows:(8)φ1z⋮φNz=1NCY1z⋮YNzwhere(9)C=c11λ1⋯cN1λ1⋮⋮c1NλN⋯cNNλN.If the firstk EEFs in (8) are truncated for model reduction of STSs, then we have(10)φ1z⋮φkz=1NC1Y1z⋮YNz,(11)φk+1z⋮φNz=1NC2Y1z⋮YNz,where k≤N and (12)C1=c11λ1⋯cN1λ1⋮⋮c1kλk⋯cNkλk,(13)C2=c1k+1λk+1⋯cNk+1λk+1⋮⋮c1NλN⋯cNNλN.However, truncation of high modes in (8) may have a great influence on performance of the EEF-based reduced model. The solution of obtained nonlinear dynamic systems will have high-sensitive dependence on perturbations, a slight perturbation (e.g., inappropriate truncation of higher modes) would lead to topological changes of the dynamic system. Thus, the modified EEFs are given by adding an extra weight matrix Q into (8)(14)ϕ1z⋮ϕkz=1NQCY1z⋮YNz,where (15)Q=Q11⋯Q1N⋮⋮Qk1⋯QkN.The modified EEFs (14) can be rewritten as follows:(16)ϕ1z⋮ϕkz=1NQ1Q2C1C2Y1z⋮YNz=1NQ1C1Y1z⋮YNz+1NQ2C2Y1z⋮YNzwhere(17)Q1=Q11⋯Q1k⋮⋮Qk1⋯Qkk,Q2=Q1k+1⋯Q1N⋮⋮Qkk+1⋯QkN.Substituting (10) and (11) into (16), the modified EEFs are obtained from the first k EEFs and truncated higher modes as follows:(18)ϕ1z⋮ϕkz=Q1φ1z⋮φkz+Q2φk+1z⋮φNz.It is clear that the weights of initial EEFs in (8) are changed from C to QC. And we can find that each modified EEFs is the linear combination of initial N EEFs as follows:(19)ϕiz=Qi1φ1z+Qi2φ2z+⋯+QiNφNz,where i=1,2,⋯,k. This relationship of modified EEFs and initial EEFs can be rewritten in the following:(20)ϕ1z,ϕ2z,⋯,ϕkz=φ1z,φ2z,⋯,φNzQTwhere k<n and Q denotes the matrix of combined coefficients. The modified EEFs is transformed from initial EEFs using (20), which number is smaller than that of initial EEFs. This indicates that new reduced model based on modified EEFs will have fewer modes compared with the set of initial N EEFs-based models. The coefficients’ matrix is derived from the nonlinear temporal dynamics of DPSs to introduce the information of higher modes to modified EEFs.


### 3.2. The Calculations of Coefficients’ Matrix

LetTP=T1,T2,⋯,Tr1;TiTTi=I,i=1,⋯,r1 be a set of r1 orthogonal N×N matrices, where r1 denotes the number of matrices for excitation/perturbation directions. Let Ms1=c1,c2,⋯,cs1;ci>0,i=1,⋯,s1 be a set of s1 positive constants, where s1 denotes the number of different excitation/perturbation sizes for each direction. Let EP=e1,e2,⋯,ep be p standard unit vectors in Rp, where p denotes the number of inputs. Given a function xt, define the mean by(21)xt=limT→∞1T∫0Txtdt.The definitions of empirical controllability and observability matrices calculated according to the dynamical system (3) are given to compute the coefficients’ matrix. The empirical controllability matrix [29–31] is defined by(22)W^C=∑l=1r1∑m=1s1∑i=1p1r1·s1·cm2∫0∞Φilmtdt,where Φilmt∈RN×N is given by (23)Φilmt=xilmt-x-ilmxilmt-x-ilmT.xilmt and x-ilm are the state and mean state of Nth dynamical system corresponding to the impulsive input ut=cmTleiδt, which is the system trajectory resulting from different excitations. The empirical controllability matrix W^c has the property that its eigenvectors corresponding to nonzero eigenvalues span a subspace which contains the set of states reachable using the chosen initial impulsive inputs.LetTN=T1,T2,⋯,Tr2;TiTTi=I,i=1,⋯,r2 be a set of r2 orthogonal N×N matrices, where r2 denotes the number of matrices for excitation/perturbation directions. Let Ms2=c1,c2,⋯,cs2;ci>0,i=1,⋯,s2 be a set of s2 positive constants, where s2 denotes the number of different excitation/perturbation sizes for each direction. Let EN=e1,e2,⋯,eN be N standard unit vectors in RN. The empirical observability matrix [29–31] is the analogue of the previous one for the output behaviors of the nonlinear system, which is defined by(24)W^O=∑l=1r2∑m=1s21r2·s2·cm2∫0∞TlΨlmtTlTdtwhere Ψlmt∈RN×N is given by (25)Ψijlmt=yilmt-y-ilmTyilmt-y-ilm.yilmt and y-ilm are the output and mean output of Nth dynamical system corresponding to the initial condition x0=cmTlel with ut=0, which is the system trajectory resulting from different perturbations in the initial conditions with the steady inputs.Each of the empirical matrices can be diagonalized by a linear coordinate transformation. This can be done for an empirical matrix by determining the eigenvectors and corresponding eigenvalues of the matrix [29–31]. If the diagonalized ones of the matrices are equal, then the system is said to be in balanced form and the transformation is called a balancing transformation. The balancing-like transformation [30, 32] is used within a Galerkin projection in order to transform the matrices W^C,W^O into the balanced form W-C,W-O. The calculations of transformation can be formulated as finding an invertible state transformation that makes two symmetric positive semidefinite matrices diagonal and equal in the states that are both controllable and observable. The proof that such a transformation exists is given by Zhou and Doyle [33]. The transformation Q- that diagonalizes the matrices and balances the states that are both observable and controllable is shown as follows:(26)Q-W^CQ-T=Q--1TW^OQ--1=diagσ1,σ2,⋯,σNwhere Q- denotes the balanced transformation matrix, σ1≥σ2≥⋯≥σN≥0, and σi’s are the Hankel singular values. The computations for the balanced transformation matrix Q- according to the matrices W^C,W^O are given by Lall [30]. The state-coordinates of the nonlinear systems can be changed and truncated using the Galerkin projection. The row of Q- may be thought of as giving the ‘modes’ of the system associated with Hankel singular values. Thus, the coefficients matrix in (20) is set to be Q=Q-:,1:kT with the MATLAB style colon notation, which is the transposition of first k columns of transformation matrix Q-.


### 3.3. Modified EEFs Based Model Reduction

The spatiotemporal variable of PDE (1) can be expanded onto the modified EEFs ϕiz with corresponding temporal coefficients x-it as follows:(27)Xz,t≈∑i=1kx-itϕiz.For (1), the following is derived:(28)∑i=1kx-˙itϕiz=A∑i=1kx-itϕiz+B∑i=1muithiz+F∑i=1kx-itϕiz,∂∑i=1kx-itϕiz∂z,⋯,∑i=1muithiz,∂∑i=1muithiz∂z,⋯.Using the Galerkin method, the following is then obtained:(29)∫Ω∑i=1kx-˙itϕizϕjzdz=∫ΩA∑i=1kx-itϕiz+B∑i=1muithizϕjzdz+∫ΩF∑i=1kx-itϕiz,∂∑i=1kx-itϕiz∂z,⋯,∑i=1muithiz,∂∑i=1muithiz∂z,⋯ϕjzdz.This will transform (29) into the following ODE system:(30)x-˙t=D-1A-x-t+D-1B-ut+D-1gx-t,utyt=C-x-twhere the D-1 denotes the inverse matrix of D,(31)x-t=x-1t,x-2t,⋯,x-ktT,ut=u1t,u2t,⋯,umtT,gx-t,ut=g1x-t,ut,⋯,gkx-t,utT,gix-t,ut=∫ΩF∑i=1kx-itϕiz,∂∑i=1kx-itϕiz∂z,⋯,∑i=1muithiz,∂∑i=1muithiz∂z,⋯ϕjzdz.LettingQi denotes the i th row of the coefficient matrix Q, then the elements of matrices D,A,B,C can be calculated as follows:(32)Dij=∫Ωϕizϕjzdz=∑l=1NQilQjl∫Ωφizφjzdz=QiQjT,A-ij=∫ΩAϕjzϕizdz=∫ΩA∑l=1NQjlφlz∑l=1NQilφlzdz=QiAQjT,B-ij=∫ΩBhjzϕizdz=∫ΩBhjz∑l=1NQilφlzdz=QiBj,C-ij=ϕjzi=Si1Si2⋯SiNQjT.For simplicity, the equation can be derived as(33)x-˙t=Alx-t+Blut+f-x-t,utyt=Clx-twhere Al=D-1A-, Bl=D-1B-, and Cl=C-,(34)f-x-t,ut=D-1gx-t,ut.


### 3.4. Model Reduction Performances

The effectiveness of the improved EEFs in our previous work was proofed in [28], and model reduction performance of modified EEFs can be given in a similar way. Suppose that there are enough sensors for measurements and YPz,t denotes the predicted spatiotemporal variable for Xz,t. Let Yj=Yz1,tj,Yz2,tj,⋯,YzM,tjT and YPj=YPz1,tj,YPz2,tj,⋯,YPzM,tjT be the measured data of spatiotemporal variable Yz,t and YPz,t at M spatial locations z1,z2,⋯,zM and the sampling times tj, respectively. To evaluate the model reduction performance at any sampling time tj by using modified EEFs and initial EEFs, the root-square error (RSE) as a performance index is introduced for comparisons(35)RSE=∑l=1Mezi,tj2,where ezl,tj=Yzl,tj-YPzl,tj. The effectiveness of model reduction with the modified EEFs is given in the following theorem.Theorem 1.
Given an coefficient matrixQ obtained in Section 3.2, then the RSE based on k modified EEFs is smaller than that based on k initial EEFs at any sampling time tj if EQQT-1QE1E2QTQQT-1E is negative semidefinite, where E is a diagonal matrix with Eii=QiQiT,i=1,2,⋯,k, and (36)E1=Ik000N-k-QTE-1Q,E2=Ik002IN-k-QTE-1Q.Proof.
From POD decomposition for the measured spatiotemporal observation at the sampling timetj, we have (37)Yj=Y-+∑i=1Nyitjφi.
The predicted output based onk initial EEFs at the sampling time tj is as follows:(38)YIEj=Y-+∑i=1kyitjφi.
The predicted output based onk modified EEFs at sampling time tj can be expressed as follows: (39)YEEj=Y-+∑i=1ky-itjϕi.
The RSE withk initial EEFs is defined as(40)GIEtj=∑l=1M∑i=k+1Nyitjφizl2.
The RSE withk modified EEFs is defined as(41)GEEtj=∑l=1M∑i=1Nyitjφizl-∑i=1ky-itjϕizl2.
To prove that0≤GEEtj<GIEtj, only the following inequality has to be proved: (42)GEEtj2<GIEtj2.
Substituting (40) and (41) into (42) yields (43)∑l=1M∑i=1Nyitjφizl-∑i=1ky-itjϕizl2<∑l=1M∑i=k+1Nyitjφizl2.
Then, we have(44)∑l=1M∑i=1Nyitjφizl-∑i=1ky-itjϕizl2-∑i=k+1Nyitjφizl2<0.
And the following inequality can be derived:(45)∑l=1M∑i=1kyitjφizl-∑i=1ky-itjϕizl·∑i=k+1Nyitjφizl+∑i=1Nyitjφizl-∑i=1ky-itjϕizl<0.
Note that(46)y-itj=Yz,t-Y-,ϕiϕi,ϕiwhere ϕi,ϕi=∑j=1NQijφj,∑j=1NQijφj=QiQiT,(47)Yz,t-Y-,ϕi=Yz,t-Y-,∑j=1NQijφj=y1t,⋯,yNtQ1i⋮QNi.
LetEk×k be the diagonal matrix and Eii=QiQiT,i=1,2,⋯,k; then (48)y-1t,⋯,y-kt=y1t,⋯,yNtQTE-1Then (49)y1t,⋯,yNt=y-1t,⋯,y-ktEQQT-1QIn inequality (45), (50)∑i=1kyitjφizl-∑i=1ky-itjϕizl=φ1zl,φ2zl,⋯,φNzlE1y1tjy2tj⋮yNtj,(51)∑i=k+1Nyitjφizl+∑i=1Nyitjφizl-∑i=1ky-itjϕizl=φ1zl,φ2zl,⋯,φNzlE2y1tjy2tj⋮yNtj,where (52)E1=Ik000N-k-QTE-1Q,E2=Ik002IN-k-QTE-1Q.
Then (50) timing (51) yields(53)GEEtj2-GIEtj2=y1tj,y2tj,⋯,yNtjE1E3E2y1tjy2tj⋮yNtj.where(54)E3=∑l=1Mφ1Tzlφ1zl⋯∑l=1Mφ1TzlφNzl⋮⋯⋮∑l=1MφNTzlφ1zl⋯∑l=1MφNTzlφNzl.
Asφ1,φ2,⋯,φN are orthogonal each other, thus(55)E3=10⋯001⋯0⋮⋮⋱⋮0001=INSubstituting (55) and (49) into (53) yields(56)GEEtj2-GIEtj2=y-1tj,y-2tj,⋯,y-ktj·EQQT-1QE1E2QTQQT-1Ey-1tjy-2tj⋮y-ktjIf EQQT-1QE1E2QTQQT-1E is negative semi-definite, then GEEtj<GIEtj because of GIEtj≥0 and GEEtj≥0. This completes the proof and the results will be used in the numerical example. This theorem can only certificate that the modified EEFs is superior to the initial EEFs under certain constrained conditions. The comparisons for the performances of modified EEFs and improved EEFs [28] are given in the numerical example.


## 3.1. The Modified Principle

It is assumed that{Yl(z)} denote the ensemble data on space location z in the method of snapshots, where l is the time sampling point. EEFs is usually calculated in the situation that NTim, the number of grid points for each snapshots, is much larger than N, the total number of the ensemble {Yl(z)}. It is computationally advantageous to transform the calculation of EEFs as a N×N eigenvalue problem. Thus, the nth EEFs of the eigenvalue problem (Appendix, (A.6)) can be reconstructed from the coefficients cjn and eigenvalue λn via(7)φnz=1λnNc1n⋯cNnY1z⋮YNz.The maximum number of the EEFs in the eigenvalue problem isN, which can be given as follows:(8)φ1z⋮φNz=1NCY1z⋮YNzwhere(9)C=c11λ1⋯cN1λ1⋮⋮c1NλN⋯cNNλN.If the firstk EEFs in (8) are truncated for model reduction of STSs, then we have(10)φ1z⋮φkz=1NC1Y1z⋮YNz,(11)φk+1z⋮φNz=1NC2Y1z⋮YNz,where k≤N and (12)C1=c11λ1⋯cN1λ1⋮⋮c1kλk⋯cNkλk,(13)C2=c1k+1λk+1⋯cNk+1λk+1⋮⋮c1NλN⋯cNNλN.However, truncation of high modes in (8) may have a great influence on performance of the EEF-based reduced model. The solution of obtained nonlinear dynamic systems will have high-sensitive dependence on perturbations, a slight perturbation (e.g., inappropriate truncation of higher modes) would lead to topological changes of the dynamic system. Thus, the modified EEFs are given by adding an extra weight matrix Q into (8)(14)ϕ1z⋮ϕkz=1NQCY1z⋮YNz,where (15)Q=Q11⋯Q1N⋮⋮Qk1⋯QkN.The modified EEFs (14) can be rewritten as follows:(16)ϕ1z⋮ϕkz=1NQ1Q2C1C2Y1z⋮YNz=1NQ1C1Y1z⋮YNz+1NQ2C2Y1z⋮YNzwhere(17)Q1=Q11⋯Q1k⋮⋮Qk1⋯Qkk,Q2=Q1k+1⋯Q1N⋮⋮Qkk+1⋯QkN.Substituting (10) and (11) into (16), the modified EEFs are obtained from the first k EEFs and truncated higher modes as follows:(18)ϕ1z⋮ϕkz=Q1φ1z⋮φkz+Q2φk+1z⋮φNz.It is clear that the weights of initial EEFs in (8) are changed from C to QC. And we can find that each modified EEFs is the linear combination of initial N EEFs as follows:(19)ϕiz=Qi1φ1z+Qi2φ2z+⋯+QiNφNz,where i=1,2,⋯,k. This relationship of modified EEFs and initial EEFs can be rewritten in the following:(20)ϕ1z,ϕ2z,⋯,ϕkz=φ1z,φ2z,⋯,φNzQTwhere k<n and Q denotes the matrix of combined coefficients. The modified EEFs is transformed from initial EEFs using (20), which number is smaller than that of initial EEFs. This indicates that new reduced model based on modified EEFs will have fewer modes compared with the set of initial N EEFs-based models. The coefficients’ matrix is derived from the nonlinear temporal dynamics of DPSs to introduce the information of higher modes to modified EEFs.


## 3.2. The Calculations of Coefficients’ Matrix

LetTP=T1,T2,⋯,Tr1;TiTTi=I,i=1,⋯,r1 be a set of r1 orthogonal N×N matrices, where r1 denotes the number of matrices for excitation/perturbation directions. Let Ms1=c1,c2,⋯,cs1;ci>0,i=1,⋯,s1 be a set of s1 positive constants, where s1 denotes the number of different excitation/perturbation sizes for each direction. Let EP=e1,e2,⋯,ep be p standard unit vectors in Rp, where p denotes the number of inputs. Given a function xt, define the mean by(21)xt=limT→∞1T∫0Txtdt.The definitions of empirical controllability and observability matrices calculated according to the dynamical system (3) are given to compute the coefficients’ matrix. The empirical controllability matrix [29–31] is defined by(22)W^C=∑l=1r1∑m=1s1∑i=1p1r1·s1·cm2∫0∞Φilmtdt,where Φilmt∈RN×N is given by (23)Φilmt=xilmt-x-ilmxilmt-x-ilmT.xilmt and x-ilm are the state and mean state of Nth dynamical system corresponding to the impulsive input ut=cmTleiδt, which is the system trajectory resulting from different excitations. The empirical controllability matrix W^c has the property that its eigenvectors corresponding to nonzero eigenvalues span a subspace which contains the set of states reachable using the chosen initial impulsive inputs.LetTN=T1,T2,⋯,Tr2;TiTTi=I,i=1,⋯,r2 be a set of r2 orthogonal N×N matrices, where r2 denotes the number of matrices for excitation/perturbation directions. Let Ms2=c1,c2,⋯,cs2;ci>0,i=1,⋯,s2 be a set of s2 positive constants, where s2 denotes the number of different excitation/perturbation sizes for each direction. Let EN=e1,e2,⋯,eN be N standard unit vectors in RN. The empirical observability matrix [29–31] is the analogue of the previous one for the output behaviors of the nonlinear system, which is defined by(24)W^O=∑l=1r2∑m=1s21r2·s2·cm2∫0∞TlΨlmtTlTdtwhere Ψlmt∈RN×N is given by (25)Ψijlmt=yilmt-y-ilmTyilmt-y-ilm.yilmt and y-ilm are the output and mean output of Nth dynamical system corresponding to the initial condition x0=cmTlel with ut=0, which is the system trajectory resulting from different perturbations in the initial conditions with the steady inputs.Each of the empirical matrices can be diagonalized by a linear coordinate transformation. This can be done for an empirical matrix by determining the eigenvectors and corresponding eigenvalues of the matrix [29–31]. If the diagonalized ones of the matrices are equal, then the system is said to be in balanced form and the transformation is called a balancing transformation. The balancing-like transformation [30, 32] is used within a Galerkin projection in order to transform the matrices W^C,W^O into the balanced form W-C,W-O. The calculations of transformation can be formulated as finding an invertible state transformation that makes two symmetric positive semidefinite matrices diagonal and equal in the states that are both controllable and observable. The proof that such a transformation exists is given by Zhou and Doyle [33]. The transformation Q- that diagonalizes the matrices and balances the states that are both observable and controllable is shown as follows:(26)Q-W^CQ-T=Q--1TW^OQ--1=diagσ1,σ2,⋯,σNwhere Q- denotes the balanced transformation matrix, σ1≥σ2≥⋯≥σN≥0, and σi’s are the Hankel singular values. The computations for the balanced transformation matrix Q- according to the matrices W^C,W^O are given by Lall [30]. The state-coordinates of the nonlinear systems can be changed and truncated using the Galerkin projection. The row of Q- may be thought of as giving the ‘modes’ of the system associated with Hankel singular values. Thus, the coefficients matrix in (20) is set to be Q=Q-:,1:kT with the MATLAB style colon notation, which is the transposition of first k columns of transformation matrix Q-.


## 3.3. Modified EEFs Based Model Reduction

The spatiotemporal variable of PDE (1) can be expanded onto the modified EEFs ϕiz with corresponding temporal coefficients x-it as follows:(27)Xz,t≈∑i=1kx-itϕiz.For (1), the following is derived:(28)∑i=1kx-˙itϕiz=A∑i=1kx-itϕiz+B∑i=1muithiz+F∑i=1kx-itϕiz,∂∑i=1kx-itϕiz∂z,⋯,∑i=1muithiz,∂∑i=1muithiz∂z,⋯.Using the Galerkin method, the following is then obtained:(29)∫Ω∑i=1kx-˙itϕizϕjzdz=∫ΩA∑i=1kx-itϕiz+B∑i=1muithizϕjzdz+∫ΩF∑i=1kx-itϕiz,∂∑i=1kx-itϕiz∂z,⋯,∑i=1muithiz,∂∑i=1muithiz∂z,⋯ϕjzdz.This will transform (29) into the following ODE system:(30)x-˙t=D-1A-x-t+D-1B-ut+D-1gx-t,utyt=C-x-twhere the D-1 denotes the inverse matrix of D,(31)x-t=x-1t,x-2t,⋯,x-ktT,ut=u1t,u2t,⋯,umtT,gx-t,ut=g1x-t,ut,⋯,gkx-t,utT,gix-t,ut=∫ΩF∑i=1kx-itϕiz,∂∑i=1kx-itϕiz∂z,⋯,∑i=1muithiz,∂∑i=1muithiz∂z,⋯ϕjzdz.LettingQi denotes the i th row of the coefficient matrix Q, then the elements of matrices D,A,B,C can be calculated as follows:(32)Dij=∫Ωϕizϕjzdz=∑l=1NQilQjl∫Ωφizφjzdz=QiQjT,A-ij=∫ΩAϕjzϕizdz=∫ΩA∑l=1NQjlφlz∑l=1NQilφlzdz=QiAQjT,B-ij=∫ΩBhjzϕizdz=∫ΩBhjz∑l=1NQilφlzdz=QiBj,C-ij=ϕjzi=Si1Si2⋯SiNQjT.For simplicity, the equation can be derived as(33)x-˙t=Alx-t+Blut+f-x-t,utyt=Clx-twhere Al=D-1A-, Bl=D-1B-, and Cl=C-,(34)f-x-t,ut=D-1gx-t,ut.


## 3.4. Model Reduction Performances

The effectiveness of the improved EEFs in our previous work was proofed in [28], and model reduction performance of modified EEFs can be given in a similar way. Suppose that there are enough sensors for measurements and YPz,t denotes the predicted spatiotemporal variable for Xz,t. Let Yj=Yz1,tj,Yz2,tj,⋯,YzM,tjT and YPj=YPz1,tj,YPz2,tj,⋯,YPzM,tjT be the measured data of spatiotemporal variable Yz,t and YPz,t at M spatial locations z1,z2,⋯,zM and the sampling times tj, respectively. To evaluate the model reduction performance at any sampling time tj by using modified EEFs and initial EEFs, the root-square error (RSE) as a performance index is introduced for comparisons(35)RSE=∑l=1Mezi,tj2,where ezl,tj=Yzl,tj-YPzl,tj. The effectiveness of model reduction with the modified EEFs is given in the following theorem.Theorem 1.
Given an coefficient matrixQ obtained in Section 3.2, then the RSE based on k modified EEFs is smaller than that based on k initial EEFs at any sampling time tj if EQQT-1QE1E2QTQQT-1E is negative semidefinite, where E is a diagonal matrix with Eii=QiQiT,i=1,2,⋯,k, and (36)E1=Ik000N-k-QTE-1Q,E2=Ik002IN-k-QTE-1Q.Proof.
From POD decomposition for the measured spatiotemporal observation at the sampling timetj, we have (37)Yj=Y-+∑i=1Nyitjφi.
The predicted output based onk initial EEFs at the sampling time tj is as follows:(38)YIEj=Y-+∑i=1kyitjφi.
The predicted output based onk modified EEFs at sampling time tj can be expressed as follows: (39)YEEj=Y-+∑i=1ky-itjϕi.
The RSE withk initial EEFs is defined as(40)GIEtj=∑l=1M∑i=k+1Nyitjφizl2.
The RSE withk modified EEFs is defined as(41)GEEtj=∑l=1M∑i=1Nyitjφizl-∑i=1ky-itjϕizl2.
To prove that0≤GEEtj<GIEtj, only the following inequality has to be proved: (42)GEEtj2<GIEtj2.
Substituting (40) and (41) into (42) yields (43)∑l=1M∑i=1Nyitjφizl-∑i=1ky-itjϕizl2<∑l=1M∑i=k+1Nyitjφizl2.
Then, we have(44)∑l=1M∑i=1Nyitjφizl-∑i=1ky-itjϕizl2-∑i=k+1Nyitjφizl2<0.
And the following inequality can be derived:(45)∑l=1M∑i=1kyitjφizl-∑i=1ky-itjϕizl·∑i=k+1Nyitjφizl+∑i=1Nyitjφizl-∑i=1ky-itjϕizl<0.
Note that(46)y-itj=Yz,t-Y-,ϕiϕi,ϕiwhere ϕi,ϕi=∑j=1NQijφj,∑j=1NQijφj=QiQiT,(47)Yz,t-Y-,ϕi=Yz,t-Y-,∑j=1NQijφj=y1t,⋯,yNtQ1i⋮QNi.
LetEk×k be the diagonal matrix and Eii=QiQiT,i=1,2,⋯,k; then (48)y-1t,⋯,y-kt=y1t,⋯,yNtQTE-1Then (49)y1t,⋯,yNt=y-1t,⋯,y-ktEQQT-1QIn inequality (45), (50)∑i=1kyitjφizl-∑i=1ky-itjϕizl=φ1zl,φ2zl,⋯,φNzlE1y1tjy2tj⋮yNtj,(51)∑i=k+1Nyitjφizl+∑i=1Nyitjφizl-∑i=1ky-itjϕizl=φ1zl,φ2zl,⋯,φNzlE2y1tjy2tj⋮yNtj,where (52)E1=Ik000N-k-QTE-1Q,E2=Ik002IN-k-QTE-1Q.
Then (50) timing (51) yields(53)GEEtj2-GIEtj2=y1tj,y2tj,⋯,yNtjE1E3E2y1tjy2tj⋮yNtj.where(54)E3=∑l=1Mφ1Tzlφ1zl⋯∑l=1Mφ1TzlφNzl⋮⋯⋮∑l=1MφNTzlφ1zl⋯∑l=1MφNTzlφNzl.
Asφ1,φ2,⋯,φN are orthogonal each other, thus(55)E3=10⋯001⋯0⋮⋮⋱⋮0001=INSubstituting (55) and (49) into (53) yields(56)GEEtj2-GIEtj2=y-1tj,y-2tj,⋯,y-ktj·EQQT-1QE1E2QTQQT-1Ey-1tjy-2tj⋮y-ktjIf EQQT-1QE1E2QTQQT-1E is negative semi-definite, then GEEtj<GIEtj because of GIEtj≥0 and GEEtj≥0. This completes the proof and the results will be used in the numerical example. This theorem can only certificate that the modified EEFs is superior to the initial EEFs under certain constrained conditions. The comparisons for the performances of modified EEFs and improved EEFs [28] are given in the numerical example.


## 4. Modified EEFs-Based Neural Modelling

Lety={yi(tj)}i=1,j=1N,L be the corresponding temporal coefficients of initial EEFs for the measured output Y={Y(zi,tj)}i=1,j=1M,L. The corresponding temporal coefficients of modified EEFs can be derived as follows:(57)Y,ϕi=Y,∑p=1NQipφp=Qiywhere Qi denotes the i th row of coefficient matrix. The corresponding temporal coefficients y~={y-i(tj)}i=1,j=1k,L can also be derived as follows: (58)y~k×L=Qk×N·yN×L.In the Galerkin method, obtaining an exact analytical description of the ODE systems is difficult and complex because of the nonlinearities in the inner product. Therefore, the neural networks can be used to identify the long-term dynamical behaviors from the inputu(t) and the corresponding temporal coefficients y~ of modified EEFs (59)y^p+1=NNy^p,up,where y^(p)=[y^1(p),y^2(p),⋯,y^k(p)]T,k<N.The advantage of the neural networks is its ability to model complex nonlinear relationships without any assumptions on the nature of these relationships. The most often used neural networks include the radial basis function networks (RBF), backpropagation (BP) neural networks [28], among others. The present study employs a feedforward BP neural network to construct low-dimensional substitute model for the dynamics of DPSs. The prediction output of nonlinear DPSs is obtained by synthesis of temporal predicted outputs and the modified EEFs(60)Ypz,p=Y-+∑i=1ky^ipϕiz.


## 5. A Numerical Example

Suppose thatY(z,t) and Yp(z,t) are the measured and predicted outputs at the M spatial locations z1,z2,⋯,zM and sampling times t1,t2,⋯,tL, respectively. For an easy comparison, the root of mean squared error (RMSE) is set up as the performance index as follows:(61)RMSE=∑i=1M∑j=1Lezi,tj2MLwhere ezi,tj=Yzi,tj-Ypzi,tj.To evaluate the performance of modified EEFs for model reduction, the rescaled Kuramoto-sivashinsky (K-S) equation [34, 35] in one-space dimension is considered. The K-S equations are one of the typical PDEs, which has been derived in 1976 by Kuramoto and Tsuzuki [34] as a model equation for interfacial instabilities in the context of angular phase turbulence for a system of a Reaction-diffusion equation that model the Belouzov-Zabotinskii reaction in three space dimensions, and independently, in 1977, by Sivashinsky [35] to model thermal diffusion instabilities observed in laminar Mame fronts in two space dimensions:(62)∂X∂t+4∂4X∂z4+α∂2X∂z2+12∂X∂z2+∑i=1mhizuit=0where α=84.25; hi(z)=δ[z+3π/4-(i-1)π/2]; m=4;(63)u1t,u2t,u3t,u4t=5cos⁡t4,5sin⁡t4,5cos⁡t2,5sin⁡t2.Equation (62) is subject to periodic boundary condition Xz,t=Xz+2π,t, and the initial condition is set to be cos⁡z. The sampling interval Δt is 0.001s and the simulation time is 0.5s. In this case, forty-one sensors uniformly distributed in the space are used for measurements. A noise dataset of 500 data is collected from (62). This size of data set used for training may be determined by the system complexity and the desired modelling accuracy. More complex system and higher accuracy requirements may need more data from the observation. The solutions of K-S Equation (62) are calculated using the same method in [28], and the performance of this method is compared with the improved EEFs proposed in [28]. After the initial EEFs are derived from the collected spatiotemporal data, a 5-order nonlinear ODE system can be obtained by time/space separation and Galerkin projection. According to the ODE system, coefficient matrix in (20) can be computed by the balanced realization and truncation for its empirical controllability and observability matrices. This computational approach is introduced by Juergen Hahn [29]. Thus, the modified EEFs can be obtained by linear combinations of initial EEFs. A new set of 100 data is collected for testing to compare the performances of two kinds of spatial basis functions (Figure 1). The spatiotemporal output of the K-S equation on testing data can be estimated from the synthesis of the temporal approximate model and spatial basis functions.Figure 1
Measured spatiotemporal output for testing.The first four initial EEFs capture over 99% of energy, the RMSEs of the approximate models based on the first four EEFs, and the first four improved EEFs, and the four modified EEFs are compared in Figure2. The values of RMSE using modified EEFs are smaller than that using the same number of improved EEFs and also much smaller than that using initial EEFs. Because the improved and modified EEFs are derived from linear combinations of initial EEFs (the number is larger than their numbers), the dynamics of the neglectful modes are used to compensate the initial EEFs-based reduced models and then the modeling performance based on improved and modified EEFs is better than initial EEFs-based models. It is worth noticing that when three or four spatial basis functions are used to model reduction for nonlinear DPSs, their model accuracy is close. The reason is that the dynamics information of the neglectful modes becomes much smaller for long-term behaviors and the compensations for the dynamics of the reduced model are restricted.Figure 2
RMSEs based on traditional, improved and modified EEFs.In order to demonstrate the performance of the modified EEFs clearly, model reduction performances by using the first two improved and modified EEFs on the testing data are compared. First, the two improved and modified EEFs are shown in Figures3 and 5, respectively. The predicted spatiotemporal outputs based on two kinds of basis functions are given in Figures 4 and 6, respectively. Compared with the testing data, the predicted distribution errors based on two kinds of basis functions are shown in Figures 7 and 8, respectively. The RMSEs of the approximate model based on the two improved and modified EEFs are 0.0728 and 0.0563, respectively.Figure 3
The first two improved EEFs in [28].Figure 4
The predicted output based on two improved EEFs based modelling.Figure 5
The first two modified EEFs.Figure 6
The predicted output based on two modified EEFs based modelling.Figure 7
Distributed error based on two improved EEFs in [28].Figure 8
Distributed error based on two modified EEFs.Remarks. Nonlinear closure modeling [24, 25] also considered the influences of higher POD modes on the stability and accuracy of the reduced-order model in turbulent, which has shown promising results for Burger’s equation and the Navier-Stokes equations. In its numerical examples, it can be found that the first several POD modes that capture over 99% of energy are used in model reduction will generate less accurate results for the cases. This illustrates the significant role of higher POD modes in the closure model, which improve the accuracy of the reduced order model. Unlike the nonlinear closure modeling approach, the method in this paper actually uses the linear combinations of initial EEFs with more modes to obtain the modified EEFs with fewer modes, and the coefficients are derived from the temporal dynamics to introduce the information of higher modes to the reduced order model. The first couple of modes with 99% of energies are employed to demonstrate the compensation effects of higher modes on the reduced model. The influences by higher modes for dynamics of the reduced model can be found from the variations of RMSEs in Figure 2. However, to further improve the accuracy of the reduced model and optimize the proposed technique, more initial EEFs can be employed to generate the modified EEFs, which can be chosen arbitrarily according to its energies (i.e., all initial modes with non-zero energies are used for linear combinations). This will require higher computational costs to enhance the numerical accuracy. The comparisons with nonlinear closure modeling approach and the utilization of sensitivity analysis [21–23] to enhance the robustness of the reduced model will be carried out in our further research work, which deserve a long-term study.


## 6. Conclusions

In this paper, modified EEFs were proposed for model reduction of the nonlinear STSs, which were derived by adding an extra weight matrix in the method of snapshots. This transforms the derivation of modified EEFs to linear combinations of initial EEFs, while the coefficients matrix was computed according to the nonlinear temporal dynamics of STSs. Thus, the effects of high modes were considered into model reduction with less computational requirements, and the dimension of reduced model is smaller under a given accuracy. The performance of modified EEFs based model reduction is proved theoretically, which are compared with that of the traditional EEFs and the improved EEFs in [28] by a numerical example.



---

*Source: 9435761-2018-11-14.xml*