Not Your Physicists' Entropy

Information Theory Applications to the Nuclear Fuel Cycle

Dr. Anthony Michael Scopatz
The University of Chicago, The Flash Center for Computational Science
UW-Madison, February 26th, 2013

What is the Nuclear Fuel Cycle?

Well, we can ask wikipedia...

...but that isn't very formal!

How do we model the fuel cycle?

  • Most simulations are formulated on premade base-case scenarios:

alt text

  • These base cases are very well studied.
  • However, what is not well known is how these sample scenarios are affected by perturbations to their initial physical parameters.
Do our intial parameter choices yield the ‘best’ solution?

Problem Space

Figure 1: Physics Modeled versus Execution Time

Modeling Approaches

ComponentRecipesEssentialTransport
Front End:RecipesPhysics ModelsPhysics Models
Reactor:RecipesRapid Burnup Code (Bright)Neutron Transport
Repository:Curve FitsRapid Code (Cyder)Neutron Transport

Motivation

  • Investigating the region of interest may provide a transformative amount of fuel cycle data.

  • However, this requires an entirely new spectrum of tools.

Figure 2: Fuel Cycle Stack

Definition

Essential physics models remain physically valid under perturbations in the locality of the region on which they are defined and do not compute extraneous parameters.

This is a form of model reduction.

What do Essential Physics Models Buy Us?

They allow us to run physically meaningful (if not economically, etc.) fuel cycle simulations.

And not just a few of them, but many.

And when we have many similar things we use (fancy?!) statistics to analyze the ensemble.

Essential Physics Fuel Cycle Facility Models

One-Group Reactor Model (R1G)

Figure 3: Reactor Model [1]

R1G Notation

NameSymbolUnits
Fluence$F = \int_0^t \phi dt^\prime = \phi \int_0^t dt^\prime = \phi \times t \cdot \frac{24 \cdot 3600 \cdot 10^3}{10^{24}}$[n/kb]
Nuclide Subscript$i$ (or $j$)[1/kg$_i$]
Region Superscript$Q$[unitless]
Mass Weights$m_i^Q = \frac{N_i^Q}{N_{\mbox{IHM}}} = \frac{n_i^Q A_i}{A_{\mbox{IHM}}} \cdot \frac{\rho^Q\cdot\mbox{MW}^F}{\rho^F\cdot\mbox{MW}^Q} \frac{V^Q}{V^F}$[kg$_i^Q$/kgIHM]
Burnup$\mbox{BU}(F) = \sum_i m_i^F \cdot \mbox{BU}_i(F)$[MWd/kgIHM]
Prod & Des Rates$p_i(F), d_i(F)$[n/s/flux/kg$_i$]
Transmutation Matrix$T_{ij}(F)$[kg$_j$/kg$_i$]
Multiplication Factor$k(F) = \frac{P(F)}{D(F)}$[unitless]

R1G Library

Starting with ORIGEN generated fluence-dependent, nuclide specific data, we know:

(a) neutron production & destruction rates,(b) burnups,(c) transmutation matrices.

Figure 4: Example Reactor Data Library for FR $^{239}$Pu

R1G Library Collapse

  • Collapsing a general library down is the central calculation of the one-group reactor model.
  • Using a two-region fuel pin cell model, denote the mass weight of each initial isotope $i$ in the fuel as $m^F_i$ and the in the coolant as $m^C_i$.
  • Then, a mass-weighted linear combination of all of the library parameters computes the full-core values.
  • For example,
$P(F) = P_{NL} \cdot p^F(F) = P_{NL} \sum^I_i m^F_i \cdot p_i(F)$

R1G

(a) burnup(b) multiplication factor

Figure 5: Linearly Combined Fast Reactor Data Library

R1G

Figure 6: Fast Reactor Discharge Composition

  • The transmutation matrices may be similarly weighted by initial composition.
  • From the criticality calculation above we also know the discharge fluence, Fd [n/kb], which will then yield the used fuel vector.

R1G Benchmark

Table 2: LWR Benchmark to OECD Burnup Credit [2]

NuclideOECD [2] [w/o]$\sigma$ of OECD [%]Results [w/o]% Difference
$^{234}$U0.01789.00.0184+3.6
$^{235}$U0.80018.10.7079-11.5
$^{236}$U0.48402.60.4781-1.2
$^{238}$U93.33330.293.6109+0.3
$^{237}$Np0.06149.40.0584-5.0
$^{238}$Pu0.022613.90.0185-18.1
$^{239}$Pu0.59917.10.5104+14.8
$^{240}$Pu0.23895.30.2591+8.5
$^{241}$Pu0.16366.90.1445-11.7
$^{242}$Pu0.06028.40.0563-6.5
$^{241}$Am0.00475.30.0032-30.8
$^{243}$Am0.014810.40.0116-21.3

