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.

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 [28]. 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 [2123], 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 [2628] 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:subject to a number of boundary and initial conditions. In (1), denotes the vector of state variable, where is the time variable, is the spatial coordinate, and only one spatial-dimension is considered here. denotes the vector of manipulated spatio-temporal inputs, where and is the th temporal signal with certain spatial distribution . and are two linear operators that involve linear spatial derivatives on the state variable and spatiotemporal input. and denote the partial derivatives of and , respectively. is a nonlinear function containing spatial derivatives for and . A scalar product is defined on spatial domain introduced in the phase space of (1), which is given by .

The spatiotemporal variable is assumed to be expanded onto the set of EEFs (the computational details are given in Appendix) with corresponding temporal coefficients

Suppose that the spatial-temporal output is measured at spatial locations. Substituting the expansion (2) into (1) and taking the inner product with (orthogonality), the integration of equation will give a finite-dimensional nonlinear ordinary differential equation (ODE) system in a general form as follows:where with is the temporal coefficient and with denotes the measured output on the location .

The matrices and the nonlinear terms in (3) are given as follows:where are the spatial locations.is the vector of nonlinear terms with

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 denote the ensemble data on space location in the method of snapshots, where is the time sampling point. EEFs is usually calculated in the situation that , the number of grid points for each snapshots, is much larger than , the total number of the ensemble . It is computationally advantageous to transform the calculation of EEFs as a eigenvalue problem. Thus, the th EEFs of the eigenvalue problem (Appendix, (A.6)) can be reconstructed from the coefficients and eigenvalue via

The maximum number of the EEFs in the eigenvalue problem is , which can be given as follows:where

If the first EEFs in (8) are truncated for model reduction of STSs, then we havewhere and

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 into (8)where

The modified EEFs (14) can be rewritten as follows:where

Substituting (10) and (11) into (16), the modified EEFs are obtained from the first EEFs and truncated higher modes as follows:

It is clear that the weights of initial EEFs in (8) are changed from to . And we can find that each modified EEFs is the linear combination of initial EEFs as follows:where . This relationship of modified EEFs and initial EEFs can be rewritten in the following:where and 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 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

Let be a set of orthogonal matrices, where denotes the number of matrices for excitation/perturbation directions. Let be a set of positive constants, where denotes the number of different excitation/perturbation sizes for each direction. Let be standard unit vectors in , where denotes the number of inputs. Given a function , define the mean by

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 [2931] is defined bywhere is given by and are the state and mean state of th dynamical system corresponding to the impulsive input , which is the system trajectory resulting from different excitations. The empirical controllability matrix 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.

Let be a set of orthogonal matrices, where denotes the number of matrices for excitation/perturbation directions. Let be a set of positive constants, where denotes the number of different excitation/perturbation sizes for each direction. Let be standard unit vectors in . The empirical observability matrix [2931] is the analogue of the previous one for the output behaviors of the nonlinear system, which is defined bywhere is given by and are the output and mean output of th dynamical system corresponding to the initial condition with , 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 [2931]. 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 into the balanced form . 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 that diagonalizes the matrices and balances the states that are both observable and controllable is shown as follows:where denotes the balanced transformation matrix, , and ’s are the Hankel singular values. The computations for the balanced transformation matrix according to the matrices are given by Lall [30]. The state-coordinates of the nonlinear systems can be changed and truncated using the Galerkin projection. The row of 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 with the MATLAB style colon notation, which is the transposition of first columns of transformation matrix .

3.3. Modified EEFs Based Model Reduction

The spatiotemporal variable of PDE (1) can be expanded onto the modified EEFs with corresponding temporal coefficients as follows:For (1), the following is derived:

Using the Galerkin method, the following is then obtained:

This will transform (29) into the following ODE system:where the denotes the inverse matrix of ,

Letting denotes the th row of the coefficient matrix , then the elements of matrices can be calculated as follows:

