Generalized Empirical Interpolation Method With H 1 Regularization: Application to Nuclear Reactor Physics

The generalized empirical interpolation method (GEIM) can be used to estimate the physical field by combining observation data acquired from the physical system itself and a reduced model of the underlying physical system. In presence of observation noise, the estimation error of the GEIM is blurred even diverged. We propose to address this issue by imposing a smooth constraint, namely, to constrain the H 1 semi-norm of the reconstructed field of the reduced model. The efficiency of the approach, which we will call the H 1 regularization GEIM (R-GEIM), is illustrated by numerical experiments of a typical IAEA benchmark problem in nuclear reactor physics. A theoretical analysis of the proposed R-GEIM will be presented in future works.


INTRODUCTION AND PRELIMINARIES
In nuclear reactor simulations, data assimilation (DA) with reduced basis (RB) enables the reconstruction of the physical field, e.g., for constructing a fast/thermal flux and power field within a nuclear core in an optimal way based on neutronic transport/diffusion model and observations (Gong et al., 2016;Argaud et al., 2017a;Gong et al., 2017;Argaud et al., 2018;. In practice, the existing methods are based on a reduced basis, however, this approach is not robust with respect to observation noise and there are some additional constraints or regularization on the low-dimensional subspaces, i.e., the related coefficients, have been proposed as possible remedies in several recent works (Argaud et al., 2017b;Gong et al., 2020;Gong et al., 2021). The idea of introducing box constraints was originally introduced in Argaud et al. (2017a) to stabilize the generalized empirical interpolation method in the presence of noise (Maday and Mula, 2013). The same idea has been applied to the POD basis and the background space (Gong et al., 2019) of the so-called parametrized-background data-weak (PBDW) data assimilation (Maday et al., 2015a). Recently, the regularization of the GEIM/POD coefficients has been studied in Gong et al. (2021). The corresponding theoretical analysis can be found in , Herzet et al. (2018), and Gong et al. (2019). This article introduces H 1 regularization schemes for the approximation, and numerical experiments in nuclear reactor physics indicate its potential to address this obstruction.
Our goal is to approximate the physical state u from a given compact set M ⊂ X (manifold), which represents the possible state of a physical system taking place in Ω. Where X is a Banach space over a domain Ω ⊂ R d (d ≥ 1) being equipped with the norm . X . In the framework of data assimilation with reduced basis, any u ∈ M can be estimated by combining two parts. The first term is a certain amount (m) of observation of u acquired directly from sensors of the underline physical system, which can be represented by a combination of linear functionals of X ′ (the dual space of X ) evaluated on u. The second term is the use of a family of (reduced) subspaces Z n of finite dimension n which is assumed to approximate well with the manifold M in a given accuracy.
The algorithms used to build the reduced subspace {Z n } n and find appropriate linear functionals have already been reported in the community of reduced modeling [see Maday and Mula, 2013;Maday et al., 2015b;Maday et al., 2015c;Maday et al., 2016]. Note that even if this is not necessary in the previous statements, the construction of the reduced spaces Z n could be recursive, i.e., we have Z n−1 ⊂ Z n . The field u ∈ M can be approximated by interpolation  or data assimilation Gong et al., 2019).
For reading convenience, let us first introduce some notations used throughout this article. We first introduces the standard L 2 (Ω) or H Hilbert space over the special domain Ω ⊂ R d equipped with an inner product (w, v) L 2 (Ω) ≡ Ω wvdx and the For a given Hilbert space U and the related dual space U′, the Riesz operator R U : U′ → U satisfies: for any given ℓ ∈ U′, we have Let us denote by u (r; μ) the solution of a parameter-dependent partial differential equation (PDE) set on Ω and on a closed parametric domain D ⊂ R p . For any given μ ∈ D, the physical field u (r; μ) belongs to U ⊂ L 2 (Ω) or H 1 (Ω), a functional space derived from the PDE. We call the set M ≡ {u(r; μ); μ ∈ D} of all parameter-dependent solution manifold. Let L M : U → R M be the vector-valued observation functional L M (u) (ℓ 1 (u), . . . , ℓ M (u)) T of u ∈ U.

FIELD RECONSTRUCTION WITH REGULARIZATION
Our goal in this work is to infer any state u ∈ U over a spacial domain Ω ∈ R d given only some corresponding noisy observations y (y obs 1 , . . . , y obs M ) T . This empirical learning problem from a limited data set is always underdetermined. In general, with observation noise, a regularization term R (u) is added to the loss function, and then a general convex model fitting problem can be written in the form: where V: R M → R is a convex loss function that describes the cost of predicting u when the observation is y. λ is a parameter which presents the importance of the regularization term. R (u) is usually a convex regularization function to impose a penalty on the complexity of u through some prior knowledge.

