Data assimilation applied to a simple inverse radiation transport problem
Data assimilation applied to a simple inverse radiation transport problem
Joshua M Hykes (Jun Li – UNC)
Dan G Cacuci
John Mattingly
Department of Nuclear Engineering North Carolina State University 2011 Winter ANS Meeting November 3, 2011
The goal is to characterize special nuclear materials without opening the container.
? Radiation detector
Unknown materials À 2 23
The available data is a measured gamma spectrum
Photopeaks Counts
Full spectrum Energy À 3 23
With no other information, 1-D spherical geometry is the simplest non-point model Aluminum Lead Uranium Void
8.741 12.9 13.2 [as suggested by Favorite in 2004] À 4 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
À 5 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
Model Model Parameters
À 6 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
Model Model Parameters
À 7 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
Model Model Parameters
À 8 23
Nonlinear least-squares methods include ~ errorbars on ~rc and (sometimes) α
~ = χ2 (α)
2 m (R (α) X i ~ − rm,i ) i=1
σr2m ,i
+
2 n ~j − α ~ j0 α X j=1
2 σα,j
,
−1 ~ ∗ )T J(α ~ ∗) Cα = J(α
À 9 23
The Bayesian posterior distribution includes all of the errorbars (and correlations) p(~ z | C) = p
1 |2πC|
1
exp − Q(~ z) 2
where ~ ≡ model variables α ~ 0 ≡ prior model values α ~ covariances Cα ≡ α ~r ≡ measurement variables ~rm ≡ measurement data Cm ≡ ~rm covariances ~ Cαr ≡ α-to-~ rm “covariances”
~ −α ~0 α ~≡ z ~r − ~rm
Cα Cαr C≡ CTαr Cm
~T C−1 z ~ Q(~ z) ≡ z [Cacuci 2010] À 10 23
Gauss-Newton optimization of Q(~ z) searches for ~ ∗ with the maximum posterior α −1 ~ α Q(~ ~ m+1 = α ~ m − λm ∇2α Q(~ α zm ) ∇ zm )
with the linearization
~ (α ~r = R ~ m ) + S(α ~ m )(α ~ −α ~ m)
where
~ i ∂R [S(α)] ~ i,j = ∂αj
~ α
[Kelley 1999]
À 11 23
~ The thicknesses of each material are put in α,
loge ~ = α
~ ∆ρ ∆ρ0
À 12 23
~ The thicknesses of each material are put in α, as well as uncertain material properties
~ ∆ρ log e ∆ρ0 ~ = α σ~t
À 13 23
We solved Favorite’s test problem 2 different ways.
1
Full spectrum, 22 energy groups • “Measurement”: SN solver • Model: SN solver
Aluminum Lead
2
Uncollided photopeaks, 4 peaks • “Measurement”: MCNP • Model: Ray tracer (written by Jun
Li)
Uranium Void
8.741 12.9 13.2
À 14 23
With no modeling discrepancy and an overdetermined solution, the multigroup full spectrum results are accurate (except void) Prior
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
Thicknesses Δρ (cm)
0.36
0.38
0.40
À 15 23
The uncollided photopeak results have (much) larger errorbars Prior 0
DA & LM 10
20
30
40
50
60
DA Uranium LM 8.5
8.0
9.0
9.5
10.0
DA Void 0
5
10
LM 15
20
25
30
DA Lead 0.44
0.46
0.48
0.50
0.52
LM 0.54
0.56
DA Aluminum LM 0.0 0.5
1.0
1.5 2.0 2.5 3.0 Thicknesses Δρ (cm)
3.5
4.0
4.5 À 16 23
The uncollided photopeak results have (much) larger errorbars ~σΔρ / Δρ Prior 0
10
20
30
40
50
60
DA Uranium LM 8.5
8.0
9.0
9.5
Void 0
5
10
LM 15
20
25
DA
0.48
0.50
0.52
LM 0.54
DA
1.0
1.5 2.0 2.5 3.0 Thicknesses Δρ (cm)
18
274%
0.12
11%
12
9.3%
0.088
1.6%
4.9
323%
0.43
1002%
2.4
0.56
Aluminum LM 0.0 0.5
0.22
0.24%
30
Lead 0.46
11%
10.0
DA
0.44
error / σΔρ
DA & LM
3.5
4.0
4.5 À 17 23
Conclusions
• More (good) information gives tighter errorbars on
the solution. • Errorbars on the least squares LM results are
qualitatively helpful. • Errorbars on the data assimilation results are
reasonably accurate. • Data assimilation provides a useful norm in which to
gauge the measurement-to-model misfit.
À 18 23
References
• Albert Tarantola, Inverse problem theory and methods for model parameter estimation, SIAM 2005. • Jeffrey Favorite, “Using the Schwinger variational functional for the solution of inverse transport problems” Nuclear Science and Engineering 146 2004. • Dan Cacuci and Mihaela Ionescu-Bujor, “Best-estimate model calibration and prediction through experimental data assimilation–I: mathematical framework” Nuclear Science and Engineering 165 2010. • CT Kelley, Iterative Methods for Optimization, SIAM 1999.
À 19 23
Extra Slides
À 20 23
Bayesian inference is a general approach to solving inverse problems p(model | data) ∝ p(data | model) · p(model) likelihood
x
×
p(x)
p(x)
=
x
prior
p(x)
posterior
x
See Inverse Problem Theory by Albert Tarantola for more details. À 21 23
Explicit formulas for gradient and Hessian of Q Partition the inverse covariance matrix: C1 C3 −1 C = CT3 C2
~αT C1 z ~α + 2~ ~r + z ~rT C2 z ~r Q(~ z) = z zαT C3 z ∂Q ~ ~α + ST C2 z ~r + C3 z ~r + ST CT3 z ~α ∇α Q(~ z) = = C1 z ~ ∂α ~ z 2Q ∂ ∂S ∂S T T T T ~α ~ ∇2α Q(~ z) = = C + S C S + C S + S C + C C z z + 1 2 3 2 r 3 ~2 ~ ~ 3 ∂α ∂ α ∂ α ~ z ∂S ~ ∂α
terms are ignored in Gauss-Newton method. ~ ∗ , the best-estimate covariance At the optimal point α be matrix C can be computed with Cacuci’s formulas, ~ ∗. provided the sensitivities are evaluated at α À 22 23
Some constraints or penalties are required to contain the early iterates During the initial iterations, the solution tends to bounce around. Sometimes a shell thickness for one of these iterates is very large, such that the transport solver breaks. One option is to use a bound-constrained optimization method, where an upper and lower bound can be ~ Although this seems supplied for each component in α. like a more elegant solution, the numerical experiments with this method were less stable than the next option. The approach we took was essentially a penalty method. If the shell thicknesses caused the transport solver to break, then zero leakages or fluxes were returned. This was sufficient penalty to return the next iterate to a more plausible region. À 23 23
Joshua M Hykes (Jun Li – UNC)
Dan G Cacuci
John Mattingly
Department of Nuclear Engineering North Carolina State University 2011 Winter ANS Meeting November 3, 2011
The goal is to characterize special nuclear materials without opening the container.
? Radiation detector
Unknown materials À 2 23
The available data is a measured gamma spectrum
Photopeaks Counts
Full spectrum Energy À 3 23
With no other information, 1-D spherical geometry is the simplest non-point model Aluminum Lead Uranium Void
8.741 12.9 13.2 [as suggested by Favorite in 2004] À 4 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
À 5 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
Model Model Parameters
À 6 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
Model Model Parameters
À 7 23
~ to “best” Inverse methods seek model values α fit computed and measured responses
Reality
?
Model Model Parameters
À 8 23
Nonlinear least-squares methods include ~ errorbars on ~rc and (sometimes) α
~ = χ2 (α)
2 m (R (α) X i ~ − rm,i ) i=1
σr2m ,i
+
2 n ~j − α ~ j0 α X j=1
2 σα,j
,
−1 ~ ∗ )T J(α ~ ∗) Cα = J(α
À 9 23
The Bayesian posterior distribution includes all of the errorbars (and correlations) p(~ z | C) = p
1 |2πC|
1
exp − Q(~ z) 2
where ~ ≡ model variables α ~ 0 ≡ prior model values α ~ covariances Cα ≡ α ~r ≡ measurement variables ~rm ≡ measurement data Cm ≡ ~rm covariances ~ Cαr ≡ α-to-~ rm “covariances”
~ −α ~0 α ~≡ z ~r − ~rm
Cα Cαr C≡ CTαr Cm
~T C−1 z ~ Q(~ z) ≡ z [Cacuci 2010] À 10 23
Gauss-Newton optimization of Q(~ z) searches for ~ ∗ with the maximum posterior α −1 ~ α Q(~ ~ m+1 = α ~ m − λm ∇2α Q(~ α zm ) ∇ zm )
with the linearization
~ (α ~r = R ~ m ) + S(α ~ m )(α ~ −α ~ m)
where
~ i ∂R [S(α)] ~ i,j = ∂αj
~ α
[Kelley 1999]
À 11 23
~ The thicknesses of each material are put in α,
loge ~ = α
~ ∆ρ ∆ρ0
À 12 23
~ The thicknesses of each material are put in α, as well as uncertain material properties
~ ∆ρ log e ∆ρ0 ~ = α σ~t
À 13 23
We solved Favorite’s test problem 2 different ways.
1
Full spectrum, 22 energy groups • “Measurement”: SN solver • Model: SN solver
Aluminum Lead
2
Uncollided photopeaks, 4 peaks • “Measurement”: MCNP • Model: Ray tracer (written by Jun
Li)
Uranium Void
8.741 12.9 13.2
À 14 23
With no modeling discrepancy and an overdetermined solution, the multigroup full spectrum results are accurate (except void) Prior
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
Thicknesses Δρ (cm)
0.36
0.38
0.40
À 15 23
The uncollided photopeak results have (much) larger errorbars Prior 0
DA & LM 10
20
30
40
50
60
DA Uranium LM 8.5
8.0
9.0
9.5
10.0
DA Void 0
5
10
LM 15
20
25
30
DA Lead 0.44
0.46
0.48
0.50
0.52
LM 0.54
0.56
DA Aluminum LM 0.0 0.5
1.0
1.5 2.0 2.5 3.0 Thicknesses Δρ (cm)
3.5
4.0
4.5 À 16 23
The uncollided photopeak results have (much) larger errorbars ~σΔρ / Δρ Prior 0
10
20
30
40
50
60
DA Uranium LM 8.5
8.0
9.0
9.5
Void 0
5
10
LM 15
20
25
DA
0.48
0.50
0.52
LM 0.54
DA
1.0
1.5 2.0 2.5 3.0 Thicknesses Δρ (cm)
18
274%
0.12
11%
12
9.3%
0.088
1.6%
4.9
323%
0.43
1002%
2.4
0.56
Aluminum LM 0.0 0.5
0.22
0.24%
30
Lead 0.46
11%
10.0
DA
0.44
error / σΔρ
DA & LM
3.5
4.0
4.5 À 17 23
Conclusions
• More (good) information gives tighter errorbars on
the solution. • Errorbars on the least squares LM results are
qualitatively helpful. • Errorbars on the data assimilation results are
reasonably accurate. • Data assimilation provides a useful norm in which to
gauge the measurement-to-model misfit.
À 18 23
References
• Albert Tarantola, Inverse problem theory and methods for model parameter estimation, SIAM 2005. • Jeffrey Favorite, “Using the Schwinger variational functional for the solution of inverse transport problems” Nuclear Science and Engineering 146 2004. • Dan Cacuci and Mihaela Ionescu-Bujor, “Best-estimate model calibration and prediction through experimental data assimilation–I: mathematical framework” Nuclear Science and Engineering 165 2010. • CT Kelley, Iterative Methods for Optimization, SIAM 1999.
À 19 23
Extra Slides
À 20 23
Bayesian inference is a general approach to solving inverse problems p(model | data) ∝ p(data | model) · p(model) likelihood
x
×
p(x)
p(x)
=
x
prior
p(x)
posterior
x
See Inverse Problem Theory by Albert Tarantola for more details. À 21 23
Explicit formulas for gradient and Hessian of Q Partition the inverse covariance matrix: C1 C3 −1 C = CT3 C2
~αT C1 z ~α + 2~ ~r + z ~rT C2 z ~r Q(~ z) = z zαT C3 z ∂Q ~ ~α + ST C2 z ~r + C3 z ~r + ST CT3 z ~α ∇α Q(~ z) = = C1 z ~ ∂α ~ z 2Q ∂ ∂S ∂S T T T T ~α ~ ∇2α Q(~ z) = = C + S C S + C S + S C + C C z z + 1 2 3 2 r 3 ~2 ~ ~ 3 ∂α ∂ α ∂ α ~ z ∂S ~ ∂α
terms are ignored in Gauss-Newton method. ~ ∗ , the best-estimate covariance At the optimal point α be matrix C can be computed with Cacuci’s formulas, ~ ∗. provided the sensitivities are evaluated at α À 22 23
Some constraints or penalties are required to contain the early iterates During the initial iterations, the solution tends to bounce around. Sometimes a shell thickness for one of these iterates is very large, such that the transport solver breaks. One option is to use a bound-constrained optimization method, where an upper and lower bound can be ~ Although this seems supplied for each component in α. like a more elegant solution, the numerical experiments with this method were less stable than the next option. The approach we took was essentially a penalty method. If the shell thicknesses caused the transport solver to break, then zero leakages or fluxes were returned. This was sufficient penalty to return the next iterate to a more plausible region. À 23 23
