- Email: [email protected]

0306-4549/89$3.00+0.00 Copyright (~ 1989PergamonPressplc

OPTIMIZATION OF BIASING PARAMETERS IN PARTICLE TRANSPORT SIMULATIONS R. INDIRA'~ Institut ffir Reaktorentwicklung der Kernforschungslange Jfilich, Postfach 1913, D-5170 Jfilich, F.R.G. (Received 14 June 1989)

Abstraet--A technique for optimization of biasing parameters in Monte-Carlo simulation of particle transport problems is described. This technique is a variant of Spanier's multi-stage procedure. It is based on stochastic evaluation of the derivatives of the second moment using small samples. Application of the technique to tracklength biasing schemes employed in deep penetration problems is discussed. Various leakage estimators are considered to demonstrate the applicability of the technique. Optimum biasing parameters are obtained for various values of scattering probability and are compared with those obtained through Spanier's technique. The advantages of our technique over Spanier's technique are discussed.

l. INTRODUCTION

function of medium properties (Levitt, 1968 ; Ponti, 1971; Dubi and Dudziak, 1979; Indira and John, Variance reduction techniques are imperative in the 1986 ; Indira, 1988b,c). These studies have shown that Monte-Carlo simulation of deep penetration probthe variance is a sensitive function of the biasing parlems (Spanier and Gelbard, 1969; Murthy, 1982). ameter and have indicated the need to optimize the These techniques usually bias one or more of the biasing parameter for practical problems. random variables in the simulation so that the variThe first attempt in this direction was made by ance of the parameter of interest is reduced. The biasing Spanier (1971). His procedure is based on estimating of the random variable is usually controlled by an the parametric dependence of variance from small arbitrary parameter in the importance function, statistical samples of the actual problem. Starting with known as the biasing parameter. The efficiency of a guess value of the parameter, simulation is carried biasing depends on the proper choice of the biasing out and variance estimated as a function of biasing parameter. Various biasing schemes employed in deep parameter. Based on this, the biasing parameter is penetration problems and the difficulties encountered modified and the procedure is repeated until the biasin making correct choice of the biasing schemes are ing parameter converges. Normally, three to four iterwell discussed in detail, see for example, Murthy ations were found to be adequate for convergence. (1980) and Indira (1988a). Spanier's procedure is very useful for practical In this paper, we shall restrict ourselves to the probproblems. However, we found that when applied to lem of optimization of biasing parameters in importance biasing schemes for particle transport problems. problems such as photon streaming in ducts, the biasUsually, a good guess of the biasing parameter is ing parameter did not converge (Murthy, 1980 ; Indira made based on physical insight or by optimization and Murthy, 1981; Murthy and Indira, 1983). Our using simplified models of the problem on hand (Mac- studies showed that this is due to the absence of a Millan, 1966; Spanier, 1970; Murthy, 1982, 1983; well-defined minimum in variance as a function of Murthy and Indira, 1986). The studies using sim- biasing parameter, when the scattering probability of plified models have been found useful and adequate the medium is low. We also observed that when there for the purpose of obtaining the general behaviour of is a flat minimum, the biasing parameter oscillates. Of variance as a function of biasing parameter and for course, this p e r se does not pose a serious problem, characterizing the efficiency of biasing schemes as a since it implies that in this particular case, variance does not depend sensitively on the biasing parameter near its optimum value. This would, however, be very t On leave from Radiation Shielding and Statistical Physics Section, Reactor Physics Division, Indira Gandhi Centre disturbing when the technique is applied in a practical for Atomic Research, Kalpakkam 603102, Tamil Nadu, problem. India. Another disadvantage of Spanier's technique is that 571

572

R. |NDIRA

the variance is estimated at discrete values of the biasing parameter. When the variance depends sensitively on the biasing parameter (i.e. when the minimum is quite steep), there is a possibility of obtaining a nonoptimal parameter, thus resulting in a less efficient simulation. The number of values at which variance is estimated cannot be increased much, as it would increase the computational time. The important point is that the shape of the parametric dependence of variance is not fully utilized in Spanier's method. Subsequently, MacMillan (1972) suggested a modification which incorporates the functional dependence, by using the first and second derivatives. These derivatives are estimated numerically using finite difference approximation. This was also found to have difficulties owing to the associated numerical errors (Indira, 1988d). In this paper, we propose a new variant to Spanier's technique. We obtain formal expressions for the first and second derivatives of the second moment of transmission with respect to the biasing parameter. These are then stochastically estimated using small samples. Two to three iterations are found adequate for convergence of biasing parameter to its optimal value. Our studies of typical transmission problems with media of varying scattering properties showed that there is a uniform convergence of the biasing parameter to its optimal value, when our technique is employed. Comparison with actual Monte-Carlo simulation with a large number of histories has confirmed the correctness of our procedure. The details of our technique are given in Section 2. In Section 3, we consider application to transport problems, where exponential biasing is employed in conjunction with various leakage estimators. In Section 4, a model deep penetration problem is considered and the results are compared with those of Spanier's technique. The principal conclusions of our studies are brought out in Section 5.

1" h2(x)f2(x)

M2(~,l~. . . .

)

= jg2(x,a,

fl,::.)g(x,a,

fl .... )dx.

(1)

Let us now consider the dependence of M2 on a single parameter, say, a. As in Spanier's technique, we start with a string of guess values of biasing parameters, a0, fl0. . . . . and we have :

M2( , 0 .... ) = I

3g(x,~o,~o .... )

h(x)f(x) , x g(x,~,flo . . . . ) gtx, COo,flu. . . . ) dx.

