Data Assimilation Applied to a Simple Inverse Radiation Transport ...
Data Assimilation Applied to a Simple Inverse Radiation Transport Problem Joshua M. Hykes, Dan G. Cacuci, and John Mattingly Department of Nuclear Engineering North Carolina State University Raleigh, NC, 27695 [email protected], [email protected], and [email protected], November 2011
INTRODUCTION Implicit methods for inverse radiation transport problems fit a transport model to an experiment to estimate properties of the source measured in the experiment. Several solution algorithms have been applied, and others are still being developed. One successful solution method is nonlinear regression, where typically the Levenberg-Marquardt ( L M ) algorithm is used to minimize a chi-squared metric of the difference between the model prediction and the experiment. In this work, we apply a data assimilation ( D A ) method to an inverse transport problem; the main advantage of D A is that it is designed to accept uncertainty in the inputs and to provide uncertainty estimates for the results. We compare L M to D A using a shielded uranium sphere as a simple test problem, first published by Favorite [1]. In Ref. [2], Favorite and Ketcheson demonstrated the use of L M for this problem. (In contrast to the current work, they used adjoint solutions to compute derivatives, and they only tracked the uncollided spectrum.) As shown in Figure 1, the uranium sphere is surrounded by concentric spherical shells of void, lead, and aluminum. Following Favorite’s example, we neglect neutron emissions by the uranium and analyze only its gamma emissions. The gammas are counted by an idealized “detector” beyond the problem’s outer boundary. The goal is to find the shell thicknesses to optimally reproduce the gamma spectrum. The unknowns we solved for are the four layer thicknesses, ~ = [8.741, 3.659, 0.5, 0.3] (all in cm). whose true values are ∆ρ
1
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
Methods Nonlinear Least Squares Procedure For comparison purposes, we implemented an L M inverse solver in a similar manner as Mattingly and Mitchell [3]. We minimize the χ2 functional ~ = χ2 (∆ρ)
£ ¤ ~ −le 2 m l c (∆ρ) X i i σ2i
i =1
+
³ ´2 0 n ∆ρ j − ∆ρ X j w 2j
j =1
,
(1)
where m = number of leakage energy groups (21 here), n = number of thickness unknowns (4 here), l ic = computed leakage in group i , l ie = experimentally measured leakage in group i , σ2i = experimental variance of l ie , ∆ρ j = j th shell thickness unknown, ∆ρ 0j = j th shell thickness prior estimate, and w j = uncertainty weighting of j th prior estimate. With this definition of χ2 , the covariance of the final parameters can be estimated as (JT J)−1 , where the Jacobian is defined as Ji j =
c 1 ∂l i σi ∂∆ρ j
.
Data Assimilation Procedure Probability theory and Bayes’ theorem provide a general approach to solve inverse problems (for example, see Tarantola [4]). Although the most general description of an inverse solution is a probability distribution function over the unknowns, this full description is Aluminum Lead Uranium Void
8.741 12.9 13.2
Figure 1: A cross section of Favorite’s shielded uranium sphere. All dimensions are in cm.
2
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
usually too complex to compute or easily comprehend. Thus, a typical simplification is to track the means and covariances of these distributions. One example of this approach is the data assimilation methods presented by Cacuci and Ionescu-Bujor [5]. Their posterior distribution is a multivariate Gaussian, written as · ¸ 1 1 p(~ z | C) = p exp − Q(~ z) , 2 |2πC| where · ¸ ~ α−~ α0 ~ z≡ , ~ r −~ rm · ¸ Cα Cαr C≡ T
Cαr
Cm
,
z , Q(~ z) ≡~ z T C−1~ and ~ α ≡ model variables, ~ α 0 ≡ prior model values,
Cα ≡ ~ α covariances, ~ r ≡ measurement variables, ~ r m ≡ measurement data,
Cm ≡~r m covariances, and Cαr ≡ ~ α-to-~ r m “covariances.” The ~ z with the maximum posterior probability is obtained by minimizing Q(~ z). Changes ~ α) are computed using a linearization in the response ~ r = R(~ ¯ ∂R i ¯¯ ~ ~ r = R(~ αk ) + S(~ αk )(~ α−~ αk ) where [S(~ αk )]i j = . ∂α j ¯~αk For a linear system, the sensitivities S do not depend on the state variable ~ α, and Q(~ z) can be minimized in one step. For a nonlinear system, S is a function of ~ α, and the minimization requires more work. DA Nonlinear Iteration One common approach for minimizing nonlinear functionals is Newton’s method, where the functional’s Hessian ∇2αQ and gradient ~ ∇αQ are combined to compute an updated step: ¡ ¢−1 ~ αk+1 = ~ αk − λk ∇2αQ(~ zk ) ~ ∇αQ(~ zk ) . (2) λk ∈ (0, 1] adjusts the step length. The initial iterate ~ α0 is set to the mean of the prior ~ α 0. At each iteration ~ αk , the sensitivities are recomputed. In the context of nonlinear least 3
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
~ squares problems, the typical approach is to ignore the second-order derivatives of R, resulting in the Gauss-Newton method [6, §2.4]. This approach is called a quasi-Newton α ∗ , the covariance matrix C can be method by Tarantola [4, §3.4.4]. At the optimal point ~ computed with Cacuci’s formulas, provided the sensitivities are evaluated at ~ α ∗. Change of Variables One final algorithmic detail is a change-of-variables for the state variables. Instead of letting αi = ∆ρ i , we let αi = loge (∆ρ i /∆ρ 0 ), where ∆ρ 0 is some reference thickness. This has two advantages. First, as noted by Tarantola [4, §1.2.4], using the logarithmic thicknesses does not favor 1 cm over 10 cm or 100 cm. This decreases the information carried by the prior P D F , a desirable result when little is initially known of the solution to an inverse transport problem. Second, the logarithmic transform disallows unphysical negative thicknesses [4, §3.2.1]. For a more equal comparison, we also performed this transformation for the L M method.
RESULTS We used the S C A L E code X S D R N P M to solve the forward transport problem and calculate group-wise leakage current on the problem’s outer boundary. For the idealized “detector” response, we used the leakage at the boundary. We also used X S D R N P M to compute a “measured” detector response for the nominal solution. Group-wise uncertainties were assigned to the leakage current based on estimated Poisson counting statistics in a detector at 1 meter from the center of the sphere. The computed “measurements” were then randomly perturbed according to a normal distribution with standard deviations from the Poisson statistics. Thus, there are no model-to-experiment inconsistencies or systematic errors. We intend to characterize the effect of systematic modeling errors in the near future. Finite differences computed the Jacobians and sensitivities. We then applied the two methods to Favorite’s test problem; the results are summarized in Figure 2. The mean of the prior and the initial iterate for both methods was ∆ρ i = e ≈ 2.718 cm for all shells. The width of the prior P D F for D A and the analogous weights w j for L M were adjusted to give equal (lack of) confidence in the prior estimate for the two methods.
CONCLUSIONS The L M and D A methods look quite similar when considering the minimization procedure, especially when the weighted χ2 cost functional is used in L M . However, the D A approach is based on a rigorous treatment of uncertainty, which allows for a firmer interpretation of the resulting P D F s. In the test problem, the D A estimates had more accurate means and smaller standard deviations for the shell thickness solution.
4
Appearing in Transactions of the American Nuclear Society, Winter 2011
Prior
preprint
DA & LM
0
10
20
30
40
50
60
DA Uranium LM 8.5
8.6
8.7
8.8
8.9
9.0
9.1
9.2
9.3
DA Void LM 0
5
10
15
20
25
DA Lead LM 0.45
0.50
0.55
0.60
0.65
0.70
0.75
DA Aluminum LM 0.26
0.28
0.30
0.32
0.34
0.36
0.38
0.40
Thicknesses Δρ (cm)
Figure 2: Shell thickness estimates with accompanying 1σ uncertainty for data assimilation ( D A ) and Levenberg-Marquardt ( L M ). The prior estimate, shown in the top subplot, was the same for each of the layers. The green vertical line indicates the true value.
5
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
References [1] J. A. FAVORITE, “Using the Schwinger variational functional for the solution of inverse transport problems,” Nuclear Science and Engineering, 146, 1, 51–70 (2004). [2] J. A. FAVORITE and D. I. KETCHESON, “Using the Levenberg-Marquardt method for the solution of inverse transport problems,” Transactions of the American Nuclear Society, 95, 527–531 (2006). [3] J. MATTINGLY and D. MITCHELL, “A framework for the solution of inverse radiation transport problems,” IEEE Transactions on Nuclear Science, 57, 6, 3734–3743 (2010). [4] A. TARANTOLA, Inverse problem theory and methods for model parameter estimation, Society for Industrial and Applied Mathematics (2005). [5] D. CACUCI and M. IONESCU-BUJOR, “Best-estimate model calibration and prediction through experimental data assimilation—I: Mathematical framework,” Nuclear Science and Engineering, 165, 1, 18–44 (2010). [6] C. T. KELLEY, Iterative methods for optimization, vol. 18, Society for Industrial and Applied Mathematics (1999).
6
INTRODUCTION Implicit methods for inverse radiation transport problems fit a transport model to an experiment to estimate properties of the source measured in the experiment. Several solution algorithms have been applied, and others are still being developed. One successful solution method is nonlinear regression, where typically the Levenberg-Marquardt ( L M ) algorithm is used to minimize a chi-squared metric of the difference between the model prediction and the experiment. In this work, we apply a data assimilation ( D A ) method to an inverse transport problem; the main advantage of D A is that it is designed to accept uncertainty in the inputs and to provide uncertainty estimates for the results. We compare L M to D A using a shielded uranium sphere as a simple test problem, first published by Favorite [1]. In Ref. [2], Favorite and Ketcheson demonstrated the use of L M for this problem. (In contrast to the current work, they used adjoint solutions to compute derivatives, and they only tracked the uncollided spectrum.) As shown in Figure 1, the uranium sphere is surrounded by concentric spherical shells of void, lead, and aluminum. Following Favorite’s example, we neglect neutron emissions by the uranium and analyze only its gamma emissions. The gammas are counted by an idealized “detector” beyond the problem’s outer boundary. The goal is to find the shell thicknesses to optimally reproduce the gamma spectrum. The unknowns we solved for are the four layer thicknesses, ~ = [8.741, 3.659, 0.5, 0.3] (all in cm). whose true values are ∆ρ
1
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
Methods Nonlinear Least Squares Procedure For comparison purposes, we implemented an L M inverse solver in a similar manner as Mattingly and Mitchell [3]. We minimize the χ2 functional ~ = χ2 (∆ρ)
£ ¤ ~ −le 2 m l c (∆ρ) X i i σ2i
i =1
+
³ ´2 0 n ∆ρ j − ∆ρ X j w 2j
j =1
,
(1)
where m = number of leakage energy groups (21 here), n = number of thickness unknowns (4 here), l ic = computed leakage in group i , l ie = experimentally measured leakage in group i , σ2i = experimental variance of l ie , ∆ρ j = j th shell thickness unknown, ∆ρ 0j = j th shell thickness prior estimate, and w j = uncertainty weighting of j th prior estimate. With this definition of χ2 , the covariance of the final parameters can be estimated as (JT J)−1 , where the Jacobian is defined as Ji j =
c 1 ∂l i σi ∂∆ρ j
.
Data Assimilation Procedure Probability theory and Bayes’ theorem provide a general approach to solve inverse problems (for example, see Tarantola [4]). Although the most general description of an inverse solution is a probability distribution function over the unknowns, this full description is Aluminum Lead Uranium Void
8.741 12.9 13.2
Figure 1: A cross section of Favorite’s shielded uranium sphere. All dimensions are in cm.
2
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
usually too complex to compute or easily comprehend. Thus, a typical simplification is to track the means and covariances of these distributions. One example of this approach is the data assimilation methods presented by Cacuci and Ionescu-Bujor [5]. Their posterior distribution is a multivariate Gaussian, written as · ¸ 1 1 p(~ z | C) = p exp − Q(~ z) , 2 |2πC| where · ¸ ~ α−~ α0 ~ z≡ , ~ r −~ rm · ¸ Cα Cαr C≡ T
Cαr
Cm
,
z , Q(~ z) ≡~ z T C−1~ and ~ α ≡ model variables, ~ α 0 ≡ prior model values,
Cα ≡ ~ α covariances, ~ r ≡ measurement variables, ~ r m ≡ measurement data,
Cm ≡~r m covariances, and Cαr ≡ ~ α-to-~ r m “covariances.” The ~ z with the maximum posterior probability is obtained by minimizing Q(~ z). Changes ~ α) are computed using a linearization in the response ~ r = R(~ ¯ ∂R i ¯¯ ~ ~ r = R(~ αk ) + S(~ αk )(~ α−~ αk ) where [S(~ αk )]i j = . ∂α j ¯~αk For a linear system, the sensitivities S do not depend on the state variable ~ α, and Q(~ z) can be minimized in one step. For a nonlinear system, S is a function of ~ α, and the minimization requires more work. DA Nonlinear Iteration One common approach for minimizing nonlinear functionals is Newton’s method, where the functional’s Hessian ∇2αQ and gradient ~ ∇αQ are combined to compute an updated step: ¡ ¢−1 ~ αk+1 = ~ αk − λk ∇2αQ(~ zk ) ~ ∇αQ(~ zk ) . (2) λk ∈ (0, 1] adjusts the step length. The initial iterate ~ α0 is set to the mean of the prior ~ α 0. At each iteration ~ αk , the sensitivities are recomputed. In the context of nonlinear least 3
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
~ squares problems, the typical approach is to ignore the second-order derivatives of R, resulting in the Gauss-Newton method [6, §2.4]. This approach is called a quasi-Newton α ∗ , the covariance matrix C can be method by Tarantola [4, §3.4.4]. At the optimal point ~ computed with Cacuci’s formulas, provided the sensitivities are evaluated at ~ α ∗. Change of Variables One final algorithmic detail is a change-of-variables for the state variables. Instead of letting αi = ∆ρ i , we let αi = loge (∆ρ i /∆ρ 0 ), where ∆ρ 0 is some reference thickness. This has two advantages. First, as noted by Tarantola [4, §1.2.4], using the logarithmic thicknesses does not favor 1 cm over 10 cm or 100 cm. This decreases the information carried by the prior P D F , a desirable result when little is initially known of the solution to an inverse transport problem. Second, the logarithmic transform disallows unphysical negative thicknesses [4, §3.2.1]. For a more equal comparison, we also performed this transformation for the L M method.
RESULTS We used the S C A L E code X S D R N P M to solve the forward transport problem and calculate group-wise leakage current on the problem’s outer boundary. For the idealized “detector” response, we used the leakage at the boundary. We also used X S D R N P M to compute a “measured” detector response for the nominal solution. Group-wise uncertainties were assigned to the leakage current based on estimated Poisson counting statistics in a detector at 1 meter from the center of the sphere. The computed “measurements” were then randomly perturbed according to a normal distribution with standard deviations from the Poisson statistics. Thus, there are no model-to-experiment inconsistencies or systematic errors. We intend to characterize the effect of systematic modeling errors in the near future. Finite differences computed the Jacobians and sensitivities. We then applied the two methods to Favorite’s test problem; the results are summarized in Figure 2. The mean of the prior and the initial iterate for both methods was ∆ρ i = e ≈ 2.718 cm for all shells. The width of the prior P D F for D A and the analogous weights w j for L M were adjusted to give equal (lack of) confidence in the prior estimate for the two methods.
CONCLUSIONS The L M and D A methods look quite similar when considering the minimization procedure, especially when the weighted χ2 cost functional is used in L M . However, the D A approach is based on a rigorous treatment of uncertainty, which allows for a firmer interpretation of the resulting P D F s. In the test problem, the D A estimates had more accurate means and smaller standard deviations for the shell thickness solution.
4
Appearing in Transactions of the American Nuclear Society, Winter 2011
Prior
preprint
DA & LM
0
10
20
30
40
50
60
DA Uranium LM 8.5
8.6
8.7
8.8
8.9
9.0
9.1
9.2
9.3
DA Void LM 0
5
10
15
20
25
DA Lead LM 0.45
0.50
0.55
0.60
0.65
0.70
0.75
DA Aluminum LM 0.26
0.28
0.30
0.32
0.34
0.36
0.38
0.40
Thicknesses Δρ (cm)
Figure 2: Shell thickness estimates with accompanying 1σ uncertainty for data assimilation ( D A ) and Levenberg-Marquardt ( L M ). The prior estimate, shown in the top subplot, was the same for each of the layers. The green vertical line indicates the true value.
5
Appearing in Transactions of the American Nuclear Society, Winter 2011
preprint
References [1] J. A. FAVORITE, “Using the Schwinger variational functional for the solution of inverse transport problems,” Nuclear Science and Engineering, 146, 1, 51–70 (2004). [2] J. A. FAVORITE and D. I. KETCHESON, “Using the Levenberg-Marquardt method for the solution of inverse transport problems,” Transactions of the American Nuclear Society, 95, 527–531 (2006). [3] J. MATTINGLY and D. MITCHELL, “A framework for the solution of inverse radiation transport problems,” IEEE Transactions on Nuclear Science, 57, 6, 3734–3743 (2010). [4] A. TARANTOLA, Inverse problem theory and methods for model parameter estimation, Society for Industrial and Applied Mathematics (2005). [5] D. CACUCI and M. IONESCU-BUJOR, “Best-estimate model calibration and prediction through experimental data assimilation—I: Mathematical framework,” Nuclear Science and Engineering, 165, 1, 18–44 (2010). [6] C. T. KELLEY, Iterative methods for optimization, vol. 18, Society for Industrial and Applied Mathematics (1999).
6