R1G Benchmark Notes

  • This benchmark compares an OECD study for a 40 MWd/kgIHM burn with no cooling.
  • This result shows the R1G matches the burnup credit to within two standard deviations for most actinides.
  • $^{241}$Am is the only nuclide showing a relatively higher departure.
  • This species is formed through $^{241}$Pu decay.
  • $^{241}$Am would thus match better if a constant flux or if zero reloading times were not assumed.
  • More benchmark results are available here and here.

Multigroup Reactor Model (RMG)

  • In the R1G, fluxes are assumed to be static over the course of the burn.
  • This restricts the types of reactors and fuel cycles which may be modeled. (MOX, IMF, Traveling Wave, etc. are unavailable to the R1G.)
  • Also in the R1G, changing reactor state variables in a way that would affect the flux spectrum required additional libraries, or inter-library interpolations.

RMG

  • To solve these issues, moving to a multi-energy group reactor model will enable automatic spectral shifts in the core as needed.
  • Rather than parameterizing burnup, transmutation, and neutron production & destruction rates as functions of nuclide and fluence, the RMG parameterizes cross sections over a number of different states.
  • The R1G parameters are explicitly solved for by the RMG model.

RMG

  • Group constants are denoted by $\sigma_{r_x,i\tau g}$ or $\sigma_{r_x,ipg}$ [b].
  • Group-to-group scattering by $\sigma_{s,i\tau g\to h}$ or $\sigma_{s,ipg\to h}$.

Table 8: Cross Section Indices

IndexSymbol
Nuclide$i$
Time$\tau$
Perturbation$p$
Incident Energy Group$g$
Exiting Energy Group$h$
Reaction Channel [5]$r_x$

Reactor State Variables

Table 9: Parameters $a_p$ which Define a Perturbation

ParameterSymbolUnits
Fuel Density$\rho_{\mbox{fuel}}$g/cm$^{3}$
Cladding Density$\rho_{\mbox{clad}}$g/cm$^{3}$
Coolant Density$\rho_{\mbox{cool}}$g/cm$^{3}$
Fuel Cell Radius$r_{\mbox{fuel}}$cm
Void Cell Radius$r_{\mbox{void}}$cm
Cladding Cell Radius$r_{\mbox{clad}}$cm
Unit Cell Pitch$\ell$cm
Number of Burn Regions$b_r$
Fuel Specific Power$p_s$MW/kgIHM
Initial Nuclide Mass Fraction$T_{i0}$kg$_i$/kgIHM
Burn Times$s$days

Multicomponent Enrichment Model

Table 10: Outer Product Perturbations

$\rho_{\mbox{fuel}}$$\rho_{\mbox{clad}}$$\rho_{\mbox{cool}}$$r_{\mbox{fuel}}$$r_{\mbox{void}}$$r_{\mbox{clad}}$$\ell$$b_r$$p_s$$T_{\mbox{$^{235}$U0}}$$s$
10.1655.870.730.410.41850.4751.3127100.040.030
10.1655.870.730.410.41850.4751.3127100.040.03100
10.1655.870.730.410.41850.4751.3127100.040.03200
10.1655.870.730.410.41850.4751.3127100.040.03300
10.1655.870.730.410.41850.4751.3127100.040.03400
10.1655.870.730.410.41850.4751.3127100.040.050
10.1655.870.730.410.41850.4751.3127100.040.05100
10.1655.870.730.410.41850.4751.3127100.040.05200
10.1655.870.730.410.41850.4751.3127100.040.05200
10.1655.870.730.410.41850.4751.3127100.040.05300
10.1655.870.730.410.41850.4751.3127100.040.05400
$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$$\vdots$

Cross Section Generation

  • Cross section generation takes a multi-tier approach.
  • For very important nuclides, $\sigma_{pg}$ were found via the Monte Carlo neutron transport code Serpent [6] (slightly modified).
  • For important species for which high-fidelity continuous cross section data is not available, physical models and/or a flux weighted collapse of 64-group CINDER data is used.
  • For species unimportant to the FC material balance, a simple interpolation between thermal and fast cross sections from KAERI is used to generate group constants.

RMG Methodology

  • Once the cross-section library has been assembled, the RMG proceeds in three main calculation routines per time step.
  • First, there is as a nearest-neighbor calculation from the given reactor state to the perturbations states in the library.
  • A multigroup spectrum-criticality calculation.
  • A burnup-transmutation calculation.

