Optimal Combination of EEFs for the Model Reduction of Nonlinear Partial Differential Equations

Abstract

Proper orthogonal decomposition is a popular approach for determining the principal spatial structures from the measured data. Generally, model reduction using empirical eigenfunctions (EEFs) can generate a relatively low-dimensional model among all linear expansions. However, the neglectful modes representing only a tiny amount of energy will be crucial in the modeling for certain type of nonlinear partial differential equations (PDEs). In this paper, an optimal combination of EEFs is proposed for model reduction of nonlinear partial differential equations (PDEs), obtained by the basis function transformation from the initial EEFs. The transformation matrix is derived from straightforward optimization techniques. The present new EEFs can keep the dynamical information of neglectful modes and generate a lower-dimensional and more precise dynamical system for the PDEs. The numerical example shows its effectiveness and feasibility for model reduction of the nonlinear PDEs.

1. Introduction

Many problems in science and engineering are reduced to a set of nonlinear partial differential equations (PDEs) [1] through a process of mathematical modeling. In this paper, we focus on the modeling of the dynamics for the following kind of nonlinear PDEs: 𝜕𝑢(𝑥,𝑡)()𝑢(),()𝜕𝑡=𝑄𝑢𝑥,𝑡+𝐹𝑥,𝑡𝜕𝑢𝑥,𝑡𝜕𝑥,…,𝑥,𝑡. (1)

Equation (1) will satisfy many types of the boundary conditions. It can use the following description to denote boundary conditions containing the Dirichlet boundary condition and Neumann boundary condition of the nonlinear PDEs (1): 𝐵₁𝑢(𝛼,𝑡)+𝐶₁𝜕𝑢(𝛼,𝑡)𝜕𝑥=𝐷₁,𝐵₂𝑢𝛽,𝑡+𝐶₂𝜕𝑢𝛽,𝑡𝜕𝑥=𝐷₂. (2) The initial conditions are also given: 𝑢(𝑥,0)=𝑢₀(𝑥), (3) where 𝑢(𝑥,𝑡) denotes the spatiotemporal state variable, 𝑥∈Ω=[𝛼,𝛽]⊂𝑅 is the spatial coordinate, and 𝑡∈[0,∞) is the time variable. 𝑄 is a linear operator involving linear spatial derivatives on the state variable. 𝐹(𝑢(𝑥,𝑡),𝜕𝑢(𝑥,𝑡)/𝜕𝑥,…,𝑥,𝑡) is a nonlinear function which contains spatial derivatives of the state variable. 𝐵₁, 𝐶₁, 𝐷₁, 𝐵₂, 𝐶₂, and 𝐷₂ are constants and 𝑢₀(𝑥) is a function of 𝑥.

Generally, two kinds of methods are mainly used to model the dynamics of nonlinear PDEs: conventional discretization approaches and advanced methods based on spatial basis function expansions [2–4]. However, conventional discretization approaches such as finite difference method (FDM) [5] and finite element method (FEM) [6] often lead to high-order dynamical models of nonlinear PDEs. The advanced methods based on spatial basis function expansions such as spectral methods [2, 3] can derive a low-order ODE system to model the dynamics of nonlinear PDEs.

Empirical eigenfunctions (EEFs) [7, 8] which are identified by proper orthogonal decomposition are one of the spatial basis functions used for the time/space separation.

With the help of the traditional identification technique [9], EEFs are widely used for PDE model identification. It is well known that the amount of energy of a system represented by the EEFs is often taken as an indication of the quality of a reduced model using those first several EEFs. However, the linear approximation for the nonlinear system cannot guarantee the assumption that minor components do not contain important information [10]. Modes representing only a tiny amount of energy can also be crucial in the generation of certain types of dynamics. For instance, the EEF-based models can have difficulties in reproducing behavior dominated by irregular transitions between different dynamical states, which have been pointed out in the past studies [11–13]. For the more precise approximation of the PDEs, it is necessary to keep the influence of the neglectful modes representing only a tiny amount of energy on retained dominated modes. Adding the dynamical information by the combination of the initial EEFs will give both qualitatively and quantitatively better results than the common EEFs-based models [8].

