Stabilization of (G)EIM in presence of measurement noise: application to nuclear reactor physics

General overview and contribution of the paper

For the sake of clarity, we shall start by formulating the goal of the paper in general terms containing statements that will be clarified in the forthcoming sections.

Under this hypothesis on the decay of $d_{n}(\mathcal{S},\mathcal{X})$, one can in principle build a sequence $\{X_{n}\}_{n}$ s.t. $\operatorname{dist}(\mathcal{S},X_{n})\coloneqq\max_{u\in\mathcal{S}}\min_{v\in X}\|u-v\|_{\mathcal{X}}\leq\varepsilon,$ where $\dim(X_{n})=n\equiv n(\varepsilon)$ is moderate.

Algorithms to build $\{X_{n}\}_{n}$ (or at least the first spaces in the sequence allowing to approximate beyond a given accuracy) and find appropriate linear functionals have been proposed in the community of reduced modeling (see MM2013 ; MPPY2014 ; MMT2016 ). Note that, even if this is not required in the previous statements, the construction of the spaces $X_{n}$ is then recursive, i.e. $X_{n-1}\subset X_{n}$. There, the approximation of $f\in\mathcal{S}$ is done by interpolation or related approximations. The methods (in practice mainly based on a greedy procedure) are however not robust with respect to noise in the measurements and this paper introduces a constrained least squares approximation for which numerical experiments indicate its potential to address this obstruction.

Mathematical setting

Let us assume that $\mathcal{M}=\overline{\hbox{span}(\mathcal{S})}$ (where the $\overline{}\mathcal{B}$ denotes the closure in $\mathcal{X}$ of the set $\mathcal{B}$) admits a Schauder basis $\{q_{i}\}_{i}$, i.e., for every $f\in\mathcal{X}$ there exists a unique sequence $\{c_{i}(f)\}$ of scalars such that $\lim_{n\to\infty}\|f-\sum_{i=1}^{n}c_{i}(f)q_{i}\|_{\mathcal{X}}=0.$ For every $n\geq 1$, we define the $n$-dimensional subspace $X_{n}\coloneqq\operatorname*{span}\{q_{1},\dots,q_{n}\}.$ Let us formulate in a different manner the hypothesis made involving the Kolmogorov $n$-width of $\mathcal{S}$ in $\mathcal{X}$: let us assume that the error in approximating the functions of $\mathcal{S}$ in $X_{n}$ is $\max_{f\in\mathcal{S}}\operatorname{dist}(f,X_{n})\leq\varepsilon_{n},$ where the sequence $(\varepsilon_{n})_{n}$ decays at a nice rate with $n$.

Let now $\{\lambda_{i}\}$ be the set of linear functionals of $\mathcal{X}^{\prime}$ (of unity norm in $\mathcal{X}^{\prime}$) such that for every $n\geq 1$, $\{\lambda_{1},\dots,\lambda_{n}\}$ and $\{q_{1},\dots,q_{n}\}$ are such that, for every $n\geq 1$ and every $1\leq j\leq n$, $\forall i,\ 1\leq j\leq i\quad\lambda_{j}(q_{i})=\delta_{i,j}.$ For any $n\geq 1$, we can now define a (generalized) interpolation operator $\mathcal{J}_{n}:\mathcal{X}\to X_{n}$ such that for all $f\in\mathcal{X}$ $\lambda_{i}(f)=\lambda_{i}\left(\mathcal{J}_{n}(f)\right),\quad i\in\{1,\dots,n\}.$ By construction, for any $n\geq 1$ and any $f\in\mathcal{X}$, $\mathcal{J}_{n}(f)=\mathcal{J}_{n-1}(f)+c_{n}(f)q_{n}.$ where $c_{n}(f)=\lambda_{n}\left(f-\mathcal{J}_{n-1}(f)\right)$ and, for notational coherence, we set $\mathcal{J}_{0}=0$.

Using $\mathcal{J}_{n}$ to approximate the functions of $\mathcal{S}$ yields the error bound