Multigroup Reactor Model Flow Diagram

RMG Metric Space

  • With $p$ as the perturbation index ($1 \le p \le n_p$), $a_p$ is a parameter value in the library and $a_r$ the value on the reactor at run time.
  • The distance from any perturbation to the current reactor is given by $d_p$:
$d_p = \sqrt{\sum_a \left(\frac{a_r - a_p}{a_{n_p}}\right)^2}$
  • The sequence $p^\star$ determines the nearest neighbors:
$p^\star = \left\{p | d_p \le d_{p+1}\right\}$

Cross Section Interpolation

  • Using the shorthand $a_p^\star$ as a substitute for $a_{p^\star}$,
  • Then $x_f$ is a unitless multidimensional linear interpolation factor,
$x_f = \sum_a \left(\frac{a_r - a_1^\star}{a_2^\star - a_1^\star}\right)$
  • Group constants for each time step are computed:
$\sigma_{r_x,i\tau g} = (\sigma_{r_x,ip_2^\star g} - \sigma_{r_x,ip_1^\star g}) \cdot x_f + \sigma_{r_x,ip_1^\star g}$

Criticality Calculation

  • The shape of the spectrum is given by an iterative solution to a matrix method.
$A_{\tau g\to h} \cdot \phi_{\tau g} = \frac{1}{k} F_{\tau g\to h} \cdot \phi_{\tau g}$
  • ...with a rescaling to obtain the magnitude,
$\phi_{\tau} = p_s \cdot \rho_{\mbox{fuel}} \cdot \frac{1}{\mbox{3.284E-14} \cdot \Sigma_{f,\mbox{fuel},\tau}}$
  • and a side calculation to obtain the multiplication factor.
$k_{\tau} = P_{\mbox{NL}} \cdot \frac{\sum_g^G V_{\mbox{fuel}} \cdot \bar{\nu}\Sigma_{f,\mbox{fuel},\tau g} \cdot \phi_{\tau g}}{\sum_g^G \left(V_{\mbox{fuel}} \cdot \Sigma_{a,\mbox{fuel},\tau g} + \zeta_{\tau g} \cdot V_{\mbox{cool}} \cdot \Sigma_{a,\mbox{cool},\tau g}\right) \cdot \phi_{\tau g}}$

Transmutation Calculation

  • ORIGEN 2.2 [7] is used as the transmutation kernel for the RMG.
  • Custom ORIGEN cross section libraries are built for each time step. The group constants for the current reactor state are simply collapsed.
$\sigma_{r_x,i\tau} = \frac{1}{\phi_{\tau}} \sum_g^G \sigma_{r_x,i\tau g} \phi_{\tau g}$

RMG Benchmark

  • The RMG was benchmarked against a Serpent run where both transport and essential physics methods were fed the OECD's LWR Burnup Credit reactor [2].
  • The benchmark examines the first year of the burn with an initial $^{235}$U enrichment of 4.5 [w/o].

RMG Benchmark

Table 11: Nearest Neighbors over Burn

days$p^\star=1$2345678910111213141516$\cdots$
01661771881119191221331020$\cdots$
40.61661771881991111221331020$\cdots$
81.11771661881992121111331020$\cdots$
1227171886161992010212313111$\cdots$
1628187171996162010313212414$\cdots$
2031881991771020616313414122$\cdots$
2431881991771020616313414122$\cdots$
2841991881020717414616313515$\cdots$
3241991020818717414515313616$\cdots$
3651020199188717515414313616$\cdots$

RMG Benchmark

Figure 20: Multiplication Factor Benchmark

RMG Benchmark

(a) BOL(b) EOL

Figure 21: Neutron Flux Spectrum Benchmarks

RMG Benchmark

(a) $^{235}$U(b) $^{246}$Cm(c) $^{137}$Cs

Figure 22: Selected Mass Fractions

RMG Benchmark

(a) $^{235}$U(b) $^{239}$Pu(c) $^{99}$Tc

Figure 23: Selected One-Group Cross-Sections

Multicomponent Enrichment Model

Here are the matched abundance ratio cascade equations:

$\frac{x_j^P}{x_j^F}\cdot\frac{P}{F} - \frac{\beta_j^{N_T+1} - 1}{\beta_j^{N_T+1} - \beta_j^{-N_P}} = 0$
$\left(\frac{x_j^F}{\frac{T}{F}} \cdot \frac{1 - \beta_i^{-N_P}}{\beta_j^{N_T+1} - \beta_j^{-N_P}} \right) - ( x_j^T\cdot\sum_i^{I} x_i^T ) = 0$
$\min\left[\frac{L}{F}\right]\to \frac{d}{dM^*} \frac{L}{F} = 0$