In this paper, an optimal combination of EEFs is proposed for model reduction of nonlinear partial differential equations (PDEs), which is obtained by the basis function transformation from the initial EEFs. High-order EEFs are initially obtained by the proper orthogonal decomposition of the measured data. Thus, a new set of basis functions is constructed by the basis function transformation of the initial EEFs. The transformation matrix is derived from straightforward optimization techniques. The present new EEFs can keep the dynamical information of neglectful modes and generate a lower dimensional and more precise dynamical system for the PDEs. The numerical example shows its effectiveness and feasibility for model reduction of the nonlinear PDEs.

2. EEFs for Model Reduction

The proper orthogonal decomposition (POD) [14], also known as the Karhunen-Loeve (K-L) decomposition, principal component analysis (PCA), and empirical orthogonal function analysis (EOF) can facilitate the modal projections of partial differential equations into reduced-order models, via a Galerkin projection [15]. The POD is a mean of extracting spatial structure from a set of time series available on a domain. The most attracting property of the POD is that it can minimize the average squared distance between the original measured signal and its reduced representation. POD allows the identification of a useful set of empirical eigenfunctions and the dimension of the subspace necessary to achieve a satisfactory approximation of nonlinear PDEs.

For simplicity, assume the measured spatial-temporal variable {𝑢(𝑥_𝑖,𝑡_𝑗)}^𝑛,𝐾_{𝑖=1,𝑡=1} (called snapshots) is uniformly sampled in the space and time coordinates. Define the inner product, norm, and ensemble average as (𝑓(𝑥),𝑔(𝑥))=∫_Ω𝑓(𝑥)𝑔(𝑥)𝑑𝑥, ‖𝑓(𝑥)‖=(𝑓(𝑥),𝑓(𝑥))^1/2, and ⟨𝑓(𝑥,𝑡)⟩=1/𝐾∑^𝐾_𝑡=1𝑓(𝑥,𝑡). Motivated by the Fourier series, the spatiotemporal variable 𝑢(𝑥,𝑡) can be expanded onto an infinite number of orthonormal spatial basis functions 𝜙_𝑖(𝑥) with temporal coefficients 𝑢_𝑖(𝑡): 𝑢(𝑥,𝑡)=^∞_𝑖=1𝑢_𝑖(𝑡)𝜙_𝑖(𝑥). (4) The spatial basis functions are orthonormal; the following are derived: 𝜙_𝑖(𝑥),𝜙_𝑗(𝑥)=_Ω𝜙_𝑖(𝑥)𝜙_𝑗(𝑥)𝑑𝑥=0,𝑖≠𝑗,1,𝑖=𝑗. (5) The temporal coefficients can be calculated from 𝑢_𝑖(𝑡)=𝜙_𝑖(𝑥),𝑢(𝑥,𝑡),𝑖=1,…,∞. (6)

If a large enough 𝑀 is given (i.e., 𝑀→∞), 𝑢(𝑥,𝑡) can be approximated by the following 𝑀-order approximation, which is truncated to be a finite dimension: 𝑢(𝑥,𝑡)≈𝑢_𝑀(𝑥,𝑡)=^𝑀_𝑖=1𝑢_𝑖(𝑡)𝜙_𝑖(𝑥), (7) where 𝑢_𝑀(𝑥,𝑡) denotes the 𝑀-order approximation.

