Use of probability tables for propagating uncertainties in neutronics

Annals of Nuclear Energy 102 (2017) 376–385

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Use of probability tables for propagating uncertainties in neutronics M. Coste-Delclaux ⇑, C.M. Diop, S. Lahaye Den-Service d’étude des réacteurs et de mathématiques appliquées (SERMA), CEA, Université Paris-Saclay, F-91191 Gif-sur-Yvette, France

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 5 August 2016 Received in revised form 23 November 2016 Accepted 26 November 2016 Available online 13 January 2017 Keywords: Random variable Probability tables Uncertainties Nuclear data Monte Carlo simulation

Probability tables are a generic tool that allows representing any random variable whose probability density function is known. In the field of nuclear reactor physics, this tool is currently used to represent the variation of cross-sections versus energy (neutron transport codes TRIPOLI4Ò, MCNP, APOLLO2, APOLLO3Ò, ECCO/ERANOS. . .). In the present article we show how we can propagate uncertainties, thanks to a probability table representation, through two simple physical problems: an eigenvalue problem (neutron multiplication factor) and a depletion problem. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction A random variable n, defined on a domain ðDÞ, whose probability density function (pdf) pðnÞ is known, can be represented by a moment-based probability table of order Nfðpi ; ni Þ; i ¼ 1; Ng where the probabilities pi and the steps (or bands) ni are calculated in order to preserve 2N moments:

Z Mn ¼

ðDÞ

nn pðnÞdn ¼

N X

pi nni

n ¼ L; . . . ; 2N þ L  1

where L is the order (positive or negative) of the first moment to be preserved. To ensure probability normalization and mean value preservation, we need to exactly perform zero-order and oneorder moments. That leads to 2  2N 6 L 6 0. In the neutronics field, when dealing with nuclear reaction cross-sections, the most accurate probability tables are obtained by preserving negative and positive moments (the value L ¼ 1  N is the one recommended by P. Ribon (Ribon and Maillard, 1986), positive moments giving a good representation of cross-section peaks and negative moments a good representation of cross-section depressions. But, in the present work, in order to calculate the moments of any random variable, we chose to preserve only positive moments (L ¼ 0Þ. As a matter of fact, negative moments of a Gaussian law do not exist (Piegorsch and Casella, 1982).

E-mail address: [email protected] (M. Coste-Delclaux). http://dx.doi.org/10.1016/j.anucene.2016.11.044 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.

 either by a Monte Carlo technique using the exact probability density function,  or by a multiband technique using the bands and the probabilities of the probability table.

ð1Þ

i¼1

⇑ Corresponding author.

The physical interpretation of such a probability table is that it provides a set of values ni ‘‘attained” by the random variable and the associated probabilities. Consequently, any random variable can be sampled:

Physical data, such as nuclear reaction cross-sections or radioactive decay constants are given in nuclear data bases by their mean values and their standard deviations. We generally make the assumption that their associated probability density functions follow a Gaussian law. When such ‘‘uncertain data” are the input of a physical equation E (slowing-down equation or depletion equation for instance), they induce a distribution on the output of equation E. Classical perturbation methods give the mean value and the standard deviation of the output of E but not the probability distribution of E itself. In this paper, we show how such a distribution can be carried out using probability tables. Section 2 describes the moment-based probability table calculation and their mathematics properties. Sections 3 and 4 deal with uncertainty propagation in the context of neutronics field, respectively for an eigenvalue and a depletion problem. This work is a continuation of a previous work done by the authors (Diop et al., 2012). It extends the concept of using probability tables for propagating uncertainties established for a fixed particle source problem to more general problems: eigenvalue and depletion problems.

377

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385 k X bj M kj ¼ ak

2. Moment-based probability table

06k6N1

ð9Þ

j¼0

2.1. Calculation

System (9) shows that: To compute the moment-based probability tables, we use a Padé approximant technique (Ribon and Maillard, 1986). It can be summarized in the following algorithm. Let us consider the function:

pðnÞ dn 1  zn

IðzÞ ¼ ðDÞ

in the neighborhood of z ¼ 0. We can give two expressions of IðzÞ. On one hand, using relations (1):

Z

pðnÞ dn ¼ 1  zn

IðzÞ ¼ ðDÞ

2N1 X

M n zn þ Oðz2N Þ

2N1 N XX

pi ðzni Þn þ Oðz2N Þ

N 2N1 X X ¼ pi ðzni Þn þ Oðz2N Þ

ð2Þ

N X pi i¼1

1 þ Oðz2N Þ 1  zni

and on another hand, using a ½N  1=N Padé approximant:

Z

pðnÞ dn ¼ 1  zn

IðzÞ ¼ ðDÞ

2N1 X n¼0

ð3Þ

½N1=N

where: N 1 X PN1 ðzÞ ¼ aj z j

ð4Þ

j¼0

and:

Q N ðzÞ ¼

N X bj z j

ð5Þ

j¼0

Eq. (3) shows that the k-order monomials for 0 6 k 6 2N  1 must be equal in the expressions:

B¼

2N1 X

!

Mn z

n¼0

n

N X bj z j

Consequently, if zero is a k -order root of Q N ðzÞ, zero is also a k -order root of PN1 ðzÞ. By supposing that the Padé approximant is reduced (no zero roots), we can always choose:

b0 ¼ 1

N X bj M kj ¼ M k

!

ð6Þ

j¼0

ð7Þ

j¼0

N1 X PN1 ðzÞ al ¼ aN þ Q N ðzÞ 1  zbl l¼1

ð8Þ

j¼0

Once the coefficients bj are known, the coefficients aj of the numerator of the Padé approximant ½N  1=N are obtained by identifying the k-order monomials for 0 6 k 6 N  1. That leads to the system of N equations, N unknowns:

ð11Þ

ð12Þ

In formulas (11) and (12), bl are the inverses of Q N ðzÞ roots. These roots are simple because they are the roots of orthogonal polynomials as shown in the continuation of this paper. By considering the equality obtained by comparing formulas (2) and (3): N X pi i¼1

1 PN1 ðzÞ þ Oðz2N Þ ¼ þ Oðz2N Þ 1  zni Q N ðzÞ

ð13Þ

and by substituting the Padé approximant by its decomposition given either by Eq. (11) or by Eq. (12), we can take:  when LQ ¼ N

ni ¼ bi ; pi ¼ ai

;

i ¼ 1; . . . ; N

ð14Þ

 when LQ ¼ N  1 and LP ¼ N  1,



The coefficients bj of the denominator of the Padé approximant ½N  1=N are obtained by identifying the k-order monomials for N 6 k 6 2N  1. All these monomials must be equal to zero. That leads to the system of N equations, ðN þ 1Þ unknowns: N X bj M kj ¼ 0 N 6 k 6 2N  1

ð10Þ

When LQ ¼ N  1 and LP ¼ N  1, which happens for symmetric distributions in a domain ðDÞ centered in zero when N takes odd values (cf. Appendix A), the Padé approximant can be decomposed into:



and N1 X A¼ aj z j

N 6 k 6 2N  1

j¼1

N PN1 ðzÞ X al ¼ Q N ðzÞ 1  zbl l¼1

M n zn þ Oðz2N Þ

PN1 ðzÞ ¼ þOðz2N Þ Q ðzÞ |fflfflfflfflN{zfflfflfflffl}

0 6 l 6 kg

Once system (10) is solved, system (9) can be solved and the exact orders of PN1 ðzÞ and Q N ðzÞ are found. Let us denote LP and LQ the respective orders of PN1 ðzÞ and Q N ðzÞ. When LQ ¼ N, which is the general case, the Padé approximant can be decomposed into:

n¼0

i¼1

0 6 l 6 kg ) fal ¼ 0;

System (8) can be rewritten as a system of N equations, N unknowns:

n¼0

n¼0 i¼1

¼

and more generally:

fbl ¼ 0;

Z

¼

b0 ¼ 0 ) a0 ¼ 0