What if we solved these equations symbolically, rather than numerically..

Multicomponent Enrichment

...and use code generation to produce a blazing solver.

Multicomponent Enrichment

We get really big speedups (20-300x)!

(a) per iteration(b) over minimzation

Figure 23: Execution Timing Comparison

See http://pynesim.org/usersguide/enrichment.html for more details.

Fuel Cycle Sensitivity Study

Sensitivty Study Overview

  • Here, the system-wide impact of physical parameter perturbations is quantified.
  • This is done by looking for linear, 1D sensitivity coefficients for each parameter to the repository capacity response, measured in [PWh].
  • Denote these sensitivity coefficients as $S_x$ for a small change in an input parameter $x$.
  • Sensitivity coefficients are presented for 10% changes in $x$ from assumed base-case values.

Sensitivity Methodology Algorithm

1. For each input parameter $x$, change its base-case value $x_0$ by $\pm10\%$.

$x \to (1 \pm 0.1)x_0$

2. If $x$ is a separation efficiency (SE), instead perturb by $\pm$ "one nine" worth of separations.

$x = 0.999 \to x \in \{0.99, 0.9999\}$

3. Meanwhile, maintain all other parameters ($y, z, \ldots$) at their base-case values.

$y, z, \ldots \to y_0, z_0, \ldots$

4. Calculate the new repository capacity response, $R$, by the perturbed cycle and record the sensitivity as the percent change from the base-case response, $R_0$.

$S_{\pm x} = \left(\frac{R}{R_0} - 1\right) \times 100$

Fuel Cycle

The essential physics models from before were used to study a Fast Reactor, Light Water Reactor symbiotic scenario analogous to Scheme 3a of a 2006 OECD report [3].

Fuel Cycle

With perturbable components, over 30 independent physical parameters may be adjusted in the fuel cycle. Sensitivity results are computed from equilibrium values derived from a full treatment of the preceding, transient cycles.

Base Case Parameter Definition

Input Parameter $x$ValueUnits
LWR Burnup50.0MWd/kgIHM
LWR Fuel to Moderator Ratio0.301
LWR UF Storage Time60years
SE of U from LWR UF0.999
SE of NP from LWR UF0.999
SE of PU from LWR UF0.999
SE of AM from LWR UF0.999
SE of CM from LWR UF0.999
SE of CS from LWR UF0.999
SE of SR from LWR UF0.999
FR Burnup140.0MWd/kgIHM
FR TRU Conversion Ratio0.50
Max Fraction of Lanthanide in FR Fuel0.0005Atoms/TRU Atom
FR UF Storage Time3years
Storage Before Disposal50years

Base Case Parameter Definition

Input Parameter $x$ValueUnits
SE of U from FR UF0.999
SE of NP from FR UF0.999
SE of PU from FR UF0.999
SE of AM from FR UF0.999
SE of CM from FR UF0.999
SE of CS from FR UF0.999
SE of SR from FR UF0.999
Density of Host Rock2580kg/m$^{3}$
Specific Heat of Host Rock840J/kg-K
Thermal Conductivity of Host Rock1.626W/m-K
Heat Loss Factor During Ventilation0.7
Drift diameter5.5m
Ventilation System On Time50years
Ambient Environment Temperature20C
Distance Between Drifts81m

Fuel Cycle Base Case Benchmark

Scheme 3aNEA [3]Model$^{1}$% DiffModel$^{2}$% Diff
Electricity Share: LWR0.6320.619459-2.02440.634907+0.4579
Electricity Share: FR0.3680.380541+3.29550.365093-0.7962
FR UF: U0.6980.713806+2.21430.715224+2.4082
FR UF: NP0.00650.00661961+1.80700.00685174+5.1335
FR UF: PU0.2660.248059-7.23270.248319-7.1204
FR UF: AM0.020.022679611.81520.0217317+7.9687
FR UF: CM0.00980.00883517-10.9200.00787319-24.4730
HLW: U0.0133240.0132681-0.42130.0134448+0.8984
HLW: NP2.26542E-052.4079E-05+5.91732.41083E-05+6.0316
HLW: PU0.0007047970.000658893-6.96680.000632361-11.4548
HLW: AM5.03426E-055.63068E-0510.59235.12344E-05+1.7405
HLW: CM2.18151E-051.99031E-05-9.60681.67094E-05-30.5563
HLW: FP0.9858760.985973+0.00980.985831-0.0046
  1. Model with initial LWR $^{235}$U enrichment of 4.2 w/o.
  2. Model with LWR discharge burnup of 50 MWd/kg.