The main problem of using proper orthogonal decomposition for time/space separation is to obtain the most characteristic spatial structure {𝜙_𝑖(𝑥)}^𝑁_𝑖=1 from the spatiotemporal snapshots {𝑢(𝑥_𝑖,𝑡_𝑗)}^𝑛,𝐾_{𝑖=1,𝑡=1}. Finding the typical {𝜙_𝑖(𝑥)}^𝑁_𝑖=1 can be performed by using the following eigenvalue problems: _Ω𝑅𝑥,𝜇𝜙_𝑖𝜇𝑑𝜇=𝜆_𝑖𝜙_𝑖(𝑥),𝜙_𝑖,𝜙_𝑖=1, (8) where 𝑅(𝑥,𝜇)=⟨𝑢(𝑥,𝑡)𝑢(𝜇,𝑡)⟩ denotes the spatial two-point correlation function.

Since the data are always discrete in space, one must solve numerically the integral equation (8). Discretizing the integral equation gives an 𝑛×𝑛 matrix eigenvalue problem. Thus, at most eigenfunctions at 𝑛 sampled spatial locations can be obtained. Then one can interpolate the eigenfunctions to locations where the data are not available.

Assume that the maximum number of nonzero eigenvalues is 𝑀≤𝑛. Let 𝜆₁>𝜆₂>⋯>𝜆_𝑀 be the eigenvalues with the corresponding eigenfunctions {𝜙₁(𝑥),𝜙₂(𝑥),…,𝜙_𝑀(𝑥)}. The eigenfunction that corresponds to the first eigenvalue is considered to be the most “energetic.” The total “energy” is defined as being the sum of the eigenvalues. To each eigenfunction, assign an “energy” percentage based on the associated eigenvalue: 𝐸_𝑖=𝜆_𝑖∑^𝑀_𝑗=1𝜆_𝑗. (9)

Usually, the sufficient number of eigenfunctions that capture 99% of the system “energy” is used to determine the value of 𝑁. Experiences show that only the first few basis functions can represent the dominant dynamics of many PDE systems.

As to the model reduction for PDE (1), the spatiotemporal variables 𝑢(𝑥,𝑡) can be approximated by expanding onto the first 𝑁 spatial-dependent EEFs {𝜙₁(𝑥),𝜙₂(𝑥),…,𝜙_𝑁(𝑥)} with the corresponding temporal coefficients 𝑢_𝑖(𝑡): 𝑢(𝑥,𝑡)≈𝑢_𝑁(𝑥,𝑡)=^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜙_𝑖(𝑥). (10) For PDE (1), the following can be derived: 𝑢̇^𝑁_𝑖=1𝑖(𝑡)𝜙_𝑖(𝑥)=𝑄^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜙_𝑖(𝑥)+𝐹^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜙_𝑖(𝑥),^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜕𝜙_𝑖(𝑥).𝜕𝑥,…,𝑥,𝑡 (11) Via a Galerkin projection, the following 𝑁-order dynamical systems are then obtained: 𝑢̇_𝑖(𝑡)=𝜆_𝑖𝑢_𝑖(𝑡)+𝑓_𝑖(𝑢(𝑡),𝑡), (12) where 𝜆_𝑖,𝑖=1,2…,𝑁 denote the eigenvalues, 𝑓_𝑖(𝑢(𝑡),𝑡) denote the nonlinear terms, and 𝑢(𝑡)=[𝑢₁(𝑡),𝑢₂(𝑡),…,𝑢_𝑁(𝑡)]^𝑇. The series of dynamical systems can be rewritten in a general form as follows: 𝑢̇(𝑡)=𝐿𝑢(𝑡)+𝑓(𝑢(𝑡),𝑡), (13) where 𝐿=diag(𝜆₁,𝜆₂,…,𝜆_𝑁);  𝑓(𝑢(𝑡),𝑡)=[𝑓₁(𝑢(𝑡),𝑡),𝑓₂(𝑢(𝑡),𝑡),…,𝑓_𝑁(𝑢(𝑡),𝑡)]^𝑇  and 𝑓_𝑖(𝑢(𝑡),𝑡) = (𝜙_𝑖(𝑥),𝐹(∑^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜙_𝑖(𝑥),∑^𝑁_𝑖=1𝑢_𝑖(𝑡)(𝜕𝜙_𝑖(𝑥)/𝜕𝑥),…,𝑥,𝑡)).