where $\Lambda_{n}\coloneqq\sup_{f\in\mathcal{X}}\|\mathcal{J}_{n}(f)\|_{\mathcal{X}}/\|f\|_{\mathcal{X}}$ is the Lebesgue constant. The value of $\Lambda_{n}$ diverges at a certain rate so the behavior of $\max_{f\in\mathcal{S}}\|f-\mathcal{J}_{n}[f]\|$ with the dimension is dictated by the trade-off between the rate of divergence of $(\Lambda_{n})$ (that is generally slow) and the convergence of $(\varepsilon_{n})$ (that is generally very fast). Also, for any $f\in\mathcal{S}$

where $\Lambda_{0}=0$ and $\varepsilon_{0}=\max_{f\in\mathcal{S}}\|f\|$.

For any $\alpha>0$, let us define the cone $\mathcal{K}_{n}(\alpha)\coloneqq\{v\in V_{n}\ :\ v=\sum_{i=1}^{n}c_{i}q_{i}\,\ |c_{i}|\leq\alpha(1+\Lambda_{i-1})\varepsilon_{i-1}\}.$ We have for any $n\geq 1$ and any $f\in\mathcal{S}$ $\mathcal{J}_{n}(f)\in\mathcal{K}_{n}(1)$. In presence of noise in the measurements, we assume that we receive values $\eta_{1}(f),\dots,\eta_{n}(f)$ such that $\eta_{i}(f)\sim\lambda_{i}(f)+\mathcal{N}(0,\sigma^{2})$ for $i\in\{1,\dots,n\}$. Interpolating from these values yields an element in $X_{n}$ denoted as $\mathcal{J}_{n}(f;\mathcal{N})$ that satisfies blurred error bound with respect to (2) that, depending on the precise definition of $\sigma$ and the norm of $\mathcal{X}$ can be

The second term of the bound diverges as $n$ increases and shows that the method is not asymptotically robust in presence of noise. An illustration of this can be found in the numerical results below. This motivates the search for other methods which would ideally yield a bound of the form $\left(1+\Lambda_{n}\right)\varepsilon_{n}+\sigma$ and for which the error is asymptotically at the level of the noise $\sigma$.

In this paper, we are running the greedy algorithm of the so-called Generalized Empirical Interpolation Method (GEIM, MM2013 ), to generate a basis $\{q_{i}\}_{i}$ and the linear functionals $\{\lambda_{i}\}_{i}$. This method is reported to have a nice behavior for the Lebesgue constant, at least in case it is trained on a set $\mathcal{S}$ with small Kolmogorov dimension (see MMT2016 ). This approach allows an empirical optimal selections of the positions of the sensors that provide (in case where no noise pollutes the measures) a stable representation of the physical system. The precise algorithm is documented elsewhere (see MM2013 and GABM2016 ). Then, we approximate any $f\in\mathcal{S}$ with the function $A_{n}(f)$ defined in (4) with $\alpha=2$. We call this scheme Constrained Stabilized GEIM (CS-GEIM).

Note that the above approach could also be used with a POD approach to provide the imbedded spaces $\{X_{n}\}_{n}$ (that are more expensive to provide than the greedy GEIM approach but are more accurate) and well chosen linear functionals $\lambda_{n,i}$ the choice of which infer on the behavior of the Lebesgue constant $\Lambda_{n}$.

Numerical results

For the physical problem that we consider in this paper, the model is the two group neutron diffusion equation : the flux $\phi$ has two energy groups $\phi=(\phi_{1},\phi_{2})$. Index 1 denotes the high energy group and 2 the thermal energy one. These are modeled by the following parameter dependent PDE model :

here $k_{\text{eff}}$ is the so-called multiplication factor and is not a data but an unknown of the problem We omit here the technical details on the meaning of $k_{\text{eff}}$ and refer to general references like Hebert2009 ., and the given parameters are

$D_{i}$ is the diffusion coefficient of group $i$ with $i\in\{1,2\}$.

$\Sigma_{a,i}$ is the macroscopic absorption cross section of group $i$.

$\Sigma_{s,1\to 2}$ is the macroscopic scattering cross section from group 1 to 2.

$\Sigma_{f,i}$ is the macroscopic fission cross section of group $i$.

$\nu$ is the average number of neutrons emitted per fission.

$\chi_{i}$ is the fission spectrum of group $i$.