Fuel Cycle Dynamics (Quasi-Statics?)

Figure 7: Input Streams to FR [kg/kgIHM]

Fuel Cycle Dynamics

Figure 8: FR Output [kg/kgd] & TRU CR

Fuel Cycle Dynamics

Figure 9: $^{237}$Np, $^{239}$Pu Input and Output [kg/kgIHM]

Fuel Cycle Dynamics

Figure 10: $^{241}$Am, $^{244}$Cm Input and Output [kg/kgIHM]

Sensitivity Results

Table 3: Top Parameters from Screening Study

Parameter $x$Base Case $x_0$$S_{-x}$$S_{+x}$
FR_SE_PU0.999-83.09+99.30
LWR_SE_PU0.999-60.90+18.28
FR_SE_AM0.999-59.53+17.55
FR_SE_CM0.999-11.15+1.28
FR_TRU_CR0.500+5.22-6.02
Rock_Density2580 [kg/m$^{3}$]-5.05+5.59
Rock_Specific_Heat840 [J/kg/K]-4.90+5.59
LWR_BUd50 [MWd/kg]+5.46-3.48
FR_BUd140 [MWd/kg]-5.44+5.32

Sensitivity Results

Table 4: Bottom Parameters from Screening Study

Parameter $x$Base Case $x_0$$S_{-x}$$S_{+x}$
LWR_SE_U0.999-0.01+0.00
FR_SE_NP0.999-0.10+0.01
FR_SE_U0.999-0.22+0.01
FR_LAN_FF_Cap0.0005 [w/o]+0.03-0.03
LWR_SE_NP0.999-0.26+0.03
LWR_SNF_Storage_Time6 [y]+0.06-0.14
LWR_SE_CM0.999-0.82+0.08
Vent_System_On_Time50 [y]-0.11+1.13
Drift_Diameter5.5 [m]+0.26+0.58

Counterintuitive Results

  • For drift diameter, both $S_{-x}$ and $S_{+x}$ are positive!
    • This arises from the assumption that the drift diameter very small as compared to the spacing between drifts.
  • Since the base-case sits near a local drift diameter minimum, linear sensitivities will never capture the full effect on the response.
  • Additionally egregious is that, with increasing LWR burnup, the repository capacity goes down!
    • This due to the overall TRU production per unit energy produced from the LWRs declining.

LWR Burnup Sensitivity

Information Theoretic Analysis

Entropy Based Methodology

  • Here again, the system-wide impact of physical parameter perturbations is quantified.
  • This is done by performing Contingency Table analysis for each parameter to the repository capacity response.
  • Denote fuel cycle responses as $R$ [GWh] for $x$ & $y$ input parameters.
  • Goal: Determine strength of association between inputs and pairs of inputs to the response in a manner that is independent of:
    • the functional form of the response to the input $R(x)$,
    • and any base-case set of initial inputs.

Entropy Based Methodology

Input Parameter Definition

Input Parameter $x$MinMaxUnitsSample Type
LWR Burnup Level30.080.0MWd/kgIHMlinear
LWR Fuel to Moderator Ratio0.280.36linear
LWR SNF Storage Time330yearslinear
SE of U from LWR UF0.990.9999nines
SE of NP from LWR UF0.90.9999nines
SE of PU from LWR UF0.90.9999nines
SE of AM from LWR UF0.90.9999nines
SE of CM from LWR UF0.90.9999nines
SE of CS from LWR UF0.90.9999nines
SE of SR from LWR UF0.90.9999nines
FR Burnup Level100.0200.0MWd/kgIHMlinear
FR TRU Conversion Ratio0.250.95linear
Max Fraction of Lanthanide in FR Fuel0.00010.005Atoms/TRU Atomlinear
FR UF Storage Time330yearslinear
Storage Before Disposal1300yearslog

Input Parameter Definition

Input Parameter $x$MinMaxUnitsSample Type
SE of U from FR UF0.990.9999nines
SE of NP from FR UF0.90.9999nines
SE of PU from FR UF0.90.9999nines
SE of AM from FR UF0.90.9999nines
SE of CM from FR UF0.90.9999nines
SE of CS from FR UF0.90.9999nines
SE of SR from FR UF0.90.9999nines
Density of Host Rock23172869kg/m$^{3}$linear
Specific Heat of Host Rock5901270J/kg-Klinear
Thermal Conductivity of Host Rock1.92043.2856W/m-Klinear
Heat Loss Factor During Ventilation0.8060.914linear
Drift diameter4.56.5mlinear
Ventilation System On Time10300yearslog
Ambient Environment Temperature12.8232.82Clinear
Distance Between Drifts56106mlinear