ni ¼ bi ; pi ¼ ai

;

i ¼ 1; . . . ; N  1

ð15Þ

and



nN ¼ 0; pN ¼ aN

ð16Þ

Some very efficient algorithms exist for calculating the roots of a polynomial, and then, the bands of the probability tables are calculated by Eqs. (14) or (15) and (16). But for the probability calculation, we prefer solving the linear system, once the bands are known:

378

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

N X pi nni ¼ M n

n ¼ 0; . . . ; N  1

ð17Þ

i¼1

This linear system can be solved without pivoting by using the known analytical inverse of a Vandermonde matrix. We obtain: N X pj ¼ ljk M k1

Z n

l

! N1 N1 X X jþ1 pðnÞdn ¼ bN1j n bN1j M lþjþ1

D

j¼0

ð26Þ

j¼0

ð18Þ

k¼1

where the factor ljk is the coefficient of the k-order monomial of the Lj ðnÞ Lagrange polynomial given by:

Lj ðnÞ ¼

The scalar product of RN with a monomial of degree l such that 0 6 l 6 N  1, is equal to:

N N Y n  nn X ¼ ljk nk1 n  nn k¼1 n¼1 j

ð19Þ

According to Eq. (8), as bN is here equal to zero, this scalar product is equal to zero for N 6 N þ l 6 2N  1, that is to say 0 6 l 6 N  1. In that case also, the RN polynomials constitute an orthogonal polynomial suite for the scalar product defined by (22). 2.3. Histogram of n distribution

n–j

This analytical approach is the one used in the APOLLO2 neutron transport code (Hébert and Coste, 2002) for calculating cross-section probability tables. 2.2. Gaussian quadrature rule According to Eq. (1), a moment-based probability table can also be seen as a N-point Gaussian quadrature rule since it leads to a quadrature rule that is exact for polynomials of degree less or equal to ð2N  1Þ. That confers to them very interesting mathematical properties. The main ones are:  The quadrature of any polynomial of order strictly less than 2N is exact.  The ni and pi values are always real  If pðnÞ remains positive on ðDÞ:  The ni values are inside ðDÞ.  More precisely, we can find a partition of ðDÞ composed of intervals Di ¼ ½ei1 ; ei  that satisfy the following properties:

ni 2 Di Z pðnÞdn ¼ pi

ð20Þ ð21Þ

Di

All these properties come from the fact that the ni values are the roots of orthogonal polynomials for the scalar product defined by:

Z

hf ; gi ¼

f ðnÞgðnÞpðnÞdn

ð22Þ

D

As a matter of fact, when LQ ¼ N, the ni values are the inverses of the roots of Q N ðzÞ, and therefore are the roots of the polynomial of degree N, RN ðzÞ defined by:

RN ðzÞ ¼

N N1 X X bNj z j ¼ bNj z j þ zN j¼0

ð23Þ

j¼0

If we calculate the scalar product of RN with a monomial of degree l such that 0 6 l 6 N  1, we obtain:

Z nl D

! N N X X bNj n j pðnÞdn ¼ bNj M lþj j¼0

ð24Þ

j¼0

And according to Eq. (8), this scalar product is equal to zero for N 6 N þ l 6 2N  1, that is to say 0 6 l 6 N  1. The RN polynomials constitute an orthogonal polynomial suite for the scalar product defined by (22). In the same way, when LQ ¼ N  1 and LP ¼ N  1, the ni values are the roots of the polynomial of degree N, RN ðzÞ defined by: N 1 N1 X X RN ðzÞ ¼ z bN1j z j ¼ bN1j zjþ1 j¼0

j¼0

ð25Þ

Properties (20) and (21) allow us to interpret a moment-based probability table as a histogram of n distribution whose bounds are ei values and whose heights are pi values. In the continuation of this paper, this specific histogram will be called ‘‘intrinsic histogram” of the probability table. The difficulty consists in performing the most precise values of the bounds ei . In that purpose, the probability density function is approximated by an analytical form allowing us to solve Eq. (21). 2.3.1. Polynomial representation A simple analytical form is a polynomial one. We approach pðnÞ by:

pðnÞ ¼

K X c k nk

ð27Þ

k¼0

The coefficients ck are calculated in order to preserve ðK þ 1Þ moments of n distribution. For a moment-based probability table at order N, ðK þ 1Þ will not exceed 2N. The coefficients are obtained by solving the linear system:

Z nn pðnÞdn  D

Z K X ck nnþk dn ¼ M n k¼0

06n6K

ð28Þ

D

Once the ck are known, we solve step by step the equation:

Z

ei

K X ck nk dn ¼ pi

ð29Þ

ei1 k¼0

That leads to find the roots of a ðK þ 1Þ-order polynomial and to choose, if they exist, the one(s) in the convenient interval. The existence and the uniqueness of the solution are not insured. 2.3.2. Piecewise constant representation Another simple analytical form of pðnÞ is a piecewise constant representation (Ribon et al., 2008):

pðnÞ ¼ qk

fk1 6 n 6 fk

16k6K

ð30Þ

The ð2K þ 1Þ unknowns fqk ; 1 6 k 6 Kg and ffk ; 0 6 k 6 Kg are calculated in order to preserve ð2K þ 1Þ moments of n distribution. For a moment-based probability table at order N, ð2K þ 1Þ will not exceed 2N. The equations to solve are: For 0 6 n 6 2K

Z nn pðnÞdn  D

Z K X qk k¼1

fk

fk1

nn dn ¼

K X fnþ1  fnþ1 k1 qk k ¼ Mn nþ1 k¼1

ð31Þ

This system reads: For 0 6 n 6 2K

q1 fnþ1 þ 0

K1 X ðqk  qkþ1 Þfnþ1 þ qK fnþ1 ¼ ðn þ 1ÞMn k K k¼1

ð32Þ

379

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

We define ð2K þ 2Þ pseudo-moments by:



It is that latter method that has been used in the following applications for representing the probability density functions.

T0 ¼ 0 T n ¼ nM n1

ð33Þ

1 6 n 6 2K þ 1

3. Eigenvalue problem

and we set:

8 > < > :

r 0 ¼ q1 r k ¼ qk  qkþ1

16k6K 1

ð34Þ

r K ¼ qK We can note that:

K X

rk ¼ 0

ð35Þ

k¼0

0 6 n 6 2K þ 1

AðnðaÞ Þ depending on K independent random variables nðaÞ 1 6 a 6 K; which probability density functions pðaÞ are known and we look for the distribution of the greatest eigenvalue. The problem to be solved is then:

AðnðaÞ ÞXðnðaÞ Þ ¼ kðnðaÞ ÞXðnðaÞ Þ where kðn Þ is the greatest eigenvalue.

ð36Þ

3.1. Multiband equation

k¼0

ðaÞ

System (36) is formally similar to system (1) and is solved in the same way. As a matter of fact, we consider the generating function:

Z

z

JðzÞ ¼ ðDÞ

ð1  znÞ2

ð1 þ xÞ ¼ 1 þ