3. Line Combination of EEFs

Considering the influence of neglectful EEFs with the tiny amount energy, a new kind of EEFs can be derived by line combinations of the initial EEFs. As to obtain the lower-order dynamical systems of nonlinear PDEs, a smaller set of new basis functions can be derived. Each new basis function is the line combination of the 𝑀 EEFs from POD: 𝜑_𝑗(𝑥)=𝑎_1𝑗𝜙₁(𝑥)+𝑎_2𝑗𝜙₂(𝑥)+⋯+𝑎_𝑀𝑗𝜙_𝑀(𝑥),𝑗=1,2,…,𝑁, (14) where 𝑁<𝑀. Equation (14) can also be rewritten as follows: 𝜑₁(𝑥),𝜑₂(𝑥),…,𝜑_𝑁(𝑥)=𝜙₁(𝑥),𝜙₂(𝑥),…,𝜙_𝑀(𝑥)𝐴, (15) where the transformation matrix 𝐴=^{𝑎11𝑎12⋯𝑎1𝑁𝑎21𝑎22⋯𝑎2𝑁}_{𝑎⋮⋮⋱⋮𝑀1𝑎𝑀2⋯𝑎𝑀𝑁}.

Using spatial basis functions expansion based on the set of new basis functions (14), the approximation of spatiotemporal variables 𝑢(𝑥,𝑡) can be expanded as follows: 𝑢(𝑥,𝑡)≈^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜑_𝑖(𝑥). (16) Combining (7), (15), and (16), the time coefficients have 𝑢(𝑡)=𝐴⁺𝑢(𝑡), (17) where 𝑢(𝑡)=[𝑢₁(𝑡),𝑢₂(𝑡),…,𝑢_𝑀(𝑡)]^𝑇;  𝑢(𝑡)=[𝑢₁(𝑡),𝑢₂(𝑡),…,𝑢_𝑁(𝑡)]^𝑇, where 𝐴⁺ denotes the generalized inverse of transformation matrix 𝐴.

For PDE (1), the following is derived: 𝑢̇^𝑁_𝑖=1‌_𝑖(𝑡)𝜑_𝑖(𝑥)=𝑄^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜑_𝑖(𝑥)+𝐹^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜑_𝑖(𝑥),^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜕𝜑_𝑖(𝑥).𝜕𝑥,…,𝑥,𝑡 (18)

Using the Galerkin method, a lower-dimensional dynamical system of nonlinear PDE (1) in a general form can be derived as follows: 𝑢̇‌(𝑡)=𝐿𝑢(𝑡)+𝑓(𝑢(𝑡),𝑡), (19) where 𝑢(𝑡)=[𝑢₁(𝑡),𝑢₂(𝑡),…,𝑢_𝑁(𝑡)]^𝑇; 𝐿=𝐿_𝑖𝑗,𝑖,𝑗=1,2,…,𝑁;𝐿_𝑖𝑗=𝜑_𝑖,𝐿𝜑_𝑗;𝑓(𝑢(𝑡),𝑡)=𝑓₁(𝑢(𝑡),𝑡),𝑓₂(𝑢(𝑡),𝑡),…,𝑓_𝑁(𝑢(𝑡),𝑡)^𝑇;𝑓_𝑖(𝑢(𝑡),𝑡)=𝜑_𝑖(𝑥),𝐹^𝑁_𝑖=1𝑢_𝑖(𝑡)𝜑_𝑖(𝑥),_𝑁𝑖=1𝑢_𝑖(𝑡)𝜕𝜑_𝑖(𝑥).𝜕𝑥,…,𝑥,𝑡 (20)

4. Optimization for the Transformation Matrix