H 1 Regularization
The goal of regularization is to prevent overfitting or to denoise in mathematics and particularly in the fields of inverse problems (Ivanov, 1976;Andreui et al., 1977;Balas, 1995;Arnold, 1998;Vladimir, 2012;Benning and Burger, 2018), by introducing additional information in order to solve ill-posed problems. From a Bayesian (James Press, 1989) point of view, regularization techniques correspond to introduce some prior distributions on model parameters. The general choice of R (u) is a norm-like form u 2 χ , where χ represents different kinds of norm depending on the underlying application. The simplest choice of R (u) is L 2 norm, say, u 2 L 2 (Ω) , which has been well studied in the literature, either from a theoretical or algorithmic point of view. The regularization is also called Tikhonov regularization (Andreui et al., 1977), which is essentially a trade-off between fitting the data and reducing the norm of the solution. For some real-world problems, there has been much interest in alternative regularization terms. For example, total variation (TV) regularization, R(u) ∇u L 1 (Ω) , is popular in image reconstruction or other domains (Rudin et al., 1992;Rudin et al., 1992;Chan et al., 1997;Chan and Tai, 2003;Wachsmuth and Wachsmuth, 2011;De los Reyes and Schönlieb, 2012). By using a L 1 norm, sharp edges would be allowed as the penalty is finite, and it also allows discontinuous controls which can be important in certain applications.
If one would like to impose a smooth control, the H 1 seminorm can be used: Examples can be found in the context of parameter estimation problems (Keung and Zou, 1998;Cai et al., 2008;Wilson et al., 2009;Barker et al., 2016), image-deblurring (Chan et al., 1997;Li et al., 2010;Cimrák and Melicher, 2012), image reconstruction (Ng et al., 2000), and flow control (Heinkenschloss, 1998;Collis et al., 2001). The work in Van Den Doel et al. (2012) shows that the proposed H 1 semi-norm regularization performs better than L 1 regularization cousin, total variation, for problems with very noisy data due to the smooth nature of controlled variables. The authors in Srikant did a comparison of H 1 and TV regularization methods and also studied the shortcomings and limitations of some of the implementations schemes, such as a Gaussian filter. H 1 regularization would perform well over uniform regions in the domain but would perform poorly over edges. Furthermore, to solve the PDE-constrained optimization problem as reported in Haber and Hanson (2007), the authors suggested a synthetic regularization functional of the form: where the parameter c can be adapted. Note that this synthetic regularization is now commonly used to solve the ill-posed inverse problems.

Generalized Empirical Interpolation Method
Recall that our goal is to estimate the state u true [μ] ∈ U of a physical system for a given parameter μ ∈ D, by using a parameterized best-knowledge model and M (potentially noisy) observations. The first step is to choose a sequence of n-dimensional subspaces {Z n } n such that the best approximation of any given u true [μ] in the space Z n converges to zero when n goes to infinity, i.e., for an acceptable tolerance ϵ Z . We further assume that the selected subspaces satisfy In other words, we choose the subspaces such that the most dominant physical system is well represented for a relatively small n. In particular, these subspaces may be constructed through the application of model reduction methods to a parameterized PDE.
The second step is to model the data acquisition procedure. Given a system in of a parameter μ ∈ D, we assume the observations are of the form where y obs m [μ] is the value of the mth observation, ℓ m is the linear functional associated with the mth sensor, and e m is the observation noise associated with the mth sensor. The detailed form of the functional ℓ m depends on the specific sensor used to acquire data. For example, if the sensor measures a local value of the state, then we may model the observation value as Gaussian convolution ) is a Gaussian distributed function to present the response of the sensor for a given physical field v, and x c m ∈ R d is the center of the sensor in the special domain Ω, and r m ∈ R >0 is the physical width of the sensor. In particular, the localized sensor is of interest in this work.
We assume that e m is independent and identically distributed (IID), and with a density of p m on R. In practice, the mean and covariance of the observation data acquired are more readily quantifiable than the distribution p m . Thus, we assume the mean and the covariance of the distribution exist and make the following assumptions on the noise term: By running the greedy algorithm of the so-called generalized empirical interpolation method [GEIM (Maday and Mula, 2013)], a set of basis {q n } n is generated and spanned the reduced space Z n span{q 1 , . . . , q n }. Then, the generalized interpolation process is well defined as follows: (2.8) With noisy observations, GEIM is, however, not robust with respect to observation noise (Argaud et al., 2017b), and in that work, a so-called constrained stabilized GEIM (CS-GEIM) by using a constrained least squares approximation was proposed to address this obstruction, where numerical experiments indicate its potential.

