Source-linked AI summary
Monte Carlo Power Iteration Algorithms: Convergence, Analysis and Application to Reactor Criticality Calculations
Alexander M. G. Cox, Tom Davis, Julian Hofstadler, Andreas E. Kyprianou
TL;DR
Monte Carlo power iteration is widely used to estimate reactor criticality, but rigorous error bounds for population-controlled methods were missing to the authors’ knowledge. The paper analyzes these algorithms through Sequential Monte Carlo under suitable physical-system and algorithmic assumptions. It proves explicit error bounds, including for a standard variant and variance-reduction schemes.
Problem
Monte Carlo power iteration is used for reactor criticality calculations, but the authors report that rigorous error analysis with explicit assumptions and bounds was missing to their knowledge.
Source: §1 · pp. 2–3Method
The authors treat population-controlled power iteration as Sequential Monte Carlo and derive bounds under assumptions on the physical system and algorithmic structure.
Source: §1 · pp. 2–3Results
The paper proves explicit error bounds with an N^-1/2 Monte Carlo term, a burn-in and iteration term, and an extinction term, and gives a simplified bound for surviving genealogies.
Source: §1 · p. 3Takeaways & Limitations
The authors state that the required conditions cover a large class of applications, including a standard power-iteration variant and variance-reduction techniques.
Source: §1 · p. 3A potentially sharper (kN)^-1/2 bound for arithmetic-mean estimators is left for future investigation.
Source: §5.2.2 · p. 19
Abstract
from arXiv · showhide
We study Monte Carlo power iteration algorithms in the context of reactor criticality simulations. Under mild conditions on the underlying physical system and the algorithmic structure, we derive explicit error bounds, which are to the best of our knowledge the first rigorous results for such methods which are ubiquitous in nuclear engineering. This is done by exploiting connections to Sequential Monte Carlo algorithms, and extending recent developments about stochastic neutron transport. We demonstrate the applicability of our results by considering a concrete version of the Monte Carlo power iteration algorithm as well as variance reduction schemes.
1 Introduction
Monte Carlo power iteration estimates the reactor criticality eigenvalue, but population counts can grow or shrink exponentially away from criticality. The paper provides rigorous error bounds by treating the algorithm as Sequential Monte Carlo and applying general SMC bounds under conditions stated in terms of the physical system and population-control step.
1 Introduction
- The keff eigenvalue characterizes neutron gain or loss across fission generations, while population control fixes each generation at N particles to avoid exponential population change away from criticality.The algorithm evolves particles until fission or exit, then selects exactly N offspring; keff is estimated from the controlled population.Source: §1 · pp. 1–2
- The paper addresses the stated lack of rigorous error analysis for a standard nuclear-engineering algorithm, whose estimates are used in criticality calculations and popular simulation tools.The authors describe earlier work as providing practical understanding, while noting that precisely formulated assumptions and error bounds were missing to their knowledge.Source: §1 · pp. 2–4
- The error bound has a typical N^-1/2 Monte Carlo term, a burn-in and iteration term, and a term for the possibility that all particles die at an iteration.The burn-in term decreases with a longer burn-in or more generations used in the estimate; the extinction contribution reflects a probability bounded by e^-cN.Source: §1 · p. 3
- A separate bound for estimates based on surviving genealogies is also proved.The supplied passage says this estimate simplifies but does not give its complete expression.Source: §1 · p. 3
- The analysis treats population-controlled Monte Carlo power iteration as a Sequential Monte Carlo method and first establishes general SMC bounds before applying them to power iteration.The underlying process is represented by an operator whose potential gives expected offspring and whose Markov operator describes neutron behavior between fission events.Source: §1 · pp. 3–4
- The authors state that their assumptions are met for a large class of applications and can be expressed through neutron-transport cross sections and population-control requirements.They cite a standard power-iteration variant and variance-reduction techniques as examples meeting the algorithmic requirements.Source: §1 · p. 3
2 Stochastic Neutron Transport
The section connects neutron transport to probability through branching processes and generational particle dynamics, then states physical conditions used for the analysis.
2.1 The probabilistic model
- The transport rates and densities induce a branching Markov process whose finite counting measures track live neutrons and their configurations.The rates and densities are σs, σf, πs, and πf; the state records the live-particle count and position–velocity configurations.Source: §2.1 · p. 5
- The continuous-time process yields a mean semigroup that solves the transport evolution equation.Its mean semigroup is defined by expected particle sums, with the transport operators specified in the section.Source: §2.1 · p. 5
- The neutron generational process follows lineages across fission events, providing a semigroup that connects the transport model to sequential Monte Carlo.The section defines generation-level measures and their expectation semigroup before introducing physical-model assumptions.Source: §2.1 · pp. 5–6
2.2 Assumptions on the physical model
- The reactor model separates fissile and non-fissile regions, with fission confined to the fissile region and that region separated from the reactor boundary.The physical assumptions also require a bounded, open, convex domain and a positive-radius ball within it.Source: §2.2 · p. 6
- The physical assumptions also impose lower bounds on scattering-related quantities and upper bounds on cross-sections and fission offspring production.These conditions are presented as assumptions on the scattering events, rates, and maximum particles released at a fission event.Source: §2.2 · p. 6
3 Mathematical setting
This section sets up an abstract positive-operator eigenvalue problem and the conditions connecting it to neutron transport. It then establishes convergence properties that support later estimators.
3.1 The abstract keff setting
- The abstract setting uses a bounded nonnegative operator Q, with a leading eigentriple whose eigenvalue dominates the moduli of other eigenvalues.In nuclear applications, Q corresponds to the expectation semigroup of the neutron generational process.Source: §3.1 · p. 7
- Factoring Q into a potential G and Markov kernel M accommodates states with zero expected offspring, while requiring positive probability of transition out of nonzero-potential states.Under the eigenvalue assumption, the eigenmeasure satisfies η(G) = k_eff.Source: §3.1 · pp. 7–8
3.2 Interacting particle systems notation
- Under the physical assumptions, the leading eigenquantities and exponential convergence assumptions are verified even when fission is restricted to a fissile region.The section adapts earlier arguments to the setting with non-fissile regions and proves the required assumptions under S, F, and B.Source: §3.2 · p. 9
- The normalized measure flow converges in total variation at a geometric rate controlled by τ.For x outside E0, the stated bound is ||Φ_n(δ_x) − η||_tv ≤ eCτ^n.Source: §3.2 · pp. 9–10
- Several population-flow quantities converge to k_eff, including η_n(G), which motivates arithmetic- and geometric-mean estimators.The section states η0Q^n(1) = η_n(G) = Φ_n(η0)(G), with η_n(G) converging to k_eff.Source: §3.2 · p. 10
4 Sequential Monte Carlo Algorithms
The section presents a sequential Monte Carlo scheme for neutron transport and relates its particle flow to power iteration. It then identifies arithmetic and geometric means as estimators of k_eff.
4 Sequential Monte Carlo Algorithms
- Sequential Monte Carlo algorithms generate samples from target distributions and provide the framework used here to analyze Monte Carlo power iteration.The paper reviews the required SMC ingredients before applying them to power iteration.Source: §4 · p. 10
4.1 SMC and power iteration for neutron transport
- The bootstrap particle filter evolves N weighted particles using resampling and particle evolution, with G representing expected offspring and M the transition to the next fission event.The scheme is connected to neutron transport through the potential and Markov kernel.Source: §4.1 · p. 11
4.2 Estimators for keff
- Arithmetic or geometric means of the particle-based estimates provide the stated approaches for estimating k_eff.The section motivates these estimators using limiting relations for the underlying population quantities.Source: §4.2 · p. 11
5 Lp error bounds
This section develops Lp error bounds for particle-measure and k_eff estimators, covering extinction and arithmetic and geometric means under stated assumptions on survival, mixing, and particle-measure stability. For geometric-mean estimators, it also establishes that the expected value is bounded away from zero.
5.1 Extinction probabilities
- If each live state has transition probability and weight at least κ, the probability of extinction from one iteration to the next is at most (1 − κ)^N.Assumption 5.1 requires M(x, E \ E0) ≥ κ and G(x) ≥ κ for x outside E0.Source: §5.1 · pp. 12–13
5.2.1 Stability of the particle measures ηN
- Under Assumptions 3.4, 5.1, and 5.3, Theorem 5.4 bounds particle-measure estimation error for bounded measurable functions across generations, and the measures also support fluence estimation by tallying.The theorem gives a constant K_p independent of t, and the authors describe tallying as the use of particle measures to compute fluence estimates.Source: §5.2.1 · pp. 13, 17
5.2.2 Arithmetic mean estimators
- Under Assumptions 3.4, 5.1, and 5.3 for some p ∈ [1, ∞), Theorem 5.9 asserts a constant Kp ∈ (0, ∞) exists for any T, k, N ∈ N.It defines a = −log(1 − κ) ∈ (0, ∞), where κ is as in Assumption 5.1.Source: §5.2.2 · p. 17
- For any ε, δ > 0, Theorem 5.11 states that sufficiently large simulation parameters can achieve the specified error tolerance and failure probability.The result is framed around choosing T, k, and N while accounting for extinction before T + k and estimator error conditional on termination.Source: §5.2.2 · p. 18
- The cost analysis gives error scaling like ε^-2 ln(ε^-1) for small ε; when the spectral-gap parameter β is small, cost grows like β^-1 and the optimal averaging length increases.The stated cost model makes total numerical cost proportional to N(T + k).Source: §5.2.2 · p. 19
5.2.3 Geometric mean estimators
- For geometric-mean estimators, Lemma 5.14 provides a q-dependent bound under Assumptions 3.4 and 5.1 for any generation t and particle count N.The supplied passage does not expose the bound’s formula.Source: §5.2.3 · p. 20
- Under Assumptions 3.4 and 5.1, the infimum over n ∈ N of ηn(G) is strictly positive.The proof uses convergence of ηn(G) to keff > 0.Source: §5.2.3 · p. 20
- Under Assumptions 3.4, 5.1, and 5.3 for some p′ > 1, Theorem 5.16 gives a bound with constant Kp,p′ for p ∈ [1, p′) and any T, k, N ∈ N.The supplied passage does not include the bound’s formula.Source: §5.2.3 · p. 20
- Corollary 5.17 gives a corresponding bound with constant Kp,p′ for p ∈ [1, p′) under Assumptions 3.4, 5.1, and 5.3.The supplied passage does not include the bound’s formula.Source: §5.2.3 · p. 23
- With essentially identical arguments and conclusions, the result can be modified to give geometric-estimator versions of Theorem 5.11 and Remark 5.13.Source: §5.2.3 · p. 23
6 Power Iteration Algorithms
The section applies the paper’s error analysis to explicit fission power iteration and to variance-reduction variants. It derives bounds for the explicit-fission algorithm and for variants whose particle estimators meet independence and unbiasedness conditions.
6.1 Explicit Fission
- Algorithm 1 initializes a particle population, resamples particles, evolves them, and generates a new bank through offspring production over T burn-in and k estimation generations.The algorithm stops with an error if the population dies out.Source: §6.1 · p. 23
6 Power Iteration Algorithms
- The section first frames explicit fission as a standard resampling-and-evolution procedure, then extends its analysis to variance-reduction changes.A common variation adds offspring after fission rather than copying the particle about to undergo fission; the paper treats it as Algorithm 1 up to one extra iteration.Source: §6 · p. 23, §6.1 · p. 24
6.1 Explicit Fission
- For explicit fission, Proposition 6.3 verifies the required condition for every p in [1, ∞), allowing the Section 5 results to yield error bounds for the arithmetic-mean estimator.The bounds are stated for systems fitting the paper’s Section 2.2 setting.Source: §6.1 · pp. 23–25
6.2 Variance reduction techniques
- Variance-reduction variants inherit error bounds when each particle’s estimator is conditionally independent and unbiased; the paper notes that broader variants can qualify if they preserve both properties.The modifications typically change evolution or resampling and are described as helping spread out particles.Source: §6.2 · pp. 25–26