(2)

Consider now the derivatives of the second moment with respect to a. Extending this to derivative with respect to all parameters is straightforward and hence we shall restrict ourselves to only one parameter ~, and drop the other parameters from the argument. Let M'2(a) and M'~(a) denote the first and second derivative of M2 with respect to c~. Then : Ih(x)f(x) h(x)f(x) g'(x, ~) M'2(a) = J g(x, ao) g(x,a) g(x,a) g(x, ao)dx,

(3)

Ih(x)f(x) h(x)f(x) M'~(a) = j g(x,c¢o) g(x,~)

×k

(g(x,~)J

g(x,a)jgt×,~0)ux. (4)

Here, g'(x, a) and g"(x, c¢) are the first and second derivatives of g(x, ~) with respect to a. In the above expressions for M~(:¢) and M'~(a), the integral over the probability space has to be carried out by actual simulation of the problem. Let us now discuss how equations (3) and (4), though formal, are used in obtaining the optimum biasing parameter. By Taylor expansion to second-order, we can write : M2(a) = M2(ao) + M'2(ao)aa+ ~M'~(Oto)(aa)2.

(5)

Dividing by 5e and rearranging the terms, we get : 2. D E R I V A T I V E S O F S E C O N D M O M E N T

Let f(x) denote formally the joint density function of all the random variables of the problem under simulation. (This is dependent on the material properties. But we shall not explicitly write them, as they do not affect our discussion.) Let g(x, c~,fl,...) denote the biased joint density, where ~, fl,.., is the string of biasing parameters. These determine the strength of biasing of each random variable. Let us denote the score function for the required parameter as h(x). Then the second moment of h(x) in the biased simulation is given by :

M'2(oO = M'2(ao) + ½M';(a0)6~.

(6)

Now, at the optimum value of the biasing parameter, M'z(~) should be zero. Hence : M~(ao)

5or = -- 2 M'~ (~0)'

(7)

c%is modified by 5a and the simulation is again carried out using the new value ofa. The procedure is repeated until convergence in c~is obtained. We note that equations (3) and (4) are only the formal expressions for M':(c0 and M';(~). Depending on the biasing schemes, the form of biasing function

Optimization of biasing parameters in particle transport simulations and the type of estimator used, one has to derive the relevant expressions for these, before they are employed in actual simulation. In the next section, we obtain the expressions for M~(~) and M'~(a) for the well-known and commonly employed exponential biasing scheme proposed by Kahn (1956). We also consider various types of estimators normally employed in transmission problems, see for example, Coveyou et al. (1967). 3. A P P L I C A T I O N TO EXPONENTIAL BIASING SCHEME