Entropy Based Methodology

  • Before attempting to comprehend the subtleties of a 30+ dimensional surface, prudence demands a check to see if all 30 variables are really needed...
  • (Hint: probably not!)
  • Thus a quantitative ranking of the associations of each input parameter to the response is desired.
  • (Note that these associations do not imply a linear response!)
  • These rankings are obtained by borrowing an analysis tool that is often used in Biology: Contingency Tables.

Contingency Tables

The $2\times 2$ table is most common:

BlondeBrunetteTotals
Female181735
Male111425
Totals293160

But doesn't this approach ignore the underlying biology?

Yes!

Contingency Tables

Contingency Tables

Contingency Tables

Contingency Tables

Contingency Tables

Contingency Tables

Fuel Cycle Contingency Tables

For example, let's construct a 2D contingency table that measures the response from a sample input: fast reactor fuel plutonium separation efficiency, FR_SE_PU.

$0.9<$SE$<0.99$$0.99<$SE$<0.999$$0.999<$SE$<0.9999$
$10^4 <$ Capacity $< 10^5$$739$$43$$27$$809$
$10^5 <$ Capacity $< 10^6$$31510$$21611$$19469$$72590$
$10^6 <$ Capacity $< 10^7$$2648$$13095$$14430$$30173$
$10^7 <$ Capacity $< 10^8$$0$$213$$1053$$1266$
$34897$$34962$$34979$$104838$

Contingency Table Statistics

There are several metrics that have been developed to measure associations with contingency tables.

The primary measure here is the Uncertainty Coefficient $U(x|R)$.

To calculate $U(x|R)$ the following are needed:

  • The Entropy $H(x)$, which is a measure of how evenly the data is spread out in the $x$ parameter.
  • The Mutual Information $I(R,x)$ which states the shared value of the $x$ together with $R$.

Contingency Table Statistics

Denote the entries in a contingency table as $N_{ij}$ for the $i$-th response bin and the $j$-th input bin.

The probability of landing in a given bin is therefore,

$p_{ij} = \frac{N_{ij}}{N}$

Thus the entropy is calculated from,

$H(x) = - \sum_j^J p_{\cdot j} \ln(p_{\cdot j})$

The mutual information is found via,

$I(R,x) = - \sum_{i,j}^{I,J} p_{ij} \ln\left(\frac{p_{ij}}{p_{i\cdot}\cdot p_{\cdot j}}\right)$

Visual Relationship

Figure 11: Entropy Relations [4]

Contingency Table Statistics

The uncertainty coefficient is then calculated from,

$U(x|R) = \frac{I(R,x)}{H(x)}$

This metric has the following useful properties:

  1. Defined on the range $[0, 1]$.
  2. $U(x|R) = 0$ implies that $I(R,x) = 0$, which indicates that the parameter is unassociated with the response.
  3. $U(x|R) = 1$ requires that $I(R,x) = H(x)$. This implies that the system response $R(x)$ is solely determined by $x$.

Input Parameters Ranked by $U(x|R)$

Rank$x$$U(x|R)$
1FR_SE_PU0.07667
2HLW_Storage_Time0.05264
3FR_SE_AM0.04148
4Heat_Loss_Factor0.01548
5LWR_SE_PU0.01343
6FR_TRU_CR0.008895
7LWR_SE_AM0.007304
8FR_BUd0.003773
9LWR_UF_Storage_Time0.003551
10Rock_Specific_Heat0.003085
11FR_SE_CM0.002522
12Rock_Thermal_Conductivity0.001954
13Ambient_Temp0.001353
14LWR_BUd0.001053
15LWR_SE_CS0.001033

Input Parameters Ranked by $U(x|R)$

Rank$x$$U(x|R)$
16LWR_SE_SR0.001024
17FR_UF_Storage_Time0.0005559
18Drift_Space0.0004421
19LWR_SE_U0.0002899
20Rock_Density0.0002718
21FR_SE_CS0.0001331
22FR_SE_SR0.0001073
23Vent_System_On_Time0.0001016
24Drift_Diameter9.648E-05
25LWR_SE_NP9.575E-05
26LWR_SE_CM7.969E-05
27LWR_Fuel2Mod7.833E-05
28FR_LAN_FF_Cap6.559E-05
29FR_SE_NP6.212E-05
30FR_SE_U6.207E-05