For simplicity, the equation can be derived as where , , and ,

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 denotes the predicted spatiotemporal variable for . Let and be the measured data of spatiotemporal variable and at spatial locations and the sampling times , respectively. To evaluate the model reduction performance at any sampling time by using modified EEFs and initial EEFs, the root-square error (RSE) as a performance index is introduced for comparisonswhere . The effectiveness of model reduction with the modified EEFs is given in the following theorem.

Theorem 1. Given an coefficient matrix obtained in Section 3.2, then the RSE based on modified EEFs is smaller than that based on initial EEFs at any sampling time if is negative semidefinite, where is a diagonal matrix with , and

Proof. From POD decomposition for the measured spatiotemporal observation at the sampling time , we have The predicted output based on initial EEFs at the sampling time is as follows:The predicted output based on modified EEFs at sampling time can be expressed as follows: The RSE with initial EEFs is defined asThe RSE with modified EEFs is defined asTo prove that , only the following inequality has to be proved: Substituting (40) and (41) into (42) yields Then, we have And the following inequality can be derived:Note that where ,Let be the diagonal matrix and ; then Then In inequality (45), where Then (50) timing (51) yieldswhereAs are orthogonal each other, thusSubstituting (55) and (49) into (53) yieldsIf is negative semi-definite, then because of and . 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

Let be the corresponding temporal coefficients of initial EEFs for the measured output . The corresponding temporal coefficients of modified EEFs can be derived as follows:where denotes the th row of coefficient matrix. The corresponding temporal coefficients can also be derived as follows:

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 input and the corresponding temporal coefficients of modified EEFs where .

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

5. A Numerical Example

Suppose that and are the measured and predicted outputs at the spatial locations and sampling times , respectively. For an easy comparison, the root of mean squared error (RMSE) is set up as the performance index as follows:where .

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:where ; ; ;

Equation (62) is subject to periodic boundary condition , and the initial condition is set to be . The sampling interval 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.

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 Figure 2. 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.

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 Figures 3 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.

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 [2123] 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.

Appendix

The Karhunen–Loève expansion (also known as principal component analysis, principal orthogonal decomposition) is to find an optimal basis from a representative set of process data. Suppose that there is an observation (called snapshots). The problem is how to compute the most characteristic structure among these snapshots .

For simplicity, assume that the observations are uniformly sampled in time and space. Defining the inner product, norm, and ensemble average as , , and .

Motivated by Fourier series, the spatiotemporal variable can be expanded onto an infinite number of orthonormal spatial basis functions with temporal coefficients :

The temporal coefficients can be computed from

denotes the -order approximation

The main procedure of using Karhunen–Loève expansion for time/space separation is computing the most characteristic spatial structure among the spatiotemporal output . This typical structure can be found by minimizing the objective function

The orthogonal constraint is imposed to ensure that the function is unique. The Lagrangian functional corresponding to this constrained optimization problem is

And the solution can be obtained as where is the spatial two-point correlation function. is the ith eigenfunction, and is the corresponding eigenvalue. Given that the covariance matrix is symmetric and positive definite, its eigenvalues are real and its eigenvectors , form an orthogonal set. Since the data are always discrete in space, one must numerically solve the integral (A.6). Discretizing the integral equation gives a matrix eigenvalue problem. Thus, at most eigenfunctions at sampling spatial locations can be obtained.

The maximum number of nonzero eigenvalues is . We arrange the eigenvalues and , in order of the magnitude of the eigenvalues. Each eigenfunction has an energy percentage which depends on the associated eigenvalues of the eigenfunctions:where denotes the sum of the matrix eigenvalues. Assuming that the eigenvalues are sorted in descending order, the eigenfunctions are ordered from most to least energetic. In general, an expansion in terms of only the first few temporal coefficientscan be used to represent the dominant dynamics of DPSs. The common model reduction can be accomplished for the nonlinear DPSs based on the initial EEFs.

Data Availability

The MATLAB programs for the data of the numerical example used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This project is supported by National Natural Science Foundation of China (Grant nos. 51775182, 51775181). Natural Science Foundation of Hunan province (Grant no. 2018JJ3170) and the Outstanding Youth Fund of the Education Department of Hunan Province (Grant no. 16B093) are also gratefully acknowledged.