they are condensed in $\mu=\{D_{1},D_{2},\Sigma_{a,1},\Sigma_{a,2},\Sigma_{s,1\to 2},\nu\Sigma_{f,1},\nu\Sigma_{f,2},\chi_{1},\chi_{2}\}$.

We assume that the parameters of our diffusion model range in, say,

then $\mathcal{D}:=[D_{1,\min},D_{1,\max}]\times\dots\times[\chi_{2,\min},\chi_{2,\max}]$ is the set of all parameters and the set of all possible states of the flux is given by

where the power $P(\mu)$ is defined from $(\phi_{1},\phi_{2})(\mu)$ as $P(\mu)(x):=\nu\Sigma_{f,1}\phi_{1}(\mu)(x)+\nu\Sigma_{f,2}\phi_{2}(\mu)(x),\quad\forall x\in\Omega$ We assume (see CD2014 for elements sustaining this hypothesis) that the Kolmogorov-width decays rapidly, hence, it is possible to approximate all the states of the flux (given by $\mathcal{S}$) with an accuracy $\varepsilon$ in well-chosen subspaces $X_{n}\subset\mathcal{X}$ of relatively small dimension $n(\varepsilon)$.

To ensure enough stability in the reconstruction and minimize the approximation error, it is necessary to find the optimal placement of the sensors in the core. The selection is done with GEIM. If we denote $\sigma(\phi_{i},x),\quad i\in\{1,2\},$ the measurement of $\phi_{i}$ at a position $x\in\Omega$ by a certain sensor, this measurement can be modeled by a local average over $\phi_{i}$ centered at $x\in\Omega$. Another possibility is to directly assume that the value $\phi_{i}(x)$ at point $x$ is $\sigma(\phi_{i},x)$In the following part of this work, we directly assume that the value $\phi_{i}(x)$ at point $x$ is $\sigma(\phi_{i},x)$ as measurement.. Note that, in principle, the measurement could depend on other parameters apart from the position. We could imagine for instance that we have sensors with different types of accuracy or different physical properties. This flexibility is not included in the current notation but the reader will be able to extrapolate from the current explanations.

A specificity of the approach here is that $\mathcal{S}$ is composed of vectorial functions $(\phi_{1},\phi_{2},P)(\mu)$. We deliberately choose to take measurements only on one of the components (say $\phi_{2}$) and thus reconstruct the whole field $\phi_{1},\ \phi_{2}$ and $P$ with the only knowledge of thermal flux measurements.

Another specificity of our approach is on the spatial location of the measurements. We consider two cases:

Case I: the sensors can be placed at any point in the domain of definition of $\phi_{2}$.

Case II: the admissible sensor locations are restricted to be deployed in a restricted part of that domain;

We have already reported in ABGMM2016 that these two specificities are well supported by the (G)EIM approach, as long as the greedy method is taught to achieve the goal of reconstructing the whole field $(\phi_{1},\phi_{2},P)(\mu)$.

Our aim here is to show that the noise can be controlled through our Constrained Stabilized (G)EIM approach.

2 Description of the PARCS 2D IAEA benchmark

We consider the classical 2D IAEA Benchmark Problem Benchmark , the core geometry which can be seen in figure 1. The problem conditions and the requested results are stated in page 437 of reference Benchmark . It is identified with the code 11-A2, and its descriptive title is Two-dimensional LWR Problem , also known as 2D IAEA Benchmark Problem . This problem represents the mid-plane $z=190~{}cm$ of the 3D IAEA Benchmark Problem, that is used by references PARCS and show in application within 2D-IAEA-Benchmark .

The reactor domain is $\Omega=\text{region}(1,2,3,4)$. The core and the reflector are $\Omega_{\text{core}}=\text{region}(1,2,3)$ and $\Omega_{\text{refl}}=\text{region}(4)$ respectively. We consider only the value of $D_{1}|_{\Omega_{\text{refl}}}$ in the reflector $\Omega_{\text{refl}}$ as a parameter (so $p=1$ and $\mu=D_{1}|_{\Omega_{\text{refl}}}$). We assume that $D_{1}|_{\Omega_{\text{refl}}}\in[1.0,3.0]$. The rest of the coefficients of the diffusion model (5) (including $D_{1}|_{\Omega_{\text{core}}}$) are fixed to the values indicated on table 3.2. In principle, one could also consider these coefficients as parameters but we have decided to focus only on $D_{1}|_{\Omega_{\text{refl}}}$ because of its crucial role in the physical estate of the core: its variation can be understood as a change in the boundary conditions in $\Omega_{\text{core}}$ which, up to a certain extent, allows to compensate the bias of the diffusion model with respect to reality. We shall report in a future paper more extended variations of the parameters.