$U(x|R)$ Rankings vs $U(R|x)$ Rankings

$U(x|R)$ is related to $U(R|x)$ via the following expression:

$U(x|R) = \frac{H(R)}{H(x)}U(R|x)$

This is effectively an entropic re-expression of Bayes' Theorem:

$p(x|R) = \frac{p(x)}{p(R)}p(R|x)$

Due to the stochastic driver, all $H(x)$ must be the same - to within error - for all $x$. Also, since there is only one repsonse, $H(R)$ is the same the for all $x$.

Hence, $U(x|R) = k U(R|x)$.

Scaling the uncertainty coefficients by a constant does not affect the relative rank ordering.

Entropy Based Covariance

  • Now that the association of $x$ to $R$ may be determined, finding the sensitivity of $R(x)$ to other inputs.
  • Call $y$ an input parameter distinct from $x$. Measuring the simultaneous effect of two inputs to a response may be performed using 3D tables!
  • Note that 'slices' of these 3D tables may be seen as 2D tables bins along a specific cut of data.
  • For 30 inputs, there are $_{30}C_2 = 435$ contingency tables to construct.
  • Different metrics will be needed to account for the added dimensionality.

FR Pu SE and HLW Storage Time

Slice for 1 < HLW_Storage_Time < 6.694

$0.9<$SE$<0.99$$0.99<$SE$<0.999$$0.999<$SE$< 0.9999$
$10^4 <$ Capacity $< 10^5$$419$$23$$14$$456$
$10^5 <$ Capacity $< 10^6$$11145$$9970$$9553$$30668$
$10^6 <$ Capacity $< 10^7$$110$$1556$$2085$$3751$
$10^7 <$ Capacity $< 10^8$$0$$0$$0$$0$
$11674$$11549$$11652$$34875$

Slice for 6.694 < HLW_Storage_Time < 44.814

$0.9<$SE$<0.99$$0.99<$SE$<0.999$$0.999<$SE$< 0.9999$
$10^4 <$ Capacity $< 10^5$$273$$19$$11$$303$
$10^5 <$ Capacity $< 10^6$$10859$$7527$$6373$$24759$
$10^6 <$ Capacity $< 10^7$$484$$4157$$5175$$9816$
$10^7 <$ Capacity $< 10^8$$0$$2$$24$$26$
$11616$$11705$$11583$$34904$

FR Pu SE and HLW Storage Time

Slice for 44.814 < HLW_Storage_Time < 300

$0.9<$SE$<0.99$$0.99<$SE$<0.999$$0.999<$SE$< 0.9999$
$10^4 <$ Capacity $< 10^5$$47$$1$$2$$50$
$10^5 <$ Capacity $< 10^6$$9506$$4114$$3543$$17163$
$10^6 <$ Capacity $< 10^7$$2054$$7382$$7170$$16606$
$10^7 <$ Capacity $< 10^8$$0$$211$$1029$$1240$
$11607$$11708$$11744$$35059$

Aggregation over HLW_Storage_Time

$0.9<$SE$<0.99$$0.99<$SE$<0.999$$0.999<$SE$< 0.9999$
$10^4 <$ Capacity $< 10^5$$739$$43$$27$$809$
$10^5 <$ Capacity $< 10^6$$31510$$21611$$19469$$72590$
$10^6 <$ Capacity $< 10^7$$2648$$13095$$14430$$30173$
$10^7 <$ Capacity $< 10^8$$0$$213$$1053$$1266$
$34897$$34962$$34979$$104838$

Covariance Statistics

  • While $U(x,y|R)$ is roughly equivalent to the concept of mutual sensitivity, 3D tables should allow for an information theoretic parallel to covariance.
  • In general for contingency tables covariance is not well defined, however here there is the special stochastic situation where $x$ and $y$ are known a priori to be independent of $R$.
  • Therefore, the following new coefficient of variation metric is proposed.
  • $c_v(x|y|R)$ is the 'sensitivity of sensitivity', or how much an input affects the response of another input.

Covariance Statistics

First, define:

$U(x|R)|y = \left\{ \left.U(x|R)\right|_{l_0}^{l_1} , \left.U(x|R)\right|_{l_1}^{l_2} , \ldots \left.U(x|R)\right|_{l_{C-1}}^{l_C} \right\}$

Then the coefficient of variation is,

$c_v(U(x|R)|y) = \frac{\sigma(U(x|R)|y)}{\mu(U(x|R)|y)}$

Finally, the symmetric $c_v$ is set as:

$c_v(x|y|R) = \frac{1}{2} \left(c_v(U(x|R)|y) + c_v(U(y|R)|x)\right)$