Let us consider a deep penetration problem of monoenergetic particle transport through a slab shield. For simplicity, we consider the medium to have unit total cross section and a scattering probability C with isotropic scattering. The source is incident in the/t = + 1 direction on the left face and we are interested in estimating the transmission through the right face. Let us consider a Monte-Carlo simulation with exponential biasing and obtain the optimum biasing parameter for this problem. Exponential biasing implies that we sample the intercollision distance x, from a biased density function g(x, ~t) of the form ~ exp ( - ~ x ) instead of from f ( x ) given by exp ( - x ) . ~ is of the form (1 - b ~ ) and our interest is to optimize b. We shall consider three types of estimators namely : (i) binomial estimator or analogue estimator, which scores the actual particle weight transmitted ; (ii) last event estimator, which scores on the occurrence of an absorption event. The score is given by the weight of the particle at the time of absorption multiplied by the probability of the particle to be transmitted uncollided through the right face ; and (iii) next event estimator, which scores at every collision event. The score from a collision is given by the probability of the particle to be transmitted uncollided through the right face multiplied by the particle weight at the collision instant. The score from a particle history is given by the sum over the collision scores. We shall now obtain the expressions for the second moment of transmission and its derivatives, M2(b), M'2(b) and M~(b), for the above estimators.

573

until the time of escape through the right face. Hence M2(b) is given by : Mz(b) =

( {. g(x,,

CN i~=1--

0,)j

C

n

11

(8) In the above and in what follows, the angular brackets denote an average over an ensemble of all possible particle histories. The first and second derivatives of M2(b) are given by : Mi(b)=

\

t~

1 1 - -

L,= l 0(x,. ~0,)J g(xi,~i)j ti~_l g(xi,~i ) db j]// ,

M~(b) =

(9)

/(x,) l

l l ~ x ~ \ c L,=,g(x,. o,)] .ffm lrf "

).

+,~=,: g(x,.~,) ~

(i0)

Substituting for f(xi), g(xi, ~0i) and g(xi, ui) and noting that ~i = 1 - b#i, we get :

m2(b)=

×L