ðaÞ

ated with the random variables nðaÞ and Dia the interval associated ðaÞ nia ,

n!

n¼0

ð1Þ

ðKÞ

Di ...Di 1

x

n

ð37Þ

Z

¼

for ða ¼ 2Þ we can write:

1 ð1  znÞ

2

¼1þ

we can write the following so-called multiband equation:

Z

1 X aða  1Þ . . . ða  n þ 1Þ

1 X ð2Þð3Þ . . . ðn þ 1Þ

n!

n¼0

ð1Þ ðKÞ Di ...Di K 1

ð38Þ

¼

X i1 ...iK ¼

ð1Þ ðKÞ Di ...Di K 1

K Y

pðaÞ ðnðaÞ ÞdnðaÞ

ð44Þ

a¼1 ðKÞ

AðnðaÞ Þ  Aðni1 ; . . . ; niK Þ ð1Þ

ðKÞ

ð45Þ

Eq. (43) becomes: ð1Þ

ð39Þ

z ð1  znÞ2

pðnÞdn ¼

2Kþ1 X

n¼0

T n zn þ Oðz2Kþ2 Þ

R

k¼0

1 þ Oðz2Kþ2 Þ 1  zfk

kðnðaÞ ÞXðnðaÞ Þ

qK ¼ r K 16k6K 1

YK kðnðaÞ Þ a¼1 pðaÞ ðnðaÞ ÞdnðaÞ K 1 YK R pðaÞ ðnðaÞ ÞdnðaÞ ð1Þ ðKÞ D ...D a¼1

ð1Þ

ð40Þ

The great difference is that the fk values are no longer the roots of orthogonal polynomials because the T n quantities do not represent moments. Consequently, there is no mathematical argument insuring that the fk values are inside ðDÞ. We deduce the qk from the rk values by:

qk ¼ r k þ qkþ1

ð1Þ ðKÞ Di ...Di K 1

K Y

pðaÞ ðnðaÞ ÞdnðaÞ

a¼1

Under the assumption that, in Eq. (46), the eigenvalue kðnðaÞ Þ can be replaced by its mean value:

n¼0

K 2Kþ1 K X X X n rk zn fk þ Oðz2Kþ2 Þ ¼ rk ¼



Z

ðKÞ

Aðni1 ; . . . ; niK ÞX i1 ...iK ¼

ð46Þ

Moreover, using Eq. (36):

k¼0

ð43Þ

(20)):

n¼0

PK ðzÞ ¼ þ Oðz2Kþ2 Þ Q Kþ1 ðzÞ

ðDÞ

pðaÞ ðnðaÞ ÞdnðaÞ

a¼1

n¼0

n¼0

JðzÞ ¼

XðnðaÞ Þ

ð1Þ

2K 2Kþ1 X X T nþ1 znþ1 þ Oðz2Kþ2 Þ ¼ T n zn þ Oðz2Kþ2 Þ

Z

K Y