Covariance Statistics

$c_v(x|y|R)$ has the following properties:

  • Defined on the range $[0, 1]$ since $0 \le U \le 1$.
  • $c_v = 0$ implies that $\sigma = 0$, which indicates that $U(x|R)$ shows no association with $y$. Thus there are no covariant effects observed.
  • $c_v = 1$ indicates that $\sigma = \mu$. This connotes that the value of $x$ solely governs the response $R$ from $y$.

Top Covariant Parameter Pair Rankings

Rank$x$$y$$c_v(x|y|R)$
1FR_SE_AMHLW_Storage_Time0.02318
2FR_SE_AMFR_SE_PU0.02316
3HLW_Storage_TimeLWR_UF_Storage_Time0.01557
4FR_SE_PUHLW_Storage_Time0.01556
5HLW_Storage_TimeLWR_SE_PU0.01155
6FR_SE_AMFR_TRU_CR0.01061
7FR_SE_PULWR_UF_Storage_Time0.008192
8FR_SE_PULWR_SE_PU0.008
9FR_BUdFR_SE_AM0.007716
10Ambient_TempFR_SE_PU0.007489
11Ambient_TempHLW_Storage_Time0.007367
12HLW_Storage_TimeLWR_SE_AM0.007248
13Ambient_TempFR_SE_CS0.007235
14Ambient_TempLWR_Fuel2Mod0.007066
15Ambient_TempHeat_Loss_Factor0.007044

Covariant Parameter Pair Example

Figure 11: Total & Top Contributors to Decay Heat [Watts/kg] of HLW for (FR_SE_AM, FR_SE_PU).

Ensuring Good Statistics

How do we know that we have picked a sufficient number of bins?

Different partitioning schemes may significantly alter our measures. Consider a a chamber the moment after a piston has been suddenly removed, as in a Carnot heat engine:

Figure 12: Abstract Piston with Representative Partitions: Blue Line for 2 Bins, Red Line for 5 Bins.

Ensuring Good Statistics

We thus want to pick the minumum number of bins after which our metrics (and rankings) have stablized.

Figure 13: Dynamic Effects of Binning Structure Example: $U(x,y|R)$ for the parameter pair (FR_SE_PU, LWR_SE_NP) to repository capacity.

Ensuring Good Statistics

But how many runs are needed to achieve this binning?

For the covariance measures (which determines the statistics), note that:

$I(x,y) \to 0$ as $N \to \infty$.

Set $t$ as a fractional statistical tolerance, ie $t=0.01$ is 1% error.

Then the mutual information goes as $I = t^d$ for $d$ indpendent inputs. Here, $d = 2$ for $x$ & $y$.

For $J$ number $x$ bins and $K$ number of $y$ bins, the total number of runs required is estimated by:

$N \approx \frac{1}{t^2}JK$

Obligatory Open Source Rant (Do it! Do it now!)

- PyNE, Cyclus, OpenMC, Fudge, Software Carpentry, etc.

- June 24th - 29th, Austin, TX! http://conference.scipy.org/scipy2013/

- http://numfocus.org/

Questions

References

  1. Anthony M. Scopatz and Erich a. Schneider. “A new method for rapid computation of transient fuel cycle material balances”. In: Nuclear Engineering and Design 239.10 (Oct. 2009), pp. 2169–2184. ISSN: 00295493. DOI: 10.1016/j.nucengdes.2009.02.022. URL: http://linkinghub.elsevier.com/retrieve/pii/S0029549309000922
  2. M Takano. Burnup Credit Criticality Benchmark - Result of Phase 1A - NEA/NSC/DOC(93)22. Tech. rep. January. Nuclear Energy Agency, 1994, pp. 1–147
  3. OECD. Advanced Nuclear Fuel Cycles and Radioactive Waste Management - NEA-5990. Tech. rep. Nuclear Energy Agency, 2006, pp. 1–248
  4. http://en.wikipedia.org/wiki/File:Entropy-mutual-information-relative-entropy-relation-diagram.svg
  5. Nuclear Energy Data Files. Los Alamos National Laboratory - Theoretical Nuclear Physics. 2009. URL: http://t2.lanl.gov/data/data.html
  6. Jaakko Lepp. PSG2 / Serpent âA ̧ a Continuous-energy Monte Carlo Reactor Physics Burnup Calculation Code. Tech. rep. 2011, p. 157
  7. A G Croff. ORIGEN 2.2 - CCC-371. Tech. rep. Oak Ridge National Laboratory, 2002, pp. 1–222