H 1 REGULARIZATION FORMULATION OF GEIM
Now, we state the H 1 regularization scheme for GEIM (R-GEIM). Given a reduced space Z N span{q 1 , . . . , q N } ⊂ U of dimension N spanned by N basis {q i } N i 1 and M measurement functionals L M ≔ (ℓ 1 , . . . , ℓ M ) T and the corresponding noisy measurements y (y obs 1 , . . . , y obs M ) T , M ≥ N, then the reconstruction problem from measurements is: find u ∈ Z N such that: where V: R M → R is loss function that evaluates the cost of estimating u giving the observation y which depends on the underlying application. The symbol "argmin · " is argument of the minimum, thus argmin u∈ZN f(u) is the value of u for which f(u) attains its minimum. The parameter ξ is a trade-off factor between the regularization term and the loss function term. Furthermore, if we have no information about the noise distribution, three proposed typical forms of V could be L(u) − y 2 l ,, where l ∞, 1, 2. For the above R-GEIM, we have the following remarks: • The basis and measurements of the chosen scheme could be based on GEIM, POD, or any other approach. • Later, we will show that l 2 corresponds to the least squares method, and l 1 corresponds to the least absolute deviations (LAD) (Bloomfield and Steiger, 1980). Compared to the traditional least squares method, the LAD is much robust and finds its applications in many areas. • By using the H 1 regularization for the reduced basis field reconstruction, the first assumption is that the field u is in H 1 space; for the most regular physical problem, this condition is satisfied automatically, and the H 1 regularization term is a kind of smooth control for the underlying field reconstruction problem.
If we have no prior information about noise, a commonly used way to formulize Eq. 3.1 is taking 2 norm for the loss function term, we have the following results: R-GEIM: Given a reduced space Z N span{q 1 , . . . , q N } ⊂ U of dimension N spanned by N basis {q i } N i 1 and M(≥ N) measurement functionals L M ≔ (ℓ 1 , . . . , ℓ M ) T and the corresponding noisy measurements y (y obs 1 , . . . , y obs M ) T , then the reconstruction problem from measurements is: find u ∈ Z N such that: Let M be an M × N full-column rank matrix with elements M i,j ℓ i (q j ), i 1, . . . , M, j 1, . . . , N and N be an N × N matrix with elements N i,j (∇q i , ∇q j ), i, j 1, . . . , N, then the algebraic form of Eq. 3.2 is: Later, we will show this is also the algebraic formulation for Gaussian noise with covariance D. If we can make use of the prior information of noise, we have the following remark: Remark. Let r(u) be the bias of the reduced model from the truth L(u) − y e + r(u). By using maximum likelihood (ML) estimation for the following common noise densities, we have: • Uniform noise, when the noise term e m is uniformly independent and identically distributed on ( − e 0 , e 0 ), then the reconstruction problem is: find u ∈ Z N such that: (3.7) If the noise bounds are different for different measurements, the constraint in Eq. 3.7 becomes |ℓ m (u) − y obs m | ≤ e m + |r m (u)|, m 1, . . . , M.
• Gaussian noise, when the noise e m is Gaussian with the zero mean and covariance matrix D, then the reconstruction problem is: find u ∈ Z N such that u argmin u∈ZN (L(u) − y) T D −1 (L(u) − y) + Mξ ∇u 2 L 2 (Ω) . (3.8) • Laplacian noise, when the noise e m is Laplacian independent, identically distributed with density p(e) 1 2e0 e −|e|/e 0 , then the reconstruction problem is: find We refer readers to Boyd and Vandenberghe (2004) for further theoretical analysis on this remark. Through this remark, the physical means of the term L(u) − y 2 l in Eq. 3.1 for different l is easier to understand. The ∞, 1, and 2 norms interpret the maximum likelihood estimation with a noise density, that is, uniform, Laplacian, and Gaussian, respectively. Considering for the most engineering problems, the noise density is Gaussian and also bounded, and thus we only present numerical results of uniform noise and Gaussian noise in this work.
Another remark is that, in this work, numerical results are illustrated based on GEIM, more precisely. The reduced basis in Eq. 3.2 is derived with GEIM, but this regularization is fit for the POD basis or the basis selected from the greedy reduced basis method (Grepl et al., 2007) without any modification.