[A <7 ']>

CN

exp

I~iXi

,

(11)

_.F f ~ exp(-bu,xl)} [[

Lli=,

--

•

~i

.~

3.1. Analogue estimator We shall consider survival biasing to be employed. This means that the particle is never absorbed. The particle weight is multiplied by the scattering probability C at every collision. As we sample the intercollision distance from g(xi, ~i), the factor f(x)/g(x, ~t) is given by the product overf(x~)/g(xi, ~g), where xi is the ith intercollision distance sampled, h(x) is given by C N, where N is the number of collisions suffered

ANE16:1~-D

c-[{,__fi

+ (,_Ll m IX i

From equations (11-13), we see that we need carry only one product term and two summation terms

574

R. INDIRA

in the course of a history, since the derivatives are estimated at b = b0 only. At the end of the history the two square bracketed terms are computed and the second moment and its derivatives are scored. After a particular number of histories are generated the modified value of b is computed employing equation (7). The procedure is repeated to get the converged value of b.

3.2. Last event estimator In this case, we do not employ survival biasing. The score h(x) is given by C e x p - ( T - x r ) , where x~. is the site of absorption event. Proceeding as in the case of binomial estimator, we find that the expressions for Mz(b), M'2(b) and M~(b) are same as equations (11 13) with C N replaced by Cexp--(T--xL) and N replaced by NL, where NL is the number of collisions suffered by the particle before it is absorbed.

3.3. Next event estimator Normally, survival biasing is employed when the next event estimator is used. Hence, we have a score of C"exp-(T-X,) at the nth collision event, where X, is the site of the nth collision event. Obviously, 2(, = E"_ ~pixy. Let N be the total number of collisions suffered by the particle at the time of escape from the medium through the right face. Then :

M=(b) =

([ ]~l C"exp- (T-X,) I]" exp (-b0#,x,) l n

i= I

O~Oi

x [ L C"exp-(T-X,,)~_.=

X

f i exp

(14)

(--c~o/~x~)][~

i=l

n

{C"exp--(T-X,) I

" exp ~o~ ][.~=,{C"exp--(T-X.) x~_~ N

x fiexp(--bl~'x')~{X#-2)(,

~-

/1i

-J-i

]~i

,

4, M O D E L D E E P P E N E T R A T I O N

PROBLEM--RESULTS

AND DISCUSSIONS

We consider particle transport through a 30-m.f.p. thick shield, of scattering probability ranging from 0.1 to 0.9. The scattering is considered to be isotropic. A unit source is incident on the left face. Transmission through the right face is the parameter of interest. Analogue, last event and next event estimators are considered. Exponential biasing is employed and the optimum biasing parameter is obtained using the procedure described in Section 3. For comparison, we obtain the optimum biasing parameter and corresponding variance in transmission using Spanier's technique also. In Spanier's procedure, the second moment is estimated at b = 0.05, 0.1, 0.15,... ,0.95. 20,000 histories were generated in each iteration of the simulations using our technique and Spanier's technique. In Table 1, we present the optimum biasing parameter and corresponding variance in transmission for the case of exponential biasing with an analogue estimator. These have been obtained by the simulation of equations (11-13) for C= 0.1, 0.3, 0.5, 0.7 and 0.9. Also given are the corresponding values obtained using Spanier's technique. We see that the optimum biasing parameter obtained by our method is close to that obtained with Table 1. C o m p a r i s o n o f an o p t i m u m biasing parameter for an exponentially biased simulation with an analogue estimator for a 30-m.f.p. thick shield

M'~(b) = { [ ~, C" e x p - ( r - x,,) (--bo,uixi)

The above equations give the expressions for the derivatives of the second moment of transmission, and these are used along with equation (7) to obtain the optimum biasing parameter. We also note that the extra terms computed for the purpose of obtaining the derivatives are linear terms. Hence our method requires much less computer time compared to Spanier's scheme, wherein the second moment has to be computed for various values of b, each computation involving an exponential term. In the next section, we obtain the optimum biasing parameter and corresponding variance in transmission for a model problem.

L ~! (16)

C

b

M2(b)

bs

M2(b~)

0. I 0.3 0.5 0.7 0.9

0.96 0.95 0.88 0.79 0.53

0.411E-24 0.104E-23 0.419E-22 0.112E-18 0.216E-I I

0.95 a 0.95 a 0.90 0.80 0.50 ~

0.414E-24 0.116E-23. 0.799E-22 0.212E-18 0.407E-11

" Indicates the nonconvergence o f biasing parameter. ,5 = o p t i m u m biasing parameter obtained with the new scheme. ,5, = o p t i m u m biasing parameter obtained with Spanier's scheme.

Optimization of biasing parameters in particle transport simulations Table 2. Comparison of an optimum biasing parameter for an exponentially biased simulation with a last event estimator for a 30-m.f.p. thick shield C

b

M2(b)

b,

M2(bO

0.1 0.3 0.5 0.7 0.9

0.98 0.95 0.94 0.78 0.67

0.380E-22 0.146E-21 0.126E-20 0.239E-15 0.177E-11

0.95~ 0.90 0.90 0.70 0.60~

0.339E-22 0.230E-21 0.243E-20 0.623E-15 0.248E-11

"Indicates the nonconvergence of biasing parameter. b = optimum biasing parameter obtained with the new scheme. bs = optimum biasing parameter obtained with Spanier's scheme.

Spanier's m e t h o d . D u e to the fine t u n i n g o f b in o u r procedure, there is also a n i m p r o v e m e n t in the variance by a factor o f 1.1-2. It is seen t h a t the biasing p a r a m e t e r does n o t converge within three to four iterations for C = 0. !, 0.3 a n d 0.9 w h e n Spanier's m e t h o d is used. In o u r m e t h o d , convergence p r o b l e m s do n o t arise. In Tables 2 a n d 3, we present the same results for a n exponentially biased simulation with last event a n d next event estimator, respectively. W e find t h a t the o p t i m u m biasing p a r a m e t e r s o b t a i n e d using the two m e t h o d s are quite close. Again, we see t h a t for low a n d very high scattering probabilities, the o p t i m u m biasing p a r a m e t e r keeps oscillating w h e n Spanier's m e t h o d is employed. V a r i a n c e reduction is i m p r o v e d by a factor o f 1.2-2 in simulations with the o p t i m u m biasing p a r a m e t e r o b t a i n e d using o u r technique. A n explicit M o n t e - C a r l o simulation was carried out for the case o f T = 30 m.f.p, a n d C = 0.5. E x p o n e n t i a l biasing was e m p l o y e d a l o n g with the three estimators described earlier. The simulations were carried o u t using the o p t i m u m biasing p a r a m e t e r s o b t a i n e d using o u r procedure. In each simulation 100,000 histories were generated. The second m o m e n t o f transmission a n d its first a n d second derivatives were c o m p u t e d . C o m p a r i n g the results o f o u r p r o c e d u r e with those o f explicit M o n t e - C a r l o simulation, it was seen t h a t the second m o m e n t was m a t c h i n g within + 5 % a n d the

Table 3. Comparison of an optimum biasing parameter for an exponentially biased simulation with a next event estimator for a 30-m.f.p. thick shield C

b

M2(b)

bs

M2(bs)

0.1 0.3 0.5 0.7 0.9

0.94 0.91 0.81 0.67 0.48

0.141E-22 0.110E-21 0.873E-20 0.448E-17 0.322E-10

0.95" 0.90 0.80 0.70 0.50a

0.199E-22 0.127E-21 0.160E-19 0.579E-17 0.620E-10

a Indicates the nonconvergence of biasing parameter. b = optimum biasing parameter obtained with the new scheme. b~ = optimum biasing parameter obtained with Spanier's scheme.

575

first a n d second derivatives were m a t c h i n g within -t-10%.

5. CONCLUSIONS A technique for o p t i m i z a t i o n o f biasing p a r a m e t e r s in i m p o r t a n c e biasing schemes has been developed. This technique, based o n stochastic evaluation o f the derivatives o f the second m o m e n t , has been f o u n d to be very effective in o p t i m i z a t i o n of biasing p a r a m e t e r s of a t r a c k l e n g t h biasing scheme in particle t r a n s p o r t problems. This technique has various a d v a n t a g e s over Spanier's t e c h n i q u e : (i) for all values o f the scattering probability, the biasing p a r a m e t e r converges uniformly to its o p t i m u m value ; a n d (ii) since the second m o m e n t a n d its derivatives are estimated at only one value o f biasing p a r a m e t e r , the c o m p u t e r time consumed is less. T h e o p t i m u m biasing p a r a m e t e r s o b t a i n e d using the new technique c o m p a r e well with those o b t a i n e d using Spanier's technique for the model problem. The results o f o u r technique have also been c o m p a r e d against a M o n t e - C a r l o simulation with a large n u m b e r o f histories for a typical case a n d have been f o u n d to agree well. Acknowledgements--The author is thankful to P. Cloth and D. Filges for their interest and to the Institut ffir Reaktorentwicklung, Kernforschungsanlage J/ilich for its hospitality. REFERENCES

Coveyou R. R., Cain V. R. and Yost K. J. (1967) Nucl. Sci. Engn 9 27, 219. Dubi A. and Dudziak D. J. (1979) Nucl. Sci. Engng 70, 1. Indira R. (1988a) A new weight moments formulation for the study of statistical errors in Monte-Carlo radiation transport. Ph.D. Thesis submitted to Madras University, India. Indira R. (1988b) Ann. nucl. Energy 15, 67. Indira R. (1988c) Ann. nucl. Energy 15, 261. Indira R. (1988d) Unpublished. Indira R. and John T. M. (1986) Nucl. Sci. Engn9 94, 323. Indira R. and Murthy K. P. N. (1981) Monte-Carlo simulation of benchmark experiment on neutron transport in thick sodium. Report RRC-48, Reactor Research Centre, India. Kahn H. (1956) Application of Monte-Carlo. Report AECU-3259, The Rand Corporation. Levitt L. B. (1968) Nucl. Sci. Engn9 31, 500. MacMillan D. B. (1966) Nucl. Sci. Engng 26, 366. MacMillan D. B. (1972) Nucl. Sci. Engn9 48, 219. Murthy K. P. N. (1980) Tracklength biasing in Monte-Carlo radiation transport. M.Sc. Thesis, Bombay University, India. Murthy K. P. N. (1982) Proc. Wkshp Monte-Carlo Method, January 1~20, BARC, Bombay, India. Murthy K. P. N. (1983) Ann. nucl. Energy 10, 375. Murthy K. P. N. and Indira R. (1983) Matscience Conf. on

576

R. INDIRA

Probability Theory and Monte-Carlo Methods, Ootacamund, India. Murthy K. P. N. and Indira R. (1986) Nucl. Sci. Enyn 9 92, 482. Ponti C. (1971 ) Angular and tracklength distribution biasing in Monte-Carlo deep penetration problems. Report ORNL-RSIC-29, ORNL.

Spanier J. (1970) S I A M J. Appl. Math. 18, 172. Spanier J. (1971) A new multistage procedure for systematic variance reduction in Monte-Carlo. Report ORNL-RSIC29, ORNL. Spanier J. and Gelbard E. M. (1969) Monte-Carlo Principles in Neutron Transport Problems. Addison-Wesley, Reading, Mass.