An error function that measures energy error between the ODE system (13) and the lower-dimensional dynamical system (19) is introduced by Error=^∞₀^𝑀_𝑖=1𝑢_𝑖(𝑡)²−^𝑁_𝑖=1𝑢_𝑖(𝑡)²²𝑑𝑡. (21) Combining the vector expressions of 𝑢(𝑡),  𝑢(𝑡) with (17), the error (21) can be strictly derived as follows: Error=^∞₀^𝑀_𝑖=1𝑢_𝑖(𝑡)²−^𝑁_𝑖=1𝑢_𝑖(𝑡)²²=𝑑𝑡_∞0𝑢(𝑡)^𝑇𝑢(𝑡)−𝑢(𝑡)^𝑇𝑢(𝑡)²=𝑑𝑡_∞0𝑢(𝑡)^𝑇𝑢(𝑡)−𝑢(𝑡)^𝑇𝐴^+𝑇𝐴⁺𝑢(𝑡)²=𝑑𝑡_∞0𝑢(𝑡)^𝑇𝐼−𝐴^+𝑇𝐴⁺𝑢(𝑡)²𝑑𝑡. (22)

The combination matrix 𝐴 is determined by minimizing the error function (22). To evaluate the error function, a temporal interval [0,𝑇_max] is considered for (22) and approximated by a finite sum, where the max integral time 𝑇_max is selected as free parameter. The temporal integral interval [0,𝑇_max] generally contains the processes from initial state to steady state of the nonlinear PDEs. The restricted condition for the minimizing problem (22) is that the combination matrix should be column orthogonal. Thus, the error functions for combination matrix can be approximated as follows.

Since Error(𝐴)≈^𝑇max₀𝑢(𝑡)^𝑇𝐼−𝐴^+𝑇𝐴⁺𝑢(𝑡)²≈𝑑𝑡_𝑖𝑡_𝑖−𝑡_𝑖−1𝑢𝑡_𝑖^𝑇𝐼−𝐴^+𝑇𝐴⁺𝑢𝑡_𝑖², (23) we consider approximate error 𝐸(𝐴)=_𝑖𝑡_𝑖−𝑡_𝑖−1𝑢𝑡_𝑖^𝑇𝐼−𝐴^+𝑇𝐴⁺𝑢𝑡_𝑖²,subjectto𝐴(∶,𝑖)^𝑇𝐴∶,𝑗=𝑘_𝑖𝛿_𝑖𝑗. (24)

The optimization of 𝐸(𝐴) (24) mainly contains several steps. Firstly, the fourth-order Runge-Kutta method is used to calculate the ODE system (13); thus 𝑡_𝑖, 𝑢(𝑡_𝑖) of Error (24) are known.

Secondly, the particle swarm optimization (PSO) algorithm [16] is used to optimize the restricted error functions (24). PSO is an efficient, robust, and simple optimization algorithm for solving many optimization problems. In nature, PSO is a stochastic optimization approach which maintains a swarm of candidate solutions, referred to as particles. The method is inspired by the movement of particles and their interactions with their neighbors in the group. The principle and details of PSO algorithm are given in [17].

For the optimization of (24) by PSO directly, the obtained combination matrix 𝐴 is not orthogonal. To keep the minimization of the error function (24), a satisfied combination matrix 𝐴 can be derived by singular value decomposition for 𝐴. Let the singular value decomposition of 𝐴 be 𝐴=𝑈𝑆𝑉^𝑇, where 𝑈,  𝑆 are orthogonal and diagonal matrices, respectively. Thus combination matrix 𝐴=𝑈𝑆 is orthogonal and it can be derived that (𝐴⁺)^𝑇𝐴⁺=𝑈𝑆(𝑆^𝑇𝑆)⁻¹(𝑆^𝑇𝑆)⁻¹𝑆^𝑇𝑈^𝑇=(𝐴⁺)^𝑇𝐴⁺, which means that the minimum of the error function (24) is retained.