NUMERICAL RESULTS AND ANALYSIS
In this section, we illustrate the performance of R-GEIM on a typical benchmark problem in nuclear reactor physics. The test example is adapted based on the classical 2D IAEA benchmark problem (Benchmark Problem Book, 1977;Gong et al., 2017); the geometry of the 2D core is shown in Figure 1. This problem represents the mid-plane z 190 cm of the 3D IAEA benchmark problem, that is used by references (Theler et al., 2011). The reactor spacial domain is Ω region (1, 2, 3, 4). The core and control regions are Ω core region (1, 2, 3) and Ω control region (3), respectively. We consider the value of Σ a,2 | Ω1 , Σ a,2 | Ω2 , and Σ a,2 | Ω 3 in the core region Ω 1,2,3 as a parameter (so p 3 and μ [Σ a,2 | Ω 1 , Σ a,2 | Ω 2 , Σ a,2 | Ω 3 ]).
We assume that Σ a,2 | Ω i ∈ [0.080, 0.150] for i 1, 2, 3. The rest of the coefficients of the diffusion model are fixed to the values indicated in Table 1 of Gong et al. (2017).
The neutronic field (fast and thermal flux, power distribution) is derived by solving two group diffusion equations. The numerical algorithm is implemented by employing the free finite elements solver FreeFem++ (Hecht, 2012). The norm · L 2 (Ω) is induced by the inner product (w, v) L 2 (Ω) Ω wvdx, and the semi-norm · H 1 (Ω) is induced by the inner product (w, v) H 1 (Ω) Ω ∇w∇vdx. The measurement we employed here is same as Eq. 2.7 with r m 1 cm and set {w i } N i 1 being the corresponding Riesz representation with H 1 inner product. Finally, we set the finite element size Frontiers in Energy Research | www.frontiersin.org January 2022 | Volume 9 | Article 804018 to d 0.1 cm, which is enough for our analysis. We refer to Argaud et al. (2017b) for detailed implementation of this problem with FreeFem++. The regularization factor ξ is essential for R-GEIM. It can significantly affect the reconstruction error of R-GEIM, and if they are incorrectly specified then the field reconstructed with R-GEIM is suboptimal. We show the variation of the errors in L ∞ , H 1 , and L 2 norms for R-GEIM with respect to different regularization factors ξ in Figure 2. The dimension of reduced basis is fixed to n 80, the number of measurements is set to m FIGURE 1 | Geometry of a 2D IAEA benchmark. Upper octant: region assignments; lower octant: fuel assembly identification [from reference (Benchmark Problem Book, 1977;Theler et al., 2011)].
FIGURE 2 | Variation of the errors in L ∞ , H 1 , and L 2 norms for R-GEIM with respect to the different regularization factor ξ. The reduced dimension is n 80, number of measurement is m 2n, and the observation noise is uniformly distributed with noise levels 10 −2 (A) and 10 −3 (B). 2n, and the observation noise is uniformly distributed, with a noise level σ 10 −2 , 10 −3 . It can be observed that the optimal ξ is different for different error metrics. For the errors evaluated in L 2 norm, the optimal ξ op ∼ 0.1, and for L ∞ norm or H 1 norm, the optimal ξ op ∼ 1. In the left of this work, we fix ξ to be the optimal value and evaluate the errors in L 2 norm and H 1 norm, which reflect the average error of the reconstructed field itself and its gradient. This section illustrates the behavior of GEIM, CS-GEIM, and R-GEIM, in case of noisy observations. We first show the variation of the errors in L 2 norm and H 1 norm for GEIM,   Frontiers in Energy Research | www.frontiersin.org January 2022 | Volume 9 | Article 804018 CS-GEIM, and R-GEIM with respect to different reduced dimensions n in Figure 3. The observation noise is assumed to be uniformly distributed, with a noise level σ 10 −2 . The number of measurements is m 2n. Figure 4 illustrates the case for the noise level σ 10 −3 . The cases with Gaussian-distributed observation noise are shown in Figure 5 for the noise level σ 10 −2 and in Figure 6 for the noise level 10 −3 .
From these figures, we can conclude that R-GEIM shows a better stability performance in the case of noisy measurement. The accuracy can be as good as CS-GEIM, but the R-GEIM algorithm is much simpler, with relatively low computational cost; the main cost for the online stage is to solve the matrix system Eq. 3.4. But for CS-GEIM, the relative complex constrained quadratic programming problem has to be solved.

CONCLUSION AND FUTURE WORKS
The traditional generalized empirical interpolation method is well studied for data assimilation in many domains. However, this reduced modeling-based data assimilation method is not robust with respect to observation noise. We propose addressing this issue by imposing a smooth constraint, namely, an H 1 semi-norm of the reconstructed field to involve some prior knowledge of the noise. The efficiency of the approach, which we call R-GEIM, is illustrated by an IAEA benchmark numerical experiment, dealing with the reconstruction of the neutronic field derived from neutron diffusion equations in nuclear reactor physics. With H 1 regularization, the behavior of the reconstruction is improved in the case of noisy observation. Further works are ongoing: i) mathematical analysis of the stable and accurate behavior of this regularization approach and ii) the regularization trade-off factor will be studied to give an outline on how to choose these factors for generic problems.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HG contributed to conceptualization, methodology, software, data curation, and writing-original draft and editing. ZC contributed to conceptualization. QL contributed to conceptualization, supervision, and review.

FUNDING
This work is supported by the National Natural Science Foundation of China (Grant No. 11905216 and Grant No. 12175220).