Figure 5(a) shows the behavior of the Lebesgue constants in both cases. We can find that 1) the Lebesgue constant increases with GEIM interpolation function dimension, 2) if detectors are limited in a domain part (Case II), the Lebesgue constant gets worse, as an effect of the extrapolation that is required here, nevertheless the increase is still moderate.

Figure 5(b) shows the coefficients upper limits described as $r_{n}(x_{n},\mu_{n})\equiv(1+\Lambda_{n-1})\varepsilon_{n-1}$ (see (2)) (so $|c_{n}|\leq r_{n}(x_{n},\mu_{n})$) for Case I and Case II, which decreases quickly with $n$.

We still focus on the 300 parameters $\mathcal{D}^{(test)}$, and compute the errors with equation (8), for each test case, we perform the interpolation process with CS-GEIM a number of times. Figure 7 shows the averaged $L^{2}(\Omega)$ norm for the decay of $e^{\text{(test)}}_{n}(\phi_{1})$, $e^{\text{(test)}}_{n}(\phi_{2})$ and $e^{\text{(test)}}_{n}(P)$, with noise amplitude $10^{-2}$ for Case I and Case II. For different measurement noise amplitude, the averaged errors in $L^{2}(\Omega)$ norm, $L^{\infty}(\Omega)$ norm and $H^{1}(\Omega)$ norm are shown in figure 8, figure 9 and figure 10 respectively, for Case I and Case II. The main conclusions are: in the noisy case, i) CS-GEIM improves the interpolation, with the error comparable to the noise input level, ii) in extrapolation case, CS-GEIM reduces the interpolation error dramatically, which extends GEIM practical use.

If we take more measurements with fixed number of interpolation functions, the ratio $n/m$ of the number of measurements $n$ to the number of interpolation functions $m$ increases, so it is expected to have the same effect than to repeat independent measure at the same point in order to measure the evaluation of the measure. We consider the analytical function $g(x,\mu)\equiv\mathcal{V}((x_{1},x_{2});(\mu_{1},\mu_{2}))\equiv((x_{1}-\mu_{1})^{2}+(x_{2}-\mu_{2})^{2})^{-1/2}$ for $x\in\Omega\equiv]0,1[^{2}$ and $\mu\in\mathcal{D}\equiv[-1,-0.01]^{2}$; we choose for $\mathcal{D}^{(training)}$ a uniform discretization sample of 400 pointsWe replace the synthetic neutron problem here by the above analytical function so as to be able to have a more thorough and extensive numerical analysis. Then we change the ratio $n/m$ of the number of measurements $n$ to the number of interpolation functions $m$ with CS-GEIM process, figure 12 also shows the error converges with $\sim n^{-\frac{1}{2}}$.

We have presented some results obtained the Empirical Interpolation Method (EIM) for the reconstruction of the whole field, solution to a simple but representative problem in nuclear reactor physics as an example of a set of parameterized functions. With EIM, a high accuracy can be got in reconstructing the physical fields, and also a better sensors deployment is proposed with which most information can be extracted in a given precision even if only part of the field (either in space or in component) is omitted in the measurement process. Then an improved Empirical Interpolation Method (CS-(G)EIM) is proposed. With CS-(G)EIM, i) the behavior of the interpolant is improved when measurements suffer from noise, ii) the error is dramatically improved in noisy extrapolation case, iii) it is possible to decrease the error by increasing the number of measurements.

Further works and perspective are ongoing: i) mathematical analysis of the stable and accurate behavior of this stabilized approach, ii) in this work, our first assumption is the model is perfect (i.e. we work on in silico solutions, a broader class of methods which couple reduced models with measured data named PBDW MPPY2015 are able to correct the bias of the model and use real data.