5. Numerical Example

In this section, the optimal combination of EEFs is applied for model reduction of a simple PDE model. Note that this example serves as an illustration of the proposed methodology rather than a demonstration of the reduction of a very large and complex system.

In order to evaluate the performance of optimal combination of EEFs for model reduction of nonlinear PDEs, the Burgers equation [18, 19] is studied. The Burgers equation is a one-dimensional spatial model of a variety of three-dimensional physical phenomena, greatly simplifying the problems while retaining many of the complex behavior characteristics. This equation contains nonlinear convection and diffusion terms and retains many of the interesting features of the Navier-Stokes equation. The governing equation may be written as follows: 𝜕𝑢(𝑥,𝑡)()𝜕𝑡+𝑢𝜕𝑢𝑥,𝑡𝜕𝜕𝑥=𝜀²𝑢(𝑥,𝑡)𝜕𝑥², (25) where the viscosity 𝜀 is 0.2 and (25) subjects to the following boundary conditions: 𝑢(0,𝑡)=𝑢(2𝜋,𝑡)=0, (26) and an initial condition: 𝑢(𝑥,0)=sin(𝑥)+sin(2𝑥)+sin(3𝑥). (27)

Thirty-two sensors uniformly distributed in the space are used for measurement. A noise-free dataset of 300 data is collected from (25). The sampling interval Δ𝑡 is 0.01 s and the simulation time is 3 s.

5.1. Comparisons with the EEFs

The EEFs are obtained from the singular value decomposition for the sampling data. The energy values of the dynamical systems based on the EEFs are shown in Table 1.

Table 1 

Energy value of the dynamical systems based on EEFs.

Using spatial basis functions expansion based on the new basis functions and the Galerkin methods, low-dimensional dynamical systems can be derived. Meanwhile, the energy functions of the obtained dynamical systems are calculated. The results in Table 2 show that the low-dimensional dynamical system based on the first several optimal combinations of EEFs can capture the dynamics of the nonlinear PDE (25) perfectly.

Table 2 

Energy value of the dynamical systems based on new basis functions.

5.2. Low-Dimensional Approximation Based on 3 Optimal Combined EEFs

Using the line combinations of the EEFs and optimization for the energy error functions of the line combination coefficients, the optimal combination of EEFs for model reduction of nonlinear PDEs is obtained, where the first 3 initial EEFs and 3 optimal combined EEFs are shown in Figures 1 and 2.

Figure 1 

The first 3 initial EEFs.

Figure 2 

The first 3 optimal combined EEFs.

A new set of data is collected for testing to demonstrate the spatial-temporal performance of the two kinds of EEFs. For practical industrial control, the low-dimensional model based on optimal combined EEFs is very useful because it cannot only capture the dominant dynamics of the nonlinear PDEs but simplify the control design as well. The spatiotemporal errors of the predicted output of two models combined with two kinds of EEFs are shown in Figures 3 and 4, respectively. As given in Figures 3 and 4, the predicted distributed error of the approximate model based on the optimal combination of EEFs is shown to be significantly smaller than that of the same order approximate model based on initial EEFs.

Figure 3 

Predicted error based on 3 initial EEFs.

Figure 4 

Predicted error based on 3 optimal combined EEFs.

6. Conclusions

In this note, an optimal combination of EEFs is proposed for model reduction of nonlinear partial differential equations (PDEs), which is obtained by the basis function transformation from the initial EEFs. The transformation matrix is determined by the optimization for the energy error functions of the dynamical systems; the present new EEFs can keep the dynamical information of neglectful modes and generate a lower-dimensional and more precise dynamical system for the PDEs. The numerical example shows its effectiveness and feasibility for model reduction of the nonlinear PDEs.

Acknowledgment

Financial support from National Natural Science Foundation of China (11271376) is gratefully acknowledged.