and for nð1Þ 2 Di1 ; . . . ; nðKÞ 2 DiK we make the approximation (cf. Eq.

z

pðnÞdn 2 ðDÞ ð1  znÞ # Z "X 2K ¼ ðn þ 1Þznþ1 nn þ Oðz2Kþ2 Þ pðnÞdn ðDÞ

kðnðaÞ ÞXðnðaÞ Þ

Z

Consequently:

Z

K Y pðaÞ ðnðaÞ ÞdnðaÞ

a¼1

K

n¼0

JðzÞ ¼

AðnðaÞ ÞXðnðaÞ Þ

We define:

ðznÞn

1 X ¼1þ ðn þ 1ÞðznÞn

ðaÞ

By denoting fðpia ; nia Þ; ia ¼ 1; N ðaÞ g the probability table associ-

to

pðnÞdn

in the neighborhood of z ¼ 0. Using the expansion: a

ð42Þ

ðaÞ

With such definitions system (32) and Eq. (35) become: K X r k fnk ¼ T n

The first application of propagating uncertainties thanks to probability tables is an eigenvalue problem. We consider a matrix

ð41Þ

Once the piecewise constant approximation of the probability density function is known, we solve step by step the equation giving the ei values. The existence of the solution is not insured but, if it exists, it is unique.

ki1 ...iK ¼

ðKÞ

Di ...Di

i1

ð47Þ

iK

we obtain: ð1Þ

ðKÞ

Aðni1 ; . . . ; niK ÞX i1 ...iK ¼ ki1 ...iK X i1 ...iK

ð48Þ

Q One has to solve a¼1;K N ðaÞ eigenvalue problems in order to obtain the distribution of the greatest eigenvalues fki1 ...iK g. The fki1 ...iK ; 1 6 i1 6 N ð1Þ ; . . . ; 1 6 iK 6 N ðKÞ g with their associated Q ðaÞ probabilities Ka¼1 pia give a probabilistic distribution of the ran-

dom variable kðnðaÞ Þ. The mean value and the variance of the random variable kðnðaÞ Þ can be approximated by Eqs. (49) and (50).

380

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

X Z

E½k ¼

ð1Þ ðKÞ Di ...Di K 1

i1 ¼1;N ð1Þ ... iK ¼1;N ðKÞ

X

¼

... iK ¼1;N ðKÞ

X

ðKÞ

1

K

ki1 ...iK

i1 ¼1;N ð1Þ ... iK ¼1;N ðKÞ

X

V½k ¼

ð1Þ

Di ...Di K Y

a¼1

3.2. Results

a¼1

Z

i1 ¼1;N ð1Þ

¼

K Y kðnðaÞ Þ pðaÞ ðnðaÞ ÞdnðaÞ

ki1 ...iK

K Y

pðaÞ ðnðaÞ ÞdnðaÞ

a¼1

ð49Þ

ðaÞ

pia

For testing the method, we consider a stationary one-point problem for neutron transport in an infinite homogeneous medium with two energy groups, epithermal (group 1) and thermal (group 2) (Duderstadt and Hamilton, 1976). The bound between the two energy groups is taken so that neutrons cannot be up-scattered from group 2 to group 1. The neutron flux U is given by the equation:

QU ¼

ðki1 ...iK  E½kÞ2

i1 ¼1;N ð1Þ ... iK ¼1;NðKÞ

M n ½k ¼

i1 ¼1;N

ð1Þ

... iK ¼1;NðKÞ

kni1 ...iK

K Y

a¼1

ð52Þ

with:

K Y ðaÞ pia

ð50Þ

a¼1

Q¼ 

The multiband calculation allows us to estimate the moments of the random variable k at any order n by the formula:

X

1 PU k

ðaÞ

pia

ð51Þ

and to build a probability table that gives a very precise and concise representation of k distribution.

P¼

Ra1 þ R1!2 s

0

R1!2 s

Ra2

!

v1 m1 Rf 1 v1 m2 Rf 2 v2 m1 Rf 1 v2 m2 Rf 2

ð53Þ  ð54Þ

Ra1 and Ra2 are the macroscopic absorption cross-sections in group 1 and group 2, R1!2 is the macroscopic transfer crosss section from group 1 to group 2, m1 Rf 1 and m2 Rf 2 are the production cross-sections in group 1 and in group 2, v1 and v2 are the fission spectra in group 1 and group 2.

Table 1 Gaussian laws describing macroscopic the cross-sections.

Mean value li Standard deviation/mean value rli

nð1Þ ¼ m1 Rf 1

nð2Þ ¼ m2 Rf 2

nð3Þ ¼ Ra1

nð4Þ ¼ Ra2

nð5Þ ¼ R1!2 s

5.9752103 4.5%

1.0694101 3%

8.8383103 4.5%

6.8707102 6%

1.8330102 9%

i

Fig. 1. Comparison between the 10-order multiband histogram and the ‘‘intrinsic histogram”.

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

381

Due to the choice of the bound between the two energy groups, we can assume that all the fission neutrons appear in epithermal group and consequently take v2 ¼ 0 and v1 ¼ 1. Eq. (52) is equivalent to the eigenvalue problem:

These values come from a homogeneous condensation of a 900 MW PWR assembly.

Q 1 PU ¼ kU

1. We carry out the reference calculation that is a Monte Carlo calculation consisting in 103 batchs of 106 elementary eigenvalue calculations in which the coefficients of the random matrix

ð55Þ

with:

0

m1 Rf 1

B Ra1 þ R1!2 s B Q 1 P ¼ B @ m1 Rf 1 R1!2 s

1

m2 Rf 2

R

Ra2 ðRa1 þ R1!2 Þ R s

1!2 C a1 þ s C C 1!2 A 2 f2 s 1!2 Þ a2 ð a1 þ s

R

mR R R

ð56Þ

R

1

The matrix Q P is singular. One eigenvalue is then equal to zero and the second one gives the infinite multiplication factor k. The flux U is the eigenvector associated with the eigenvalue:

k¼

m 1 Rf 1 Ra1 þ R

1!2 s

þ

m2 Rf 2 R1!2 s

ð57Þ

Ra2 ðRa1 þ R1!2 Þ s

We consider the five independent random variables:

8 ð1Þ > > >n > > ð2Þ > > > > > nð4Þ > > > : ð5Þ n

¼ m1 Rf 1 ¼ m2 Rf 2 ¼ Ra1

ð58Þ

¼ Ra2 ¼ R1!2 s

and we assume that they respectively follow the normal distributions Gðli ; ri Þ with the mean values li and the standard deviations ri given in Table 1.

The stages of the test are as follows.

are sampled using for the random variables nðaÞ the Gaussian laws given in Table 1. It results, from this calculation, a reference distribution for the infinite multiplication factor, allowing us to calculate its three first centered moments and their associated statistical uncertainties. 2. We carry out the multiband calculation by describing the Gaussian laws on nðaÞ with probability tables at order 10. That leads to a multiband distribution, called 10-order multiband distribution, for the infinite multiplication factor described by 105 values and their associated probabilities. 3. The 10-order multiband distribution is a discrete distribution characterized by its moments. It can be represented itself by a moment-based probability table. We perform then the various moments of the 10-order multiband distribution and we look for the most precise moment-based probability table. The highest order that can be obtained is a 6-order. 4. We look for the ‘‘intrinsic histogram” associated to the 6-order probability table representing the 10-order multiband distribution. The grid of this histogram has been calculated by taking a piecewise constant approximation with 5 constant values for the probability density function. 5. We calculate the histogram of the 10-order multiband distribution (105 values with their associated probabilities) in the grid

Fig. 2. Comparison between the reference histogram and the ‘‘intrinsic histogram”.

382

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

defined by the ‘‘intrinsic histogram” and we compare this histogram, called 10-order multiband histogram, with the ‘‘intrinsic histogram” (Fig. 1). 6. We calculate the histogram of the reference distribution (109 values) in the grid defined by the ‘‘intrinsic histogram” and we compare this histogram, called reference histogram, with the ‘‘intrinsic histogram” (Fig. 2). 7. We compare the 10-order multiband histogram and the reference histogram (Fig. 3). 3.2.1. Results on moments We first study globally the results by comparing the three first centered moments of the reference and of the 10-order multiband (MB) distributions. Table 2 gives in the first column, the three first centered moments of the reference distribution, in the second column their associated one standard deviation and in the third column, the absolute discrepancy between the moments calculated on the 10-order multiband distribution and those calculated on the reference distribution. We see that these moments are very precisely reproduced by the multiband calculation, the absolute discrepancies versus the reference calculation being always less or equal to one standard deviation.

For eigenvalues between 1.092 and 1.434, the relative discrepancies do not exceed 6.1%. In Fig. 3, by comparing directly the histograms of the 10-order multiband and of the reference distributions within the grid of the ‘‘intrinsic histogram”, we check the agreement between the two distributions. It is excellent since, for 99.96% of the distribution, the relative discrepancies do not exceed 1.9%. One has to keep in mind that for the reference calculation, 109 eigenvalue problems are solved while for the multiband calculation only 105 eigenvalue problems are taken into account. 4. Depletion problem The second application of propagating uncertainties thanks to probability tables is a depletion problem. Let us solve the stochastic differential equation:

dNðt; kÞ ¼ kNðt; kÞ dt

ð59Þ

where k is a random variable whose probability density function pðkÞ is known. We suppose that:

Nð0; kÞ ¼ N0

ð60Þ

The theoretical solution of this equation is: 3.2.2. Results on distributions Fig. 1 gives the comparison between the ‘‘intrinsic histogram” and the 10-order multiband histogram, absolute discrepancy in the high part and relative discrepancy in the low part. For eigenvalues between 1.092 and 1.434 (96% of the 10-order multiband distribution), the relative discrepancies do not exceed 5.4%, which proves the efficiency of the probability table representation and the correctness of the probability density function approximation. The same kind of comparison was also performed by using the reference distribution (109 values) instead of the 10-order multiband distribution. The results, given in Fig. 2, are very similar.

Nðt; kÞ ¼ N 0 expðktÞ

ð61Þ

Consequently, if k follows a Gaussian law Gðm; rÞ, Nðt; kÞ follows a Log-Normal law whose mean value and variance are respectively given by:

E½Nðt; kÞ ¼ N0 expðtm þ

t2 r2 Þ 2

ð62Þ

and:

V½Nðt; kÞ ¼ N20 ðexpð2tm þ t 2 r2 ÞÞðexpðt 2 r2 Þ  1Þ

Fig. 3. Comparison between the reference and the 10-order multiband histograms. Table 2 Three first centered moments of the eigenvalue distributions for the reference and the multiband calculations.

1-order moment (Mean value) 2-order moment 3-order moment

Reference distribution moments

One standard deviation on the reference distribution moments

Absolute discrepancy between the reference and the 10-order multiband distribution moments

1.27307642 5.845358103 1.405776104

2.37106 2.80107 4.45108

3.04107 ± 2.37106 2.28107 ± 2.80107 4.46108 ± 4.45108

ð63Þ

383

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

4.1. Multiband equation

V½Nðt; kÞ ¼

i¼1

By denoting fðpi ; ki Þ; i ¼ 1; Kg the probability table representing the random variables k and Di the interval associated to ki , we can write the multiband equation derived from Eq. (59):

Z Di

dNðt; kÞ pðkÞdk ¼  dt

Z kNðt; kÞpðkÞdk

ð64Þ

Di

We define:

Nðt; kÞpðkÞdk

ð65Þ

Di

and for k 2 Di , we make the approximation (cf. Eq. (20)):

k  ki

dNi ðtÞ ¼ ki Ni ðtÞ dt

Z ð69Þ

Di

The quantities

Ri ðtÞ ¼ N0 expðki tÞ

ð70Þ

associated to the probabilities pi give the probabilistic distribution of the random variable Nðt; kÞ. The mean value and the variance of Nðt; kÞ can be approximated by Eqs. (71) and (72). K Z X i¼1

¼ N0

The moments of the random variable Nðt; kÞ at any order n are given by the formula:

Nðt; kÞpðkÞdk ¼

Di

K X Ni ðtÞ i¼1

K X pi expðki tÞ ¼ pi Ri ðtÞ

K X i¼1

ð71Þ

i¼1

K X pi ðRi ðtÞÞn

ð73Þ

i¼1

allowing us to build a probability table corresponding to the Nðt; kÞ distribution.

ð67Þ

ð68Þ

Nð0; kÞpðkÞdk ¼ N0 pi

ð72Þ

i¼1

We treat here the decay rate of a radionuclide whose disintegration constant is denoted k. We suppose that k follows the Gaussian law G(0.00385, 0.0004) and that its initial concentration is equal to unit. We study Nðt; kÞ distribution for various times t:

with:

E½Nðt; kÞ ¼

K X pi ðRi ðtÞ  E½Nðt; kÞÞ2

4.2. Results

Differential Eq. (67) admits the solution:

Ni ðtÞ ¼ Ni ð0Þexpðki tÞ

ðNðt; kÞ  E½Nðt; kÞÞ2 pðkÞdk Di

ð66Þ

Eq. (64) becomes:

Ni ð0Þ ¼

¼

Mn ½Nðt; kÞ ¼

Z Ni ðtÞ ¼

K Z X

T T 2T ð30sÞ; ð60sÞ; ð120sÞ; Tð180sÞ; 2Tð360sÞ; 5Tð900sÞ 6 3 3 where T is the half-life corresponding to the mean value of k. For each time t, we carry out a reference calculation by a Monte Carlo method consisting in 103 batchs of 106 elementary calculations (Eq. (61)) in which the disintegration constant is sampled following its Gaussian law. We also perform a multiband calculation in which the Gaussian law followed by k is described by a 24-order probability table (24 elementary calculations (Eq. (68))). At the various times t, Table 3 gives in the first column, the two first centered moments of the Log-Normal law described by Eqs. (62) and (63), in the second column, the two first centered moments of the reference distribution, in the third column the associated one standard deviation and in the fourth column, the absolute discrepancy between the moments calculated on the 24-order multiband distribution and those calculated on the reference distribution. The results are all in excellent agreement. The absolute discrepancies between the 24-order multiband distribution and the reference distribution are always less or equal to two standard deviations as shown in Table 3.

Table 3 Two first centered moments of the various distributions at various times. One standard deviation on the reference distribution moments

Absolute discrepancy between the reference and the 24-order multiband distribution moments

t ¼ T6 0.890985152 1.1432184 104

3.35 107 4.86 109

4.31 107 ± 3.35 107 1.35 109 ± 4.86 109

0.793968096 3.632065 104

t ¼ T3 0.793968863 3.632022 104

5.98 107 1.55 108

-7.68 107 ± 5.98 107 4.40 109 ± 1.55 108

1-Order moment (Mean value) 2-Order moment

0.630748544 9.176887 104

t ¼ 2T 3 0.630749759 4 9.176767 10

9.50 107 3.95 108

-1.21 106 ± 9.50 107 1.20 108 ± 3.95 108

1-Order moment (Mean value) 2-Order moment

0.50137147 1.30650297 103

t¼T 0.50137291 1.30648437 103

1.13 106 5.71 108

-1.44 106 ± 1.13 106 1.86 108 ± 5.71 108

1-Order moment (Mean value) 2-Order moment

0.25267985 1.3377555 103

t ¼ 2T 0.25268128 1.3377301 103

1.15 106 6.28 108

-1.43 106 ± 1.15 106 2.54 108 ± 6.28 108

1-Order moment (Mean value) 2-Order moment

0.033366596 1.540547 104

t ¼ 5T 0.033367025 1.540473 104

3.90 107 1.04 108

-4.29 107 ± 3.90 107 7.39 109 ± 1.04 108

Log-normal law moments

Reference distribution moments

1-Order moment (Mean value) 2-Order moment

0.890984721 1.1432317 104

1-Order moment (Mean value) 2-Order moment

384

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

Distribuon of the concentraons at various mes Ref t=30 MB t=30 Ref t=60 MB t=60 Ref t=120 MB t=120 Ref t=180 MB t=180 Ref t=360 MB t=360 Ref t=900 MB t=900

600000 500000 400000 300000 200000 100000 0 0

0.2

0.4

0.6

0.8

1 N/N0

Fig. 4. Comparison between the reference histogram and the ‘‘intrinsic histogram” at various times (s).

At each time t, the distribution obtained by the 24-order multiband calculation is represented by a probability table that gives an ‘‘intrinsic histogram”. The histogram of the reference distribution is then calculated on the grid associated to this ‘‘intrinsic histogram”. It is called reference histogram. The 24-order multiband distribution (normalized to the reference distribution) and the reference distributions are in a good agreement for all times t as shown in Fig. 4. The relative discrepancy between the central parts of the two distributions does not exceed 10% for each studied time. Once again, this example shows the power of the probability table representation, 24 ‘‘chosen” calculations giving a very good approximation of one billion calculations.

System (10) reads: N X bj M kj ¼ M k

2N0 þ 1 6 k 6 4N0 þ 1

ð75Þ

j¼1 0

For odd values of k ¼ 2k þ 1, only bj terms with odd j have a non-zero coefficient: N0 X 0 b2j0 þ1 M2k0 2j0 ¼ 0 N0 6 k 6 2N0

ð76Þ

j0 ¼0 0

while for even values of k ¼ 2k , only bj terms with even j have a non-zero coefficient: N0 X b2j0 M 2k0 2j0 ¼ M 2k0

0

N0 þ 1 6 k 6 2N0

ð77Þ

0

j ¼0

5. Conclusion Probability tables appear as a generic and powerful tool for representing and propagating uncertainties, giving the complete probability distribution of the quantity of interest, unlike the classical perturbation methods that only give the mean value and the variance. It has been successfully applied to a fixed source problem (Diop et al., 2012), an eigenvalue problem and a depletion problem. The strength of this tool is that it can be used to describe, in a very accurate and concise way, either the input random data, whatever the probability law is, or the output data, which allows to chain calculations. The next step will be to extend this work to correlated random variables. One method would be to use correlated moment-based probability tables calculated by using multidimensional Pade approximants but related theoretical problems have to be solved prior.

System (10) splits into two systems, one for bj terms with odd j and one for bj terms with even j. System (76) is a homogeneous system with non-zero determinant. Consequently, all bj terms with odd j are equal to zero. That shows that, in such a case, the order of the polynomial Q N ðzÞ is ðN  1Þ and Q N ðzÞ only contains even monomials. Consequently, if zi is a root of Q N ðzÞ, (- zi Þ is also a root of Q N ðzÞ. The steps of the N -order probability table are zero and 2N 0 symmetric values. In the same way, system (9) reads: k X bj M kj ¼ ak

0 6 k 6 2N0

ð78Þ

j¼0 0

For odd values of k ¼ 2k þ 1, only bj terms with odd j have a non-zero coefficient: N0 X b2j0 þ1 M2k0 2j0 ¼ a2k0 þ1

0

0 6 k 6 N0  1

ð79Þ

0

j ¼0

Acknowledgment

0

The authors are grateful to Dr. J.M. Martinez and Dr. F. Bachoc for their helpful advice. They also thank Dr. C. Jouanne and Dr. S. Mengelle for providing the numerical data.

Let us consider a N-order probability table, N being odd, representing a symmetric distribution in a domain ðDÞ centered in zero. All odd moments of this distribution are equal to zero. We set:

N ¼ 2N þ 1

N0 X b2j0 M 2k0 2j0 ¼ a2k0

0

0 6 k 6 N0

ð80Þ

j0 ¼0

Appendix A.

0

while for even values of k ¼ 2k , only bj terms with even j have a non-zero coefficient:

ð74Þ

Eq. (79) shows that all ak terms with odd k are equal to zero. The order of polynomial P N1 ðzÞ is ðN  1Þ and PN1 ðzÞ only contains even monomials. Eq. (12) can be written: 0

0

N N X X a0l PN1 ðzÞ al þ ¼ aN þ Q N ðzÞ 1  zbl l¼1 1 þ zbl l¼1

ð81Þ

M. Coste-Delclaux et al. / Annals of Nuclear Energy 102 (2017) 376–385

As:

PN1 ðzÞ PN1 ðzÞ ¼ ¼ aN þ Q N ðzÞ Q N ðzÞ

References N0 X l¼1

al

1 þ zbl

þ

N0 X l¼1

a0l

1  zbl

ð82Þ

we can deduce:

al ¼ a0l 1 6 l 6 N0

385

ð83Þ

In the N-order probability table, the probability of symmetric steps are the same. In the same way it is very easy to show that a N-order probability table, N being even, representing a symmetric distribution in a domain ðDÞ centered in zero has N=2 symmetric steps with equal probabilities.

Ribon, P., Maillard, J.M., 1986. Les tables de probabilité: Application au traitement des sections efficaces pour la neutronique. Note CEA-N-2485. Piegorsch, W.W., Casella, G., 1982. The Existence of Negative Moments of Continuous Distributions. Biometrics Unit, Cornell University, Ithaca, NY 14853. Diop, C.M., Coste-Delclaux, M., Lahaye, S., 2012. Cross-section uncertainty propagation in the multigroup slowing down equations using probability table formalism. Nucl. Sci. Eng. 170, 87–97. Hébert, A., Coste, M., 2002. Computing moment-based probability tables for selfshielding calculations in lattice codes. Nucl. Sci. Eng. 142, 245–257. Ribon, P., Coste-Delclaux, M., Jouanne, C., Diop, C.M., 2008. Angular anisotropy representation by probability tables. PHYSOR 2008, Interlaken, Switzerland. September 14–19. Duderstadt, J.J., Hamilton, L.J., 1976. Nuclear Reactor Analysis. John Wiley & Sons.