Access the full text.
Sign up today, get DeepDyve free for 14 days.
Loading next page...
/lp/springer-journals/model-order-reduction-for-bayesian-approach-to-inverse-problems-rZ7ujcJFtm
- Publisher
- Springer Journals
- Copyright
- Copyright © 2014 by Nguyen et al.; licensee Springer.
- Subject
- Engineering; Theoretical and Applied Mechanics; Computational Science and Engineering; Classical Continuum Physics; Mathematical Applications in the Physical Sciences
- eISSN
- 2196-1166
- DOI
- 10.1186/2196-1166-1-2
- Publisher site
- See Article on Publisher Site
Abstract
mpennh@nus.edu.sg Singapore-MIT Alliance, National This work presents an approach to solve inverse problems in the application of water University of Singapore, Singapore quality management in reservoir systems. One such application is contaminant 117576, Singapore Full list of author information is cleanup, which is challenging because tasks such as inferring the contaminant location available at the end of the article and its distribution require large computational efforts and data storage requirements. In addition, real systems contain uncertain parameters such as wind velocity; these uncertainties must be accounted for in the inference problem. The approach developed here uses the combination of a reduced-order model and a Bayesian inference formulation to rapidly determine contaminant locations given sparse measurements of contaminant concentration. The system is modelled by the coupled Navier-Stokes equations and convection-diffusion transport equations. The Galerkin finite element method provides an approximate numerical solution-the ’full model’, which cannot be solved in real-time. The proper orthogonal decomposition and Galerkin projection technique are applied to obtain a reduced-order model that approximates the full model. The Bayesian formulation of the inverse problem is solved using a Markov chain Monte Carlo method for a variety of source locations in the domain. Numerical results show that applying the reduced-order model to the source inversion problem yields a speed-up in computational time by a factor of approximately 32 with acceptable accuracy in comparison with the full model. Application of the inference strategy shows the potential effectiveness of this computational modeling approach for managing water quality. Keywords: Bayesian; Convection-diffusion equation; Navier-Stokes equations; Markov chain Monte Carlo; Inverse problem; Proper orthogonal decomposition; Reduced-order model Background Hydrodynamic processes such as contaminant transport in lakes and reservoirs have a direct impact on water quality. The contaminants will appear, spread out, and decrease in concentration, etc. because of some processes such as convection, diffusion, time rate release of contaminants, and distance of travel. To simulate such processes, a coupled sys- tem of partial differential equations (PDEs) including the Navier-Stokes equations (NSEs) and contaminant transport equations needs to be solved. A better understanding of these processes is important in managing water resources effectively. The direct or forward problems compute the distribution of contaminant directly from given input information such as contaminant location, contaminant properties, fluid flow © 2014 Nguyen et al.; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 2 of 17 http://www.apjcen.com/content/1/1/2 properties, boundary conditions, initial conditions, etc. On the contrary, the inverse prob- lems infer the unknown physical parameters, boundary conditions, initial conditions, or geometry given a set of measured data. These known data can be obtained experimen- tal or computational. In realistic applications, data are not perfect because of error due to sensor noise. In addition, the model may contain some uncertain parameters such as wind velocity. Thus, inverse problems are generally stochastic problems. A Bayesian inference approach to inverse problems is one way to deal with the stochastic problem. This approach takes into account both the information of the model parameters with uncertainty and the inaccuracy of data in terms of a probability distri- bution [1,2]. The Bayesian approach provides a general framework for the formulation of wide variety of problems such as climate modeling [3], contaminant transport modeling [4-6], and heat transfer [7]. Under the Bayesian framework, simulated solutions need to be evaluated repeatedly over different samples of the input parameters. There are available sampling strategies associated with Bayesian computation such as Markov chain Monte Carlo (MCMC) methods [8-11]. However, if we use traditional PDE discretization meth- ods, such as finite element or finite volume methods, the resulting numerical models describing the system will be very large and expensive to solve. This paper presents an efficient computational approach to solve the statistical inverse problem. The approach uses the combination of a reduced-order model and a Bayesian inference formulation. The Galerkin finite element method [12] provides an approximate numerical solution - the ‘full model’ or ‘forward model’. The MCMC method is applied to solve the Bayesian inverse problems. In particular, we are interested in inferring an arbi- trary source location in a time-dependent convection-diffusion transport equation, given a velocity field and a set of measured data of the concentration field at sparse spatial and temporal locations. We obtain the velocity field by solving the Navier-Stokes equations. Solution of large-scale inverse problems can be accelerated by applying model order reduction [13-15]. The idea is to project the large-scale governing equations onto the sub- space spanned by a reduced-space basis, yielding a low-order dynamical system. Proper orthogonal decomposition (POD) is the most popular method to find the reduced basis for a given set of data. The snapshot POD method, which was proposed by Sirovich [16], provide an efficient way for determining the POD basis vectors. Snapshots are solutions of a numerical simulation at selected times or choices of the parameters. The choice of snapshots should ensure that the resulting POD basis captures the most important characteristics of the system. Since the convection term in contaminant transport equations is a velocity-dependent term, we cannot apply the standard POD reduction framework. The evaluation of the velocity-dependent convection term in the reduced system still depends on the full finite element dimension and has the same complexity as the full-order system. Thus, a sys- tem with velocity-dependent convection term needs an additional treatment to obtain an efficient reduced-order model. This paper is outlined as follows. In the section ‘Mathematical model and numerical methods’, we introduce the mathematical model, the finite element approximate tech- nique, and the Bayesian inference formulation and its solution using an MCMC method. In section ‘Model order reduction’, we describe the model order reduction technique and the treatment of the velocity-dependent convection term. In the section ‘Numerical example’, we use the numerical example to demonstrate the solution of the statistical 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 3 of 17 http://www.apjcen.com/content/1/1/2 inverse problem and the reduced-order model performance. We provide some concluding remarks in the final section. Methods In this section, the mathematical model, finite element approximation technique, Bayesian inference formulation, and MCMC method are briefly introduced. Transport model Consider the fluid flow through a physical domain D ⊂ R , d = 1, 2, 3 with boundary . The contaminant concentration c :≡ c(x, t; θ) satisfies the dimensionless parabolic PDEs, boundary conditions, and initial conditions as follows: ∂c 1 + u ·∇c − c = f in D×[t , t ], (1) 0 f ∂t Pe c = g on ×[t , t ], (2) D 0 f ∂c = 0on ×[t , t ], (3) N 0 f ∂n c(x, t ) = c (x) in D,(4) 0 0 where x ∈ D denotes the spatial coordinates, t ∈ [t , t ] denotes time, and c is the given 0 f 0 initial condition. The inlet boundary is subjected to a Dirichlet condition, while the remainder of the boundary = \ satisfies Neumann conditions. The external N D source f (x, t; θ) in Equation 1 is used here as a input parameter, where θ ∈ D is the source location. In this study, we suppose that the contaminant source drives the trans- uL port process, given a velocity field u ∈ R . Thus, the Péclet number Pe = ,where κ is the diffusivity and L is the characteristic length. The velocity field is a function of x and t, u(x, t), which can be obtained by solving the Navier-Stokes equations, as follows: ∂u 1 + u ·∇ u − ∇ u +∇p = f in D× [t , t ], (5) c 0 f ∂t Re ∇· u =0in D× [t , t ], (6) 0 f u = u on × [t , t ], (7) D 0 f u(x, t ) = u (x) in D,(8) 0 0 where Re is the Reynolds number, p is the pressure field, and f is the body force. So, for any given f and u, one can solve for the solution c(x, t; θ). Finite element approximations The finite element method [12] associated with the stabilized second-order fractional- step method is applied to solve the incompressible Navier-Stokes equations (5 to 8). The detailed discussion and derivation of this method is beyond the scope of this paper. For more details, please refer to [17,18]. We also use the finite element method for spatial discretization of the contaminant transport equations in Equations 1 to 4. This leads to the semi-discrete equations Mc ˙ + C(u) + K c = F(t; θ),(9) c(t ) = c , (10) 0 0 with the output of contaminant solution at sensor locations in the domain y = Bc. (11) 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 4 of 17 http://www.apjcen.com/content/1/1/2 Here, c(t) ∈ R is the discretized approximation of c(x, t) and contains N state unknowns where N is the number of grid points in the spatial discretization. c ˙ is the derivative of c with respect to time. The vector y ∈ R contains the N outputs of the system related to N ×N N ×N the state through the matrix B ∈ R , N N. The matrices M, K, C(u) ∈ R and the vector F(t; θ) ∈ R are defined in the finite element space as follows: M = φ (x)φ (x) dx, (12) i,j i j K = ∇φ (x) ·∇φ (x) dx, (13) i,j i j Pe C(u) = u ·∇φ (x)φ (x) dx, (14) i j i,j F(t; θ) = f (t; θ)φ (x) dx. (15) Here φ with i = 1, ··· , N is the finite element basis and θ ∈ R is thecoordinatevector of the source location. Bayesian inference to inverse problems The Bayesian formulation of an inference problem is derived from Bayes’ theorem [1]: p(y|θ)p(θ ) p(θ |y) = , (16) p(y) N d where y ∈ R is the observed data and θ ∈ R ⊆ D is the vector of input parameters. p(θ ) denotes the prior distribution of the parameters, p(y|θ) is the likelihood function, and p(θ |y) is the posterior distribution of parameters given the observed data. p(y) is the prior predictive distribution of y and may be set to an unknown constant that is not needed in the MCMC solution method. In this case, we write Bayes’s formula as p(θ |y) ∝ p(y|θ)p(θ ). (17) A simple model for the likelihood function can be obtained from the following relation- ship, y = G(θ ) + η, (18) where η is the noise which we assume to be white Gaussian noise η ∼ N (0, σ I). The input-output relative in Equations 1 to 4 is denoted as the forward model G(θ ). Subsequently, the likelihood can be written as 1 1 p(y|θ) = √ exp − y − G(θ ) . (19) 2σ σ 2π There are many approaches to define the prior information, such as Gaussian Markov random fields or uniform distribution models, etc. In this work, we assume that our only prior information on the source location is given by the bounds on the domain. Thus, using the principle of maximum entropy [19], we take our prior to be a uniform distribution. Therefore, Equation 17 becomes p(θ |y) ∝ exp − (y − G(θ ) ) (y − G(θ ) ) . (20) i i i i 2σ i=1 Here, K is the number of time steps over which we collect the output data. 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 5 of 17 http://www.apjcen.com/content/1/1/2 Markov chain Monte Carlo for posterior sampling MCMC simulation provides a sampling strategy from a proposal distribution to the target distribution using the Markov chain mechanism. Hastings introduced the Metropolis- Hastings (MH) algorithm [10], which is a generalized version of the Metropolis algorithm [8,9], in which we can apply both symmetric and asymmetric proposal distributions. In this work, the MH algorithm is used to solve the Bayesian inverse problem. The MH algorithm is summarized as follows: 1. Initialize the chain θ and set n = 0 2. Repeat n = n + 1 ∗ ∗ Generate a proposal point θ ∼ q(θ |θ) Generate U from a uniform U(0, 1) distribution n+1 Update the state to θ as θ , if β< U n+1 θ = θ , otherwise 3. Until n = N → stop. mcmc Here, β is the acceptance-rejection ratio, given by ∗ n−1 ∗ π (θ ) q θ |θ β = min 1, . (21) n−1 ∗ n−1 π θ q (θ ) |θ N is the total number of samples, π(θ) is the target distribution, q(θ |θ) is a proposal mcmc distribution, and U is a random number. The discretized model in the form of ordinary differential equations (ODEs) (Equations 9 to 11) may have very large dimensions and be expensive to solve. The MCMC method requires evaluating repeatedly the solution of the forward model (many thousands or even millions of times). Hence, these simulations can be a computation- ally expensive undertaking. In such situations, the reduced-order model is needed to approximate the large-scale model, which allows efficient simulations. Model order reduction This section presents the model order reduction framework. This includes the reduction via projection, the proper orthogonal decomposition method, the additional treatment for the velocity-dependent convection term, and the error estimation. Reduction via projection A reduced-order system approximating the ODEs (9 to 11) can be obtained by approx- imating the full state vector c as a linear combination of m basis vectors as follows: c ≈ Vc , (22) r 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 6 of 17 http://www.apjcen.com/content/1/1/2 m N ×m where c ∈ R is the reduced order state and V = [v v ··· v ] ∈ R is an orthonor- r 1 2 m mal basis, i.e., V V = I. The low-order model can be derived by projecting the governing system (9 to 11) onto the reduced space formed by the column span of basis V. This yields M c ˙ + C (u) + K c = F (t; θ), (23) r r r r r r c (t ) = c , (24) r 0 0r y = B c , (25) r r r where M = V MV , (26) K = V KV , (27) C (u) = V C(u)V , (28) F (t; θ) = V F(t; θ), (29) B = BV , (30) c = V c . (31) 0r 0 The model reduction task is to find a suitable basis V so that m N and the reduced- order model yields accurate results. This study will consider POD as the method to compute the basis. Proper orthogonal decomposition POD provides a method to compute the reduced-order basis V and construct the low- order system. Here, we briefly describe the general POD method (more details may be found in [16]). N ×Q Supposed that we have the snapshot matrix X ∈ R ,where Q is the number of snapshot solutions, which are built from the instantaneous solution c (t ) of the forward model, corresponding to s = 1, ··· , S random input source locations, picked up at k = 1, ··· , T different time steps (so that Q = ST). The POD basis vectors are the m left singular vectors of X corresponding to the largest singular values (m ≤ Q).Let σ , i = 1, 2, ··· , Q be the singular values of X in decreasing order. We choose the number of POD basis vectors to retain in the reduced-order model as m ≤ Q so that m Q 2 2 σ / σ ≥ , (32) i j i=1 j=1 where (%) is the required amount of energy, typically taken to be 99% or higher. After obtaining the POD basis vectors for the contaminant, we have defined our reduced-order system. However, as mentioned earlier, the dimension of the system of Equations 23 to 25 is reduced only in state (concentration). The reduced convection matrix (C (u)) in Equation 28 still depends on the full dimension of the velocity field as in Equation 14. This means that we have to re-compute the full convection matrix (C(u)) at each time-step before projecting it onto the reduced-space basis to obtain the reduced convection matrix (C (u)). Hence, we need an additional treatment to avoid this computational cost. 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 7 of 17 http://www.apjcen.com/content/1/1/2 Linear expansion techniques for velocity-dependent term We proceed by representing the velocity field as a linear combination of POD velocity basis vectors. The velocity-dependent convection matrix can then be expanded based on this velocity expansion. In particular, let u(x, t) be a given velocity field and u(x, t ) k=1 be a set of ‘snapshots’ of velocities which are obtained from Equations 5 to 8 and taken over T time steps. The velocity field is decomposed as u(x, t) = u (x) + u (x, t), (33) 1 k where u (x) = u(x, t ) is the mean velocity field and u (x, t) is the fluctuating T k=1 velocity field. The fluctuating velocity field is then represented by u (x, t) = α (t) (x), (34) k k k=1 where (x) is the kth POD velocity basis and α (t) is the corresponding time dependent k k amplitude. We now consider the expansion of the velocity field as u(x, t) = u (x) + α (t) (x), (35) m k k k=1 where N T is thenumberofPOD basisvectors usetorepresent thevelocity. Figure 1 The domain of interest. (a) Physical domain. (b) Computational domain. 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 8 of 17 http://www.apjcen.com/content/1/1/2 Galerkin projection of the time-dependent velocity term in Equation 1 has the form as u ·∇c, V = u + α ·∇c, V i m k k i k=1 = [u ·∇c, V ] + α ·∇c, V . (36) m i k k i k=1 ThediscreteformofEquation36thenbecomes C (u) = V , u ·∇V + α V , ·∇V r i m j i j k k k=1 = C u + α C . (37) rm m k rk k k=1 Here, C u (x) and C (x) , k = 1, ··· , N are the reduced-order forms of the full- rm m rk k u order convection matrices C u (x) and C (x) , respectively. They are computed only m k once in the offline step. Thus, the first case of the reduced-order model (with Equation 28) is only reduced in state (concentration), but the second case (with Equation 37) is now reduced in both state (concentration) and parameter (velocity). Error estimation In order to estimate the efficiency of the reduced model relative to the full model, we use the relative errors of state solutions and relative errors of outputs. These errors are defined as follows: t = 4.00 t = 8.00 ab 1.5 1.5 0.5 0.5 2 2 Z axis 0 Z axis X axis 1 X axis t = 32.00 cd t = 16.00 0.8 0.3 0.6 0.2 0.4 0.1 0.2 Z axis Z axis 0 1 X axis X axis Figure 2 Contaminant field of forward model at specific times. Figures are taken at (a) t = 4.0s, (b) t = 8.0s, (c) t = 16.0s and (d) t = 32.0s. 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 9 of 17 http://www.apjcen.com/content/1/1/2 Eigenvalues of velocity field 0 100 200 300 400 500 Relative error of velocity field 0 20 40 60 Number of POD modes Figure 3 Properties of POD velocity basis vectors: (a) eigenvalues and (b) relative error analysis. k k c t ; θ − Vc t ; θ r L (D) ε t , θ = , (38) c t ; θ L (D) k k y t ; θ − y t ; θ r L (D) ε t , θ = . (39) y t ; θ L (D) k k k Here, c t ; θ , c t ; θ ,1 ≤ k ≤ T are the full and reduced solutions, while y t ; θ and y (t ; θ),1 ≤ k ≤ T are the full and reduced outputs of interest, respectively. Results and discussion We consider a numerical example based on a 2D mathematical model. The velocity field in the reservoir is obtained by solving a 2D laterally averaged Navier-Stokes model. The Bayesian formulation of the inverse problem is then solved to determine an uncertain con- taminant source location. We compare the effectiveness of solving the inverse problem using both the full model and reduced model techniques. Model setup The physical domain is illustrated in Figure 1a, which represents a simplified model of a 2D laterally averaged reservoir system. The reservoir system includes a main reservoir σ / Σ σ ||u u || /||u|| i i r 2 2 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 10 of 17 http://www.apjcen.com/content/1/1/2 section and the river connections. We assume that the contaminant is present within the main reservoir section and that the contaminant transport processes are mainly affected by the inflow and wind velocity. We assume that other factors (such as heat exchange, etc.) have little influence on the processes. The spatial domain is discretized with a finite element mesh, as shown in Figure 1b. The total number of grid points is N = 1, 377 and the total number of elements is N = 2, 560. The simulation is run from initial time t = 0 E 0 to final time t = 40, with time-step size t = 0.08. Thus, the number of time steps is T = 500. A time-dependent velocity field is obtained from the 2D laterally averaged system, T T which is given in Equations 5 to 8, where the body force f = [ f ; f ] = [0; g] with g c cx cz the gravitational acceleration. The boundary conditions are set up as follows (refer to Figure 1a for the boundary definitions): (u, w) = (1, 0) on , (40) (u, w) = (−16 ∗ (2.0 + z) ∗ (1.5 + z),0) on , (41) (u, w) = (−16 ∗ (1.0 + z) ∗ (0.5 + z),0) on , (42) (u, w) = (0.03V ,0) on , (43) a 11 p = 0on . (44) The velocity on remaining boundaries are set to zero. Here, V is the wind speed at 10 m above the water surface. In this example, we assumed that V = 2m/s forthe entire simulation time. The Reynolds number is set to Re = 5.0e ,and amixinglength turbulence model is used [20]. The full model system is given in Equations 9 to 10. The Crank-Nicolson method [21] is used to discretize the system in time. −0.5 −0.5 −1 −1 −1.5 −1.5 −2 −2 −1 0 1 2 −1 0 1 2 X axis X axis c d −0.5 −0.5 −1 −1 −1.5 −1.5 −2 −2 −1 0 1 2 −1 −0.5 0 0.5 1 1.5 2 2.5 X axis X axis Figure 4 The first four POD basis vectors of the velocity field. (a) First POD velocity basis vector. (b) Second POD velocity basis vector. (c) Third POD velocity basis vector. (d) Fourth POD velocity basis vector. Z axis Z axis Z axis Z axis 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 11 of 17 http://www.apjcen.com/content/1/1/2 U profile at x = 0.5 V profile at x = 0.5 ab 0.5 0.5 1.5 1.5 Ur(mor) Vr(mor) U(full) V(full) 2 2 0.5 0 0.5 1 1.5 0.2 0.1 0 0.1 0.2 0.3 x x cd U profile at x = 1 V profile at x = 1 0.5 0.5 1 1 1.5 1.5 Vr(mor) Ur(mor) U(full) V(full) 2 2 0.5 0 0.5 1 0.3 0.2 0.1 0 0.1 Figure 5 Velocity profile taken at t = t at specific locations. Figures are captured (a) U profile at x = 0.5, (b) V profile at x = 0.5, (c) U profile at x = 1and (d) V profile at x = 1. Here U(full) and V(full) are the original velocity while Ur(mor) and Vr(mor) are the approximate velocity using a POD expansion with N = 12 modes. The setup for the contaminant transport part of the problem is as follows. We use a contaminant source that is defined as the superposition of Gaussian sources, each one active on the time interval t ∈ [t , t ]and centered at θ ∈ D,withstrength h and 0 off 0k k k width σ : sk h |θ − x| k k f (x, t; θ) = exp − δ(t − t ). (45) 0k 2 2 2πσ 2σ sk sk k=1 Here, to simplify the problem, we choose the number of sources to be n = 1, located at θ = (x , z ),withstrength h = 0.2 and width σ = 0.05. The active time of the 1 c c 1 s1 source is t ∈ [0, t ]with t = 10. The inflow boundary and other solid boundaries 01 off off satisfy a homogeneous Dirichlet condition on the contaminant concentration; the outflow boundaries and free surface boundary satisfy a homogeneous Neumann condition. The diffusivity coefficient is assumed to be constant, κ = 0.005. Thus, the Péclet number uL Pe = = 100, where the length of the inflow section is used as the characteristic Table 1 The properties of various reduced-order models t t full full N m ε ε s ROM−1 ROM−2 t t ROM−1 ROM−2 10 153 3.25e-1 4.32e-1 0.801 111 20 220 9.81e-2 9.74e-2 0.799 56 30 257 8.03e-3 6.75e-3 0.795 32 40 284 7.53e-3 7.08e-3 0.785 26 50 308 3.39e-3 2.12e-2 0.776 22 60 328 5.83e-3 6.08e-3 0.767 19 70 342 1.41e-2 1.21e-2 0.764 18 80 355 9.09e-3 8.21e-3 0.763 15 z 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 12 of 17 http://www.apjcen.com/content/1/1/2 length L = 0.5. The contaminant concentration is assumed to be zero at initial time t = 0. With this setup, we have specified our forward problem with the input parameter θ ∈ R in the range := [0, 2] × [−2, 0]. Figure 2 shows the contaminant solution c(x, t; θ) of the forward model with θ = (0.5, −0.5) at specific times. The contaminant field grows while the source is active. After the shutoff time of the source, the contaminant moves away, spreads out, and decreases in concentration due to convection and diffusion, until it flows out of the domain. The outputs of interest are the values of contaminant solution c(x, t; θ) at selected sen- sor locations in the computational domain. These sensors are located on a 4 × 4uniform grid covering the reservoir domain as shown in Figure 1b. Model order reduction To perform the model order reduction, we need to obtain the POD expansion of velocity field. Figure 3a shows all eigenvalues of the velocity field. As the number of eigenvalues increases, the significant of eigenvalues decreases, respectively. Figure 3b shows the rela- tive errors of representing the velocity field with various numbers of POD velocity modes. Eigenvalues of contaminant field 0 500 1000 1500 2000 2500 Relative error of outputs of interest 0 100 200 300 400 500 Number of samples Figure 6 Properties of POD contaminant basis vectors: (a) eigenvalues and (b) relative error analysis. σ / Σ σ ||y y || /||y|| i i r 2 2 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 13 of 17 http://www.apjcen.com/content/1/1/2 The goal of this step is to find a suitable N that can provide a fast computation of the velocity-dependent convection term and to ensure that the relative errors are acceptable. In this example, we choose the number of terms in the velocity expansion to be N = 12 −3 yields the relative error around 5 × 10 . Figure 4 illustrates the first four POD basis vectors of the velocity field. Clearly, they represent the most important characteristics of the velocity field. A comparison between the original velocity and the approximate velocity representa- tion at final simulating time is given in Figure 5. We observe that the approximate velocity profiles at x = 0.5 (Figures 5a,b) and x = 1.0 (Figures 5c,d) are able to represent well the characteristics of the original velocity profiles. To generate the snapshots needed for the POD basis to represent the contaminant con- centration, we choose S location samples in the computational domain. In this example, we generate random input locations θ ∈ D ⊂ R for source term f (x, t; θ). Table 1 shows the relative errors and computational time ratio of the reduced-order models (ROM) for different sample sizes of snapshot matrices. Here, we use the nota- tion ROM-1 to denote the reduced model mentioned in the section ‘Proper orthogonal decomposition’ and ROM-2 for the one mentioned in the section ‘Linear expansion tech- niques for velocity-dependent term’. The POD basis is chosen based on the energy as in Equation 32, with = 99.999%. We observe that while the relative errors in the out- puts are similar for both cases, the ratio of computational time differs significantly. The case ROM-2 shows speedups from 15 to 111 times, while the case ROM-1 is in fact more expensive to solve than the full forward model (t ≈ 285 seconds) .Thisisbecause full the first case consumes a lot of time in re-computing the velocity-dependent convection 1.5 ab 1.5 0.8 0.8 1 1 0.6 0.5 0.4 0.6 0.2 0.5 0.5 0.4 1 1 0.2 0.5 1 1 Z axis Z axis 0 1 X axis X axis cd 1.5 0.5 0 0 0.5 0.5 1 1 2 2 0.5 1 1 Z axis Z axis 1.5 0 0 2 2 1 1 X axis X axis Figure 7 The first 4 POD basis vectors of the contaminant field. (a) First POD contaminant basis vector. (b) Second POD contaminant basis vector. (c) Third POD contaminant basis vector. (d) Fourth POD contaminant basis vector. 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 14 of 17 http://www.apjcen.com/content/1/1/2 3 3 y y s s 2 2 y y r r 1 1 0 0 1 1 0 20 40 0 20 40 Time Time 3 3 y y s s 2 2 y y r r 1 1 0 0 1 1 0 20 40 0 20 40 Time Time 3 3 y y ideal ideal 2 2 y y real real 1 1 0 0 0 20 40 0 20 40 Time Time 3 3 y y ideal ideal 2 2 y y real real 1 1 0 0 0 20 40 0 20 40 Time Time Figure 8 A comparison of the full model (N = 1377) and ROM (m = 257) output of interest. (a) Output comparison and (b) a set of generated data at sensors 6, 7, 10 and 11.. matrix and then projecting it onto the reduced space at each time step. We choose as our reduced-order model the case ROM-2 with sample size S = 30. Figure 6a presents the eigenvalues of contaminant snapshot matrix for the sample test case S = 30. We observe that POD eigenvalues of the contaminant field are decayed very fast, and they are significantly important only at the first few hundred values. The first four POD contaminant basis vectors are shown in Figure 7. Figure 6b shows relative errors of ROMs for different numbers of POD basis vectors. With the energy taken to 99.99999%, the ROM resulted in 434 POD basis vectors with rel- −3 ative error around 4 × 10 . In the trade-off between accuracy results and computational Table 2 Estimated inverse solutions for different numbers of POD basis (%) m θ θ ε t (s) E 1 2 θ mcmc 99.0 61 0.788 -0.595 5.31e-1 1.75e+3 90.9 111 0.484 -0.491 6.95e-2 2.74e+3 99.99 178 0.451 -0.471 1.05e-2 7.55e+3 99.999 257 0.443 -0.467 2.42e-3 8.71e+4 99.9999 343 0.446 -0.466 0 5.11e+5 y y 10 6 y y 10 6 y y 11 7 y y 11 7 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 15 of 17 http://www.apjcen.com/content/1/1/2 Pairs plot parameters credible region 0.4 pdf of θ pdf of θ 0.45 0.5 0.55 0.3 0.4 0.5 0.6 Ini2 0.5 Ini3 Ini1 Ini4 1.5 1 0.5 0 0.5 1 1.5 2 2.5 Figure 9 Credible regions and historical trajectory plots of the MCMC convergence. (a) A parameter plot showing estimated source location with credible regions. (b) A contour plot of the posterior pdf. efforts (please refer to Table 1), we can choose our ROM with 99.999% of energy captured which resulting in the size of m = 257 POD basis vectors and an acceptably small error. Figure 8a shows the outputs of the full model (y) and ROM-2 (y ) with size m = 257 at selected sensor locations. The selected sensors (sensor 6, 7, 10, and 11) capture more information of the contaminant than the others for the case studied. We observe that ROM-2 is able to capture well the behavior of the full model at the sensor locations. This reduced-order model will be utilized as the forward solver in solving the inverse problem. Inverse problems MCMC is now used to solve the inverse problem for a variety of source locations using the ROM-2 solver above. By accounting for both measured noise and uncertain information in model parameters, Bayesian inference to inverse problem leads to a well-posed problem resulting in a posterior distribution of the unknown [1,2]. Once obtaining the posterior Table 3 Estimated inverse solutions for different MCMC starting points Initial θ θ θ θ t (s) 1−ini 2−ini 1−M 2−M mcmc Ini1 0.821 -0.955 0.451 -0.467 8.71e+4 Ini2 0.322 -0.451 0.458 -0.463 4.53e+4 Ini3 1.222 -0.451 0.452 -0.464 9.80e+4 Ini4 0.323 -1.313 0.456 -0.465 1.05e+5 Z 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 16 of 17 http://www.apjcen.com/content/1/1/2 distribution sample, all statistical questions related to the unknown could be answered with sample averages. To represent the behavior of uncertain variables such as wind velocity, we generate noisy data at sensor locations. The ‘real’ data (y ) is generated synthetically by adding noise real to the ideal data (y ) as in Equation 18. The noise is assumed to be Gaussian η ∼ ideal N (0, σ I). Figure 8b shows the noisy data at sensor locations with σ = 0.3. We use the snapshots with S = 30 samples above to generate several reduced-order models with different energies. To estimate the relative error of the inverse problem solu- tions (ε ), we choose the solution corresponding to the largest order as a ‘truth’ solution. We do the MCMC simulation with the starting point θ = (0.821; −0.955) .The total ini number of MCMC samples is set to N = 5000. The initial burn-in period is set to mcmc N = 500. After this stage, data is saved to compute summary statistics of source burnin locations. Table 2 shows the results of the inverse source problems with different numbers of POD basis. For the case with m = 257, the computational time is around 24 h to obtain the solution. If we use the full forward solver, with the speedup factor as given in Table 1, the computational time is estimated around 773 h or approximately 32 days. Figure 9a shows the credible regions of the source locations θ. In this figure, both a pairwise scatter plot and 1D marginal distributions are displayed. The dash line (in θ axis) and dash-dot line (in θ axis) show the probability density function of each param- eter while the solid blue regions show the 95% and 99% credible regions in which the parameters appeared the most (the blue dots). To assess the convergence of the MCMC approach, we generate four different initial points for the MH sampler. Figure 9b shows the contour plot of the posterior probability density function (the contour plot corresponds to the result of Init1) and the historical trajectories of MCMC samplers for each starting point. We observe that for all four initial points, the MCMC simulations converge to the same final solution. However, depending on the position of the initial points, the computational time varies, as given in Table 3. Conclusion This study has applied successfully the combination of a model order reduction tech- nique based on the POD and a Bayesian inference approach to solve an inverse problem that seeks to identify an uncertain contaminant location. Applying an additional POD expansion to approximate the velocity results in a reduced-order model in both state (concentration) and parameter (velocity). The resulting reduced-order model is efficient for solution of the forward problem. Solution of the Bayesian formulation of the inverse problem using the reduced-order solver is much more rapid than using the full model, yielding the probability density of the source location in reasonable computational times. The computational time of using a reduced-order model with m = 257 degrees of free- dom is about a factor of 32 times lower than using the full model with size N = 1, 377. This reduction is important in real-time water quality management applications because it reduces time cost and storage requirements. Endnote The simulations were performed on a personal computer (PC) with processor Intel(R) Core(TM)2 Duo CPU E8200 @2.66GHz 2.66GHz, RAM 3.25GB, 32-bit Operating System. 2014, 1:2 Nguyen et al. Asia Pacific Journal on Computational Engineering Page 17 of 17 http://www.apjcen.com/content/1/1/2 Author details 1 2 Singapore-MIT Alliance, National University of Singapore, Singapore 117576, Singapore. Department of Mechanical Engineering, National University of Singapore, Singapore 117576, Singapore. Massachusetts Institute of Technology, Cambridge, MA 02139-4307, USA. Received: 9 October 2013 Accepted: 23 November 2013 Published: 29 April 2014 References 1. Tarantola A (2004) Inverse problem theory and methods for model parameter estimation. SIAM, Philadelphia 2. Sivia DS, Skilling J (2006) Data analysis: a Bayesian tutorial, Second edition. Oxford University Press Inc., New York 3. Jackson C, Sen MK, Stoffa PL (2004) An efficient stochastic Bayesian approach to optimal parameter and uncertainty estimation for climate model predictions. J Climate 17: 2828–2841 4. Snodgrass MF, Kitanidis PK (1997) A geostatistical approach to contaminant source identification. Water Resour Res 34(4): 537–546 5. Woodbury AD, Ulrych TJ (2000) A full-Bayesian approach to the groundwater inverse problem for steady state flow. Water Resour Res 36(1): 159–171 6. Marzouk YM, Najm HN, Rahn LA (2007) Stochastic spectral methods for efficient Bayesian solution of inverse problems. J Comput Phys 224(2): 560–586 7. Wang J, Zabaras N (2005) Hierarchical Bayesian models for inverse problems in heat conduction. Inverse Probl 21: 183–206 8. Metropolis N, Ulam S (1949) The Monte Carlo method. J Am Stat Assoc 44: 335–341 9. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E (1953) Equations of state calculations by fast computing machines. J Chem Phys 21: 1087–1092 10. Hastings WK (1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57: 97–109 11. Gelman A, Carlin J, Stern H, Rubin D (2003) Bayesian data analysis. Second Edition, Chapman & Hall/CRC Texts in Statistical Science 12. Zienkiewicz OC, Morgan K (1983) Finite elements and approximation. Wiley, New York 13. Antoulas AC (2005) Approximation of large-scale dynamical systems. SIAM Advances in Design and Control, Philadelphia 14. Lieberman C, Willcox K, Ghattas O (2010) Parameter and state model reduction for large-scale statistical inverse problems. SIAM J Sci Comput 32(5): 2523–2542 15. Lieberman C, Willcox K (2012) Goal-oriented inference: approach, linear theory, and application to advection-diffusion. SIAM J Sci Comput 34(4): 1880–1904 16. Sirovich L (1987) Turbulence and the dynamics of coherent structures. Part 1: coherent structures. Q Appl Math 45(3): 561–571 17. Codina R, Blasco J (2000) Stabilized finite element method for the transient Navier-Stokes equations based on a pressure gradient projection. Comput Methods Appl Mech Engrg 182: 277–300 18. Ma X, Zabaras N (2008) A stabilized stochastic finite element second-order projection method for modeling natural convection in random porous media. J Comput Phys 227: 8448–8471 19. Jaynes ET (1957) Information theory and statistical mechanics. Phys Rev 104(4): 620–630 20. Wilcox CD (1993) Turbulence Modeling for CFD, DCW Industries, La Cañada, California 21. Crank J, Nicolson P (1996) A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Adv Comput Math 6(1): 207–226 doi:10.1186/2196-1166-1-2 Cite this article as: Nguyen et al.: Model order reduction for Bayesian approach to inverse problems. Asia Pacific Journal on Computational Engineering 2014 1:2. Submit your manuscript to a journal and benefi t from: 7 Convenient online submission 7 Rigorous peer review 7 Immediate publication on acceptance 7 Open access: articles freely available online 7 High visibility within the fi eld 7 Retaining the copyright to your article Submit your next manuscript at 7 springeropen.com
Journal
Asia Pacific Journal on Computational Engineering – Springer Journals
Published: Apr 29, 2014