- Email: [email protected]
A new power mapping method based on ordinary kriging and determination of optimal detector location strategy
Annals of Nuclear Energy 68 (2014) 118–123
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A new power mapping method based on ordinary kriging and determination of optimal detector location strategy Xingjie Peng a,b,⇑, Kan Wang a, Qing Li b a b
Department of Engineering Physics, Tsinghua University, Beijing, China Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, China
a r t i c l e
i n f o
Article history: Received 16 September 2013 Received in revised form 27 December 2013 Accepted 3 January 2014 Available online 2 February 2014 Keywords: Power mapping Ordinary kriging Optimal detector location Simulated annealing
a b s t r a c t The Ordinary Kriging (OK) method is presented that is designed for a core power mapping calculation of pressurized water reactors (PWRs). Measurements from DayaBay Unit 1 PWR are used to verify the accuracy of the OK method. The root mean square (RMS) reconstruction errors are kept at less than 0.35%, and the maximum reconstruction relative errors (RE) are kept at less than 1.02% for the entire operating cycle. The reconstructed assembly power distribution results show that the OK method is fit for core power distribution monitoring. The quality of power distribution obtained by the OK method is partly determined by the neutron detector locations, and the OK method is also applied to solve the optimal neutron detector location problem. The spatially averaged ordinary kriging variance (AOKV) is minimized using simulated annealing, and then, the optimal in-core neutron detector locations are obtained. The result shows that the current neutron detector location of DayaBay Unit 1 reactor is near-optimal. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Core power distribution monitoring in operating power reactors is very important in core surveillance, the 3-D power distribution is one of the basic operation parameters which can determine many other important parameters such as power peaking factor, enthalpy rising factor and quadrant tilt ratio used to evaluate the operation condition of reactor and the safe margin. The economy of reactor could be optimized if the real time 3-D power distribution is well obtained and used for surveillance and regulation. Most commercial power reactors in operation are equipped with fixed or movable in-core neutron detectors to obtain power distribution information. Many kinds of on-line monitoring systems, such as BEACON (Boyd and Miller, 1996) and GNF-ARGOS (Tojo et al., 2008), have been developed to estimate in-core power distributions using fixed in-core detectors. The detector results at certain locations reflect the actual reactor flux or power can be applied to improve the results of the only diffusion calculations. Several computational methods have been proposed for power or flux mapping. The CANDU on-line flux mapping system (Tang et al., 1978) converts the 102 vanadium detector signals to thermal fluxes at the detector sites and then maps out the 3D flux distribution by a process of least-squares fitting of the measured
⇑ Corresponding author at: Department of Engineering Physics, Tsinghua University, Beijing, China. E-mail addresses: [email protected], [email protected] (X. Peng). 0306-4549/$ - see front matter Ó 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2014.01.002
thermal fluxes to a linear expansion of pre-calculated flux modes. Combustion Engineering (CE) nuclear power plants use the Combustion Engineering Core (CECOR) Karlson, 1995 method to estimate the power distributions, and the pre-calculated twodimensional coupling coefficients are used. Jang et al. (2004) proposed a three-dimensional coupling coefficients method and Webb and Brittingham (2000) proposed a Lagrange multiplier method, which both can be regarded as an improved version of the CECOR method. Lee and Kim (2003) proposed a least-squares method by combining the coarse mesh finite difference (CMFD) form of the fixed-source diffusion equation and the detector response equation to form an over-determined linear equation. Zhong et al. (2010) proposed a mixed method based on the least-squares method and the harmonics synthesis method. Li et al. (2013) proposed three methods, namely, weight coefficient method, polynomial expand method and thin plane spline method to fit the deviation between measured and predicted results for two-dimensional radial plane. The quality of maps obtained by interpolation of observations of a target spatial variable at a restricted number of locations, is partly determined by the spatial pattern of the detector locations. In order to get the best possible results from the power mapping reconstruction procedure, the detector locations need to be optimized. The algorithm for optimal placement of detectors appear to be very similar regardless of the application. All can be posed as selecting a subset of locations from a large set of candidate locations. A natural approach to solve this problem is to minimize the prediction error variance by using effective combinatorial
119
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
optimization methods like genetic algorithm and simulated annealing algorithm. The ordinary kriging method of interpolation in geostatistical analysis is applied to many different subjects (Pokhrel et al., 2013), such as precipitation, elevation, air temperature, soil properties and electrical conductivities. In nuclear field, Lockwood and Anitescu (2012) investigate the issue of providing a statistical model for the response of a computer model of nuclear engineering system for use in uncertainty propagation. In this paper, we introduce the ordinary kriging method into nuclear reactor power distribution monitoring. A new power mapping method based on ordinary kriging method is presented, and then, the spatially averaged ordinary kriging variance is minimized using the well known simulated annealing algorithm (Brus et al., 2013) based on a Metropolis–Hastings one to obtain the optimal detector locations.
Varðr j Þ ¼ Cov fv~ j v~ j g 2Cov fv~ j v^ j g þ Cov fv^ j v^ j g
ð5Þ
The first term can further be simplified as follow,
Cov fv~ j v~ j g ¼ Varfv~ j g ¼ Var
( ) N X wi v i i¼1
¼
N X N X
wi wk Cov fv i v k g ¼
i¼1 k¼1
N X N X b ik wi wk C
ð6Þ
i¼1 k¼1
The second term can be written as,
(
! ) N N X X wi v i v^ j ¼ wi Cov fv i v^ j g
Cov fv~ j v^ j g ¼ Cov
i¼1
i¼1
N X b i0 wi C ¼
ð7Þ
i¼1
2. Methodology
Finally, the third term can be expressed as
2.1. Ordinary kriging method
b 00 Cov fv^ j v^ j g ¼ C
ð8Þ
Substituting from Eqs. (6)–(8) in Eq. (5), Kriging is a group of geostatistical techniques (Pokhrel et al., 2013) to interpolate the value of a random field at an unobserved location from observations of its value at nearby locations. It belongs to a family of linear least squares estimation algorithms that are used in several geostatistical applications. Linear kriging is divided into simple kriging (known mean), ordinary kriging (unknown but constant mean) and universal kriging (the mean is an unknown linear combination of known functions), depending on the mean value specification. We shall restrict the discussion here to ordinary kriging. Ordinary kriging is widely used because it is statistically the best linear unbiased estimator. Ordinary kriging is linear because its estimates are linear combination of the available data. It is unbiased because it attempts to keep the mean residual to be zero. Finally, it is called best because it tries to minimize the residual variance. Let the data be sampled at N locations (x1, . . . , xN), and the corre^ j at an unknown location sponding values be (v1, . . . , vN). The value v ^ xj is estimated as a weighted linear combination of sampled values, given by, N X v~ j ¼ wi v i
ð1Þ
i¼1
Here, v~ j is the estimate, and let v^ j be the actual value (unknown) at ^xj . To find the weights, the values (vi and v ^ j ) are assumed to be stationary random functions,
E½v i ¼ E½v^ j ¼ Eðv Þ
ð2Þ
Varðr j Þ ¼
i¼1 k¼1
E½v^ j E½v~ j ¼ 0 " # N X wi v i ¼ 0 E½v^ j E i¼1
N X N N N X X X b ik 2 wi C b i0 þ C b 00 þ 2k JðwÞ ¼ wi wk C wi 1 i¼1 k¼1
E½v^ j
N X @J b ik 2 C b i0 þ 2k ¼ 2 wk C @wi k¼1 ! N X @J ¼2 wi 1 @k i¼1
ð11Þ ð12Þ
Setting this to zero, we get the following system to solve, to obtain the weights w and k,
0
b 11 C B B .. B . B Bb @ C N1 1
b 1N C .. . b NN C
1
.. .
10 1 0b 1 w1 1 C 10 C C C B . .. CB B B .. C B ... C C .C C B ¼ CB C B C@ w C Bb C A 1A N @ C N0 A k 0 1
ð13Þ
Substituting Eq. (13) into Eq. (9), we can derive the kriging variance, N X N N X X b ik 2 wi C b i0 þ C b 00 wi wk C i¼1 k¼1
i¼1
i¼1
¼
! b 00 þC
k¼1
N X b i0 k 2 C b i0 þ C b 00 wi C i¼1
b 00 ¼C
i¼1 N X wi ¼ 1
N N N X X X b i0 b 00 b i0 k wi C wi k ¼ C wi C i¼1
i¼1
ð14Þ
i¼1
Eq. (13) needs to be solved at each unknown location ^ xj .
i¼1
Let the residue be rj
Therefore, the residual variance is given by,
ð10Þ
i¼1
N N X X b ik 2 C b i0 ¼ wi wk C
N X wi E½v ¼ 0
r j ¼ v~ j v^ j
i¼1
!
with 2k the Lagrange multiplier. Taking derivatives of J with respect to w and k,
Varðr j Þ ¼
ð3Þ
ð9Þ
i¼1
For ordinary kriging, it is required to find w by minimizing Var(rj) with respect to w subject to the constraint Eq. (3). This can be written as the minimization of the penalized cost function,
For unbiased estimates, we can get
E½v i ¼ E½v^ j ¼ Eðv Þ E½v^ j v~ j ¼ 0
N X N N X X b ik 2 wi C b i0 þ C b 00 wi wk C
2.2. Covariance functions
ð4Þ
b ij , either by a There are two ways of specifying the covariance C standard function or by evaluating it empirically at each location. A
120
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
functional form of covariance is preferred in this study. The covariance function is generally chosen to reflect prior information. The second order auto-regressive function (Bertrand et al., 2012) is used to prescribe the covariance. In such a function, the amount of covariance depends from the Euclidean distance between spatial points. The correlation length L has different values, which means we are dealing with a global pseudo Euclidean distance. The used function can be expressed as
b ij ¼ r2 1 þ kxi xj k exp kxi xj k C L L
ð15Þ
where r is the standard deviation, and there is a assumption that all the random variables have the same standard deviation. The value of r does not have an effect on the weights w and the interpolation result, but it does have a direct effect on the kriging variance, i.e. the uncertainty of the interpolation result. 2.3. Application in power mapping In this study, we just use the ordinary kriging method to reconstruct the radial power distribution. The axial power distribution can be reconstructed by the cubic spline synthesis method (Wang et al., 1991). The measured power of instrumented assembly i is obtained from
Pmea ¼ Ii W i i
ð16Þ
P mea i
where is the measured power, Ii is the integrated detector signal, Wi is the power-to-signal ratio, and i labels the location of the instrumented assembly. The power-to-signal ratio Wi can be calculated from fine-mesh, two-dimensional, multigroup diffusion theory (Webb and Brittingham, 2000). The model calibration factor is defined by
Di ¼ Pmea =Pcal i i
ð17Þ
Pcal i
where is the predicted power distribution calculated by neutronics design code SMART in SCIENCE code package (Zheng, 2003), Pmea is the measured power distribution provided by the nui clear power plant, and i labels the location of the instrumented assembly. The model calibration factors are defined in the instrumented assemblies. The calibration factors of the uninstrumented assemblies are obtained using the ordinary kriging interpolation technique. Then the best estimate power distribution, i.e. reconstructed power distribution, can be obtained:
e j Pcal Precon ¼D j j
where is the reconstructed power distribution of the uninstre j is the calibration factor of the uninstruumented assemblies, D mented assemblies, and Pcal is the predicted power distribution j calculated by neutronics design code SMART, and j labels the location of the uninstrumented assembly. 2.4. CECOR method descriptions The purpose of this paper is to discuss the research investigating the ordinary kriging method for estimating power, and how it compares to the CECOR method, so some brief descriptions of the CECOR are needed. The CECOR method (Jang et al., 2004) determine the powers in the uninstrumented assemblies through the use of radial coupling coefficients. The average coupling coefficient of an assembly is the ratio of the average power in the surrounding boxes to the power in the assembly itself:
PNi
cal j¼1 P j ðiÞ cal Ni P ðiÞ
hCCiðiÞ Ni PðiÞ
ð19Þ
X X Pj ðiÞ ¼ Pj ðiÞ j2U
ð20Þ
j2I
where U represents the group of uninstrumented assemblies abut on the assembly i, and I is the opposite instrumented group. Eq. (20) is applied to all the uninstrumented assemblies and can be expressed as a matrix–vector form:
AP ¼ s
ð21Þ
If we solve Eq. (21), we can obtain all the uninstrumented assembly powers throughout the whole core. 2.5. Calculating the optimal detector locations From Eq. (14), it shows that the kriging variance depends only on the distances between spatial points, the values of the observations have no influence. In other words, if the covariance function is known, the kriging variance can be predicted for any proposed detector location strategy prior to the actual survey. Kriging is, therefore, an ideal tool for designing optimal detector location strategies. The objective is to find the optimal detector location that has the smallest AOKV of unknown locations. We consider the case where the number of detectors and allowed positions are fixed. The only freedom we have is the choice of locations of detectors. Now, we want to obtain a set of optimal detector locations from allowed positions. We could introduce the parameter qj, j = 1, . . . , M such that
qj ¼
1j2I 0j2U
M X qj ¼ N
ð22Þ
j¼1
where N is the number of detectors, U represents the group of uninstrumented assemblies, and I is the opposite instrumented group. The AOKV of the unknown locations can be expressed as:
gðqÞ ¼
ð18Þ
P recon j
hCCiðiÞ ¼
where assembly j represents the assembly abut on the assembly i, and Ni represents the total number of the surrounding i assemblies. The coefficients are precalculated using power distributions obtained from neutronics design code SMART. After obtaining the coupling coefficients, the power of the uninstrumented assemblies can be solved by Eq. (20):
1 X Varðr j Þ M N j2U
ð23Þ
Then the optimal detector location problem can be expressed as:
min gðqÞ 1j2I qj ¼ 0j2U
M X
qj ¼ N
ð24Þ
j¼1
For practical applications, the number of combinations of qj will be formidable which means that an exhaustive search over all possible patterns is prohibitive. For example, the DayaBay Unit 1 reactor has 157 assemblies and 50 detectors, then the optimal detector location problem is to find the optimal location of 50 instruments within the 157 allowed positions. The pure brute force approach is impossible to solve this problem and some efficient search algorithm has to be employed. Here we use simulated annealing algorithm. Simulated annealing (Brus et al., 2013) is an iterative, combinatorial optimization algorithm in which a sequence of combinations is generated by deriving a new combination from slightly and randomly changing the previous combination. Each time a new combination is generated, the AOKV is evaluated and compared
121
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
with the value of the previous combination. The new combination is accepted if the AOKV has decreased by the change. However, the annealing algorithm also accepts some of the changes that increase the AOKV. This is to avoid being trapped in a local optimum. The probability of accepting a worse combination is given by
Df P ¼ exp T
ð25Þ
where Df is the increase in the AOKV and T is a control parameter which, by analogy with the original application of simulating the cooling of a metal into a minimum energy crystalline structure, is known as the system temperature. This temperature is gradually decreased during the optimization. Eq. (25) shows that, given the temperature, the larger the increase of the AOKV, the smaller the probability of acceptance. The temperature remains constant during a fixed number of transitions and is then decreased, which makes that the acceptance probability gradually decreases as the iteration continues. We choose the classical simple inverse cooling scheme as follows,
Ti ¼
T0 iþ1
ð26Þ
where i is the iteration step, T0 is the initial temperature. Besides these parameters, there are mainly 2 parameters to set in the present case. The first one is the number of transitions for a given temperature, which can be chosen as fixed. The second one is the stopping criterion, which may simply be a fixed number of total iterations.
The assemblies without detector are marked in white, and the assemblies with detector are marked in blue.1 The reconstruction error can be defined as following:
3. Numerical results
Error ¼
3.1. Power distribution reconstruction
where Preconstruction represents the reconstructed power distribution, and Preference represents the reference power distribution which is provided by PUISSANCE ESTIMEE code in SCIENCE code package. PUISSANCE ESTIMEE code uses the measurements of 50 detectors to give a two-dimensional whole core power distribution. The comparison of the 1/4 reactor power distribution reconstruction error of different power mapping method is represented in Fig. 2. The result shows that the ordinary kriging method provides credible results and performs better than the CECOR method. The maximum relative error and the root mean square error are defined as following:
Eighteen measurements of different burn-up cases of DayaBay Unit 1 reactor was used to verify the accuracy of the ordinary kriging power mapping method, and comparison was made between this method and the CECOR method (Karlson, 1995). DayaBay Unit 1 reactor has 157 assemblies and 50 detectors as shown in Fig. 1.
Fig. 2. Relative reconstruction error distribution (burn-up is 2.576 GWd/tU).
RMS ¼
Preconstruction Preference P reference
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 uP157 Pireconstruction Pi reference u t i¼1 P ireference
157 P i P ireference MRE ¼ max reconstruction ði ¼ 1; . . . ; 157Þ i Pi
ð27Þ
ð28Þ ð29Þ
reference
In Fig. 3, the maximum relative errors (Maximum RE) between the reconstruction and the measurement results with different burn-up are presented. The root mean square (RMS) errors of the diffusion calculation and reconstruction results are also shown in the same figure. From Fig. 3, we can see that RMS errors of the reconstruction results are lower than those of the diffusion calculation. Compared with the diffusion results, the RMS reconstruction errors are more stable, and are kept at less than 0.35%, and the maximum RE reconstruction errors are kept at less than 1.02%. This proves that the ordinary kriging method is applicable for the entire operating cycle. The ordinary kriging provide a natural way to estimate the
Fig. 1. Current detector location of Unit 1 reactor of DayaBay Nuclear Power Plant.
1 For interpretation of color in Figs. 1 and 6, the reader is referred to the web version of this article.
122
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
6 5
0.8
4
0.6
Errors (%)
Errors (%)
1
Maximum RE of Diffusion calculation RMS errors of Diffusion calculation Maximum RE of Reconstruction RMS errors of Reconstruction
3 2 1
0.4 0.2 0
0 0
2
4
6
8
10
12
14
16
18
20
-0.2
Burn-up (GWd/tU) -0.4 0
Fig. 3. Relative errors with different burn-up.
5
10
15
20
25
30
35
40
45
50
Instrumented assembly number
uncertainty of the interpolation result, and the kriging variance of the interpolation result is proportional to the variance r2. If we use the model calibration factor to interpolate, then the variance r2 includes variance due to error, i.e. the residuals between calculation and measurement. There are 107 assemblies without detector, and Fig. 4 shows the standard deviation (SD) of the interpolation result. The value of r is usually under 3%, so most of the standard deviations of different interpolation results can be less than 1%. If there is less errors between calculation and measurement, the interpolation results will have less uncertainties. To perform a more direct verification of our method, we use a subset of the instrumented assembly power to predict the other instrumented assembly power. Here, the ‘‘leave-one-out’’ validation strategy is used. We suppose that the instruments are failed sequentially; the powers are then calculated for that failed instrument, i.e., treated as non-instrumented; and last, restored as an instrumented assembly. In this case study, only one instrument is failed at a time. The reconstruction errors of different failed instrument are represented in Fig. 5. It shows that the ordinary kriging method can reconstruct the failed instrument power accurately. The assemblies with large reconstruction error all are located in the periphery of the core. High calculation speed is a basic requirement for on-line power distribution monitoring systems. The weights wi of the ordinary kriging interpolation can be calculated off-line, so the power mapping procedure is fast and meet the requirement of on-line monitoring.
Fig. 5. Reconstruction errors of failed instrumented assembly.
3.2. Neutron detector location optimization In-core detector optimal locations searching method based on ordinary kriging and simulated annealing was tested. Objective function, i.e. the average ordinary kriging covariance corresponding to the detector locations was formulated by the ordinary kriging, and its minimum was searched by simulated annealing algorithm, and then the optimal detector locations were obtained as shown in Fig. 6, the red marked assemblies only belong to the current detector location strategy, the black marked assemblies only belong to the optimal detector location strategy, and the blue marked assemblies belong to both location strategies. The number of transitions for a given temperature was chosen as 5, the number of the outer iterations was chosen as 3000, and the initial temperature was chosen as 100. A random detector locations were generated as shown in Fig. 7, the black marked assemblies only belong to the random detector location strategy. The comparison of the AOKV of three location patterns were shown in Table 1, where r is the standard deviation in Eq. (15).
0.35
0.3
SD/
0.25
0.2
0.15
0.1 0
20
40
60
80
100
120
Uninstrumented assembly number Fig. 4. Standard deviation of interpolation result.
Fig. 6. Comparison of the optimal detector location with the current detector location.
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
123
(1) The reconstruction results show that the ordinary kriging method provides credible results and performs better than the CECOR method. The RMS reconstruction errors are kept at less than 0.35%, and the maximum RE reconstruction errors are kept at less than 1.02% for the entire operating cycle. (2) The ordinary kriging method for power distribution monitoring is fast enough for on-line monitoring. (3) The ordinary kriging method can be used together with combinatorial optimization algorithm, such as simulated annealing, to determine the optimal neutron detector locations. More testing needs to be done to ensure that the ordinary kriging performs good enough in terms of power distribution reconstruction. And also, there is a potential for expanding the ordinary kriging method to a three-dimensional approach that could produce even better results. References
Fig. 7. Comparison of the random detector location with the current detector location.
Table 1 Comparison of the AOKVs. Locations pattern
AOKV (r2)
Current locations Random locations Optimal locations
0.0319 0.0510 0.0309
From Table 1, we can see that the current detector locations of Unit 1 reactor of DayaBay Nuclear Power Plant has a near optimal AOKV, and the random detector locations cannot give a good AOKV. 4. Conclusion and future work This paper presents the concept of using the ordinary kriging method to reconstruct reactor power distributions and determine the optimal neutron detector locations. The DayaBay Unit 1 PWR is taken as an example to verify this method. The following conclusions are drawn from the study:
Bertrand, B. et al., 2012. Robustness of nuclear core activity reconstruction by data assimilation. J. Power Eng. Syst. 6 (2), 289–301. Boyd, A., Miller, W., 1996. The BEACON on-line core monitoring system: functional upgrades and applications. In: In-Core Instrumentation and Reactor Core Assessment. Mito-shi, Japan. Brus, D. et al., 2013. Optimization of sample patterns for universal kriging of environmental variables. Geoderma 138 (1–2), 86–95. Jang, J. et al., 2004. On-line generation of three-dimensional core power distribution using in-core detector signals to monitor safety limits. J. Korean Nucl. Soc. 36 (6), 528–539. Karlson, C., 1995. Continuing advancements in in-core power distribution measurement methods using SIMULATE-3 and CECOR-3.4. Nucl. Sci. Eng. 121, 57–66. Lee, K., Kim, C., 2003. The least-squares method for three-dimensional core power distribution monitoring in pressurized water reactors. Nucl. Sci. Eng. 143 (3), 268–280. Li, W., et al., 2013. Development of a new flux map processing code for moveable detector system in PWR. In: (M&C 2013). Sun Valley, America. Lockwood, B., Anitescu, M., 2012. Gradient-enhanced universal kriging for uncertainty propagation. Nucl. Sci. Eng. 170, 168–195. Pokhrel, R. et al., 2013. A kriging method of interpolation used to map liquefaction potential over alluvial ground. Eng. Geol. 152 (1), 26–37. Tang, T., et al., 1978. Analytical Design of the CANDU-600 On-line Flux Mapping System. TDAI-152, Atomic Energy of Canada Limited. Tojo, M., et al., 2008. Application of the new core monitoring system ‘GNF-ARGOS’ to the start-up tracking calculation of ABWR. In: International Conference on the Physics of Reactors. Interlaken, Switzerland. Wang, K. et al., 1991. Application of cubic spline synthesis in on-line core axial power distribution monitoring. J. Korean Nucl. Soc. 23 (3), 316–320. Webb, R., Brittingham, J., 2000. Comparison of CECOR algorithm to Lagrange multiplier method to estimate reactor power distributions. Nucl. Technol. 132, 206–213. Zheng, S., 2003. GNPS 1&2 18 Month Fuel Cycles Reload Study Procedure. TFJN/DC/ 2322, Framatome ANP. Zhong, B., et al., 2010. A new method on flux mapping by function expansion. In: International Conference on Nuclear Engineering (ICONE18). Xi’an, China.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
A new power mapping method based on ordinary kriging and determination of optimal detector location strategy Xingjie Peng a,b,⇑, Kan Wang a, Qing Li b a b
Department of Engineering Physics, Tsinghua University, Beijing, China Science and Technology on Reactor System Design Technology Laboratory, Nuclear Power Institute of China, Chengdu, China
a r t i c l e
i n f o
Article history: Received 16 September 2013 Received in revised form 27 December 2013 Accepted 3 January 2014 Available online 2 February 2014 Keywords: Power mapping Ordinary kriging Optimal detector location Simulated annealing
a b s t r a c t The Ordinary Kriging (OK) method is presented that is designed for a core power mapping calculation of pressurized water reactors (PWRs). Measurements from DayaBay Unit 1 PWR are used to verify the accuracy of the OK method. The root mean square (RMS) reconstruction errors are kept at less than 0.35%, and the maximum reconstruction relative errors (RE) are kept at less than 1.02% for the entire operating cycle. The reconstructed assembly power distribution results show that the OK method is fit for core power distribution monitoring. The quality of power distribution obtained by the OK method is partly determined by the neutron detector locations, and the OK method is also applied to solve the optimal neutron detector location problem. The spatially averaged ordinary kriging variance (AOKV) is minimized using simulated annealing, and then, the optimal in-core neutron detector locations are obtained. The result shows that the current neutron detector location of DayaBay Unit 1 reactor is near-optimal. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Core power distribution monitoring in operating power reactors is very important in core surveillance, the 3-D power distribution is one of the basic operation parameters which can determine many other important parameters such as power peaking factor, enthalpy rising factor and quadrant tilt ratio used to evaluate the operation condition of reactor and the safe margin. The economy of reactor could be optimized if the real time 3-D power distribution is well obtained and used for surveillance and regulation. Most commercial power reactors in operation are equipped with fixed or movable in-core neutron detectors to obtain power distribution information. Many kinds of on-line monitoring systems, such as BEACON (Boyd and Miller, 1996) and GNF-ARGOS (Tojo et al., 2008), have been developed to estimate in-core power distributions using fixed in-core detectors. The detector results at certain locations reflect the actual reactor flux or power can be applied to improve the results of the only diffusion calculations. Several computational methods have been proposed for power or flux mapping. The CANDU on-line flux mapping system (Tang et al., 1978) converts the 102 vanadium detector signals to thermal fluxes at the detector sites and then maps out the 3D flux distribution by a process of least-squares fitting of the measured
⇑ Corresponding author at: Department of Engineering Physics, Tsinghua University, Beijing, China. E-mail addresses: [email protected], [email protected] (X. Peng). 0306-4549/$ - see front matter Ó 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.anucene.2014.01.002
thermal fluxes to a linear expansion of pre-calculated flux modes. Combustion Engineering (CE) nuclear power plants use the Combustion Engineering Core (CECOR) Karlson, 1995 method to estimate the power distributions, and the pre-calculated twodimensional coupling coefficients are used. Jang et al. (2004) proposed a three-dimensional coupling coefficients method and Webb and Brittingham (2000) proposed a Lagrange multiplier method, which both can be regarded as an improved version of the CECOR method. Lee and Kim (2003) proposed a least-squares method by combining the coarse mesh finite difference (CMFD) form of the fixed-source diffusion equation and the detector response equation to form an over-determined linear equation. Zhong et al. (2010) proposed a mixed method based on the least-squares method and the harmonics synthesis method. Li et al. (2013) proposed three methods, namely, weight coefficient method, polynomial expand method and thin plane spline method to fit the deviation between measured and predicted results for two-dimensional radial plane. The quality of maps obtained by interpolation of observations of a target spatial variable at a restricted number of locations, is partly determined by the spatial pattern of the detector locations. In order to get the best possible results from the power mapping reconstruction procedure, the detector locations need to be optimized. The algorithm for optimal placement of detectors appear to be very similar regardless of the application. All can be posed as selecting a subset of locations from a large set of candidate locations. A natural approach to solve this problem is to minimize the prediction error variance by using effective combinatorial
119
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
optimization methods like genetic algorithm and simulated annealing algorithm. The ordinary kriging method of interpolation in geostatistical analysis is applied to many different subjects (Pokhrel et al., 2013), such as precipitation, elevation, air temperature, soil properties and electrical conductivities. In nuclear field, Lockwood and Anitescu (2012) investigate the issue of providing a statistical model for the response of a computer model of nuclear engineering system for use in uncertainty propagation. In this paper, we introduce the ordinary kriging method into nuclear reactor power distribution monitoring. A new power mapping method based on ordinary kriging method is presented, and then, the spatially averaged ordinary kriging variance is minimized using the well known simulated annealing algorithm (Brus et al., 2013) based on a Metropolis–Hastings one to obtain the optimal detector locations.
Varðr j Þ ¼ Cov fv~ j v~ j g 2Cov fv~ j v^ j g þ Cov fv^ j v^ j g
ð5Þ
The first term can further be simplified as follow,
Cov fv~ j v~ j g ¼ Varfv~ j g ¼ Var
( ) N X wi v i i¼1
¼
N X N X
wi wk Cov fv i v k g ¼
i¼1 k¼1
N X N X b ik wi wk C
ð6Þ
i¼1 k¼1
The second term can be written as,
(
! ) N N X X wi v i v^ j ¼ wi Cov fv i v^ j g
Cov fv~ j v^ j g ¼ Cov
i¼1
i¼1
N X b i0 wi C ¼
ð7Þ
i¼1
2. Methodology
Finally, the third term can be expressed as
2.1. Ordinary kriging method
b 00 Cov fv^ j v^ j g ¼ C
ð8Þ
Substituting from Eqs. (6)–(8) in Eq. (5), Kriging is a group of geostatistical techniques (Pokhrel et al., 2013) to interpolate the value of a random field at an unobserved location from observations of its value at nearby locations. It belongs to a family of linear least squares estimation algorithms that are used in several geostatistical applications. Linear kriging is divided into simple kriging (known mean), ordinary kriging (unknown but constant mean) and universal kriging (the mean is an unknown linear combination of known functions), depending on the mean value specification. We shall restrict the discussion here to ordinary kriging. Ordinary kriging is widely used because it is statistically the best linear unbiased estimator. Ordinary kriging is linear because its estimates are linear combination of the available data. It is unbiased because it attempts to keep the mean residual to be zero. Finally, it is called best because it tries to minimize the residual variance. Let the data be sampled at N locations (x1, . . . , xN), and the corre^ j at an unknown location sponding values be (v1, . . . , vN). The value v ^ xj is estimated as a weighted linear combination of sampled values, given by, N X v~ j ¼ wi v i
ð1Þ
i¼1
Here, v~ j is the estimate, and let v^ j be the actual value (unknown) at ^xj . To find the weights, the values (vi and v ^ j ) are assumed to be stationary random functions,
E½v i ¼ E½v^ j ¼ Eðv Þ
ð2Þ
Varðr j Þ ¼
i¼1 k¼1
E½v^ j E½v~ j ¼ 0 " # N X wi v i ¼ 0 E½v^ j E i¼1
N X N N N X X X b ik 2 wi C b i0 þ C b 00 þ 2k JðwÞ ¼ wi wk C wi 1 i¼1 k¼1
E½v^ j
N X @J b ik 2 C b i0 þ 2k ¼ 2 wk C @wi k¼1 ! N X @J ¼2 wi 1 @k i¼1
ð11Þ ð12Þ
Setting this to zero, we get the following system to solve, to obtain the weights w and k,
0
b 11 C B B .. B . B Bb @ C N1 1
b 1N C .. . b NN C
1
.. .
10 1 0b 1 w1 1 C 10 C C C B . .. CB B B .. C B ... C C .C C B ¼ CB C B C@ w C Bb C A 1A N @ C N0 A k 0 1
ð13Þ
Substituting Eq. (13) into Eq. (9), we can derive the kriging variance, N X N N X X b ik 2 wi C b i0 þ C b 00 wi wk C i¼1 k¼1
i¼1
i¼1
¼
! b 00 þC
k¼1
N X b i0 k 2 C b i0 þ C b 00 wi C i¼1
b 00 ¼C
i¼1 N X wi ¼ 1
N N N X X X b i0 b 00 b i0 k wi C wi k ¼ C wi C i¼1
i¼1
ð14Þ
i¼1
Eq. (13) needs to be solved at each unknown location ^ xj .
i¼1
Let the residue be rj
Therefore, the residual variance is given by,
ð10Þ
i¼1
N N X X b ik 2 C b i0 ¼ wi wk C
N X wi E½v ¼ 0
r j ¼ v~ j v^ j
i¼1
!
with 2k the Lagrange multiplier. Taking derivatives of J with respect to w and k,
Varðr j Þ ¼
ð3Þ
ð9Þ
i¼1
For ordinary kriging, it is required to find w by minimizing Var(rj) with respect to w subject to the constraint Eq. (3). This can be written as the minimization of the penalized cost function,
For unbiased estimates, we can get
E½v i ¼ E½v^ j ¼ Eðv Þ E½v^ j v~ j ¼ 0
N X N N X X b ik 2 wi C b i0 þ C b 00 wi wk C
2.2. Covariance functions
ð4Þ
b ij , either by a There are two ways of specifying the covariance C standard function or by evaluating it empirically at each location. A
120
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
functional form of covariance is preferred in this study. The covariance function is generally chosen to reflect prior information. The second order auto-regressive function (Bertrand et al., 2012) is used to prescribe the covariance. In such a function, the amount of covariance depends from the Euclidean distance between spatial points. The correlation length L has different values, which means we are dealing with a global pseudo Euclidean distance. The used function can be expressed as
b ij ¼ r2 1 þ kxi xj k exp kxi xj k C L L
ð15Þ
where r is the standard deviation, and there is a assumption that all the random variables have the same standard deviation. The value of r does not have an effect on the weights w and the interpolation result, but it does have a direct effect on the kriging variance, i.e. the uncertainty of the interpolation result. 2.3. Application in power mapping In this study, we just use the ordinary kriging method to reconstruct the radial power distribution. The axial power distribution can be reconstructed by the cubic spline synthesis method (Wang et al., 1991). The measured power of instrumented assembly i is obtained from
Pmea ¼ Ii W i i
ð16Þ
P mea i
where is the measured power, Ii is the integrated detector signal, Wi is the power-to-signal ratio, and i labels the location of the instrumented assembly. The power-to-signal ratio Wi can be calculated from fine-mesh, two-dimensional, multigroup diffusion theory (Webb and Brittingham, 2000). The model calibration factor is defined by
Di ¼ Pmea =Pcal i i
ð17Þ
Pcal i
where is the predicted power distribution calculated by neutronics design code SMART in SCIENCE code package (Zheng, 2003), Pmea is the measured power distribution provided by the nui clear power plant, and i labels the location of the instrumented assembly. The model calibration factors are defined in the instrumented assemblies. The calibration factors of the uninstrumented assemblies are obtained using the ordinary kriging interpolation technique. Then the best estimate power distribution, i.e. reconstructed power distribution, can be obtained:
e j Pcal Precon ¼D j j
where is the reconstructed power distribution of the uninstre j is the calibration factor of the uninstruumented assemblies, D mented assemblies, and Pcal is the predicted power distribution j calculated by neutronics design code SMART, and j labels the location of the uninstrumented assembly. 2.4. CECOR method descriptions The purpose of this paper is to discuss the research investigating the ordinary kriging method for estimating power, and how it compares to the CECOR method, so some brief descriptions of the CECOR are needed. The CECOR method (Jang et al., 2004) determine the powers in the uninstrumented assemblies through the use of radial coupling coefficients. The average coupling coefficient of an assembly is the ratio of the average power in the surrounding boxes to the power in the assembly itself:
PNi
cal j¼1 P j ðiÞ cal Ni P ðiÞ
hCCiðiÞ Ni PðiÞ
ð19Þ
X X Pj ðiÞ ¼ Pj ðiÞ j2U
ð20Þ
j2I
where U represents the group of uninstrumented assemblies abut on the assembly i, and I is the opposite instrumented group. Eq. (20) is applied to all the uninstrumented assemblies and can be expressed as a matrix–vector form:
AP ¼ s
ð21Þ
If we solve Eq. (21), we can obtain all the uninstrumented assembly powers throughout the whole core. 2.5. Calculating the optimal detector locations From Eq. (14), it shows that the kriging variance depends only on the distances between spatial points, the values of the observations have no influence. In other words, if the covariance function is known, the kriging variance can be predicted for any proposed detector location strategy prior to the actual survey. Kriging is, therefore, an ideal tool for designing optimal detector location strategies. The objective is to find the optimal detector location that has the smallest AOKV of unknown locations. We consider the case where the number of detectors and allowed positions are fixed. The only freedom we have is the choice of locations of detectors. Now, we want to obtain a set of optimal detector locations from allowed positions. We could introduce the parameter qj, j = 1, . . . , M such that
qj ¼
1j2I 0j2U
M X qj ¼ N
ð22Þ
j¼1
where N is the number of detectors, U represents the group of uninstrumented assemblies, and I is the opposite instrumented group. The AOKV of the unknown locations can be expressed as:
gðqÞ ¼
ð18Þ
P recon j
hCCiðiÞ ¼
where assembly j represents the assembly abut on the assembly i, and Ni represents the total number of the surrounding i assemblies. The coefficients are precalculated using power distributions obtained from neutronics design code SMART. After obtaining the coupling coefficients, the power of the uninstrumented assemblies can be solved by Eq. (20):
1 X Varðr j Þ M N j2U
ð23Þ
Then the optimal detector location problem can be expressed as:
min gðqÞ 1j2I qj ¼ 0j2U
M X
qj ¼ N
ð24Þ
j¼1
For practical applications, the number of combinations of qj will be formidable which means that an exhaustive search over all possible patterns is prohibitive. For example, the DayaBay Unit 1 reactor has 157 assemblies and 50 detectors, then the optimal detector location problem is to find the optimal location of 50 instruments within the 157 allowed positions. The pure brute force approach is impossible to solve this problem and some efficient search algorithm has to be employed. Here we use simulated annealing algorithm. Simulated annealing (Brus et al., 2013) is an iterative, combinatorial optimization algorithm in which a sequence of combinations is generated by deriving a new combination from slightly and randomly changing the previous combination. Each time a new combination is generated, the AOKV is evaluated and compared
121
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
with the value of the previous combination. The new combination is accepted if the AOKV has decreased by the change. However, the annealing algorithm also accepts some of the changes that increase the AOKV. This is to avoid being trapped in a local optimum. The probability of accepting a worse combination is given by
Df P ¼ exp T
ð25Þ
where Df is the increase in the AOKV and T is a control parameter which, by analogy with the original application of simulating the cooling of a metal into a minimum energy crystalline structure, is known as the system temperature. This temperature is gradually decreased during the optimization. Eq. (25) shows that, given the temperature, the larger the increase of the AOKV, the smaller the probability of acceptance. The temperature remains constant during a fixed number of transitions and is then decreased, which makes that the acceptance probability gradually decreases as the iteration continues. We choose the classical simple inverse cooling scheme as follows,
Ti ¼
T0 iþ1
ð26Þ
where i is the iteration step, T0 is the initial temperature. Besides these parameters, there are mainly 2 parameters to set in the present case. The first one is the number of transitions for a given temperature, which can be chosen as fixed. The second one is the stopping criterion, which may simply be a fixed number of total iterations.
The assemblies without detector are marked in white, and the assemblies with detector are marked in blue.1 The reconstruction error can be defined as following:
3. Numerical results
Error ¼
3.1. Power distribution reconstruction
where Preconstruction represents the reconstructed power distribution, and Preference represents the reference power distribution which is provided by PUISSANCE ESTIMEE code in SCIENCE code package. PUISSANCE ESTIMEE code uses the measurements of 50 detectors to give a two-dimensional whole core power distribution. The comparison of the 1/4 reactor power distribution reconstruction error of different power mapping method is represented in Fig. 2. The result shows that the ordinary kriging method provides credible results and performs better than the CECOR method. The maximum relative error and the root mean square error are defined as following:
Eighteen measurements of different burn-up cases of DayaBay Unit 1 reactor was used to verify the accuracy of the ordinary kriging power mapping method, and comparison was made between this method and the CECOR method (Karlson, 1995). DayaBay Unit 1 reactor has 157 assemblies and 50 detectors as shown in Fig. 1.
Fig. 2. Relative reconstruction error distribution (burn-up is 2.576 GWd/tU).
RMS ¼
Preconstruction Preference P reference
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 uP157 Pireconstruction Pi reference u t i¼1 P ireference
157 P i P ireference MRE ¼ max reconstruction ði ¼ 1; . . . ; 157Þ i Pi
ð27Þ
ð28Þ ð29Þ
reference
In Fig. 3, the maximum relative errors (Maximum RE) between the reconstruction and the measurement results with different burn-up are presented. The root mean square (RMS) errors of the diffusion calculation and reconstruction results are also shown in the same figure. From Fig. 3, we can see that RMS errors of the reconstruction results are lower than those of the diffusion calculation. Compared with the diffusion results, the RMS reconstruction errors are more stable, and are kept at less than 0.35%, and the maximum RE reconstruction errors are kept at less than 1.02%. This proves that the ordinary kriging method is applicable for the entire operating cycle. The ordinary kriging provide a natural way to estimate the
Fig. 1. Current detector location of Unit 1 reactor of DayaBay Nuclear Power Plant.
1 For interpretation of color in Figs. 1 and 6, the reader is referred to the web version of this article.
122
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
6 5
0.8
4
0.6
Errors (%)
Errors (%)
1
Maximum RE of Diffusion calculation RMS errors of Diffusion calculation Maximum RE of Reconstruction RMS errors of Reconstruction
3 2 1
0.4 0.2 0
0 0
2
4
6
8
10
12
14
16
18
20
-0.2
Burn-up (GWd/tU) -0.4 0
Fig. 3. Relative errors with different burn-up.
5
10
15
20
25
30
35
40
45
50
Instrumented assembly number
uncertainty of the interpolation result, and the kriging variance of the interpolation result is proportional to the variance r2. If we use the model calibration factor to interpolate, then the variance r2 includes variance due to error, i.e. the residuals between calculation and measurement. There are 107 assemblies without detector, and Fig. 4 shows the standard deviation (SD) of the interpolation result. The value of r is usually under 3%, so most of the standard deviations of different interpolation results can be less than 1%. If there is less errors between calculation and measurement, the interpolation results will have less uncertainties. To perform a more direct verification of our method, we use a subset of the instrumented assembly power to predict the other instrumented assembly power. Here, the ‘‘leave-one-out’’ validation strategy is used. We suppose that the instruments are failed sequentially; the powers are then calculated for that failed instrument, i.e., treated as non-instrumented; and last, restored as an instrumented assembly. In this case study, only one instrument is failed at a time. The reconstruction errors of different failed instrument are represented in Fig. 5. It shows that the ordinary kriging method can reconstruct the failed instrument power accurately. The assemblies with large reconstruction error all are located in the periphery of the core. High calculation speed is a basic requirement for on-line power distribution monitoring systems. The weights wi of the ordinary kriging interpolation can be calculated off-line, so the power mapping procedure is fast and meet the requirement of on-line monitoring.
Fig. 5. Reconstruction errors of failed instrumented assembly.
3.2. Neutron detector location optimization In-core detector optimal locations searching method based on ordinary kriging and simulated annealing was tested. Objective function, i.e. the average ordinary kriging covariance corresponding to the detector locations was formulated by the ordinary kriging, and its minimum was searched by simulated annealing algorithm, and then the optimal detector locations were obtained as shown in Fig. 6, the red marked assemblies only belong to the current detector location strategy, the black marked assemblies only belong to the optimal detector location strategy, and the blue marked assemblies belong to both location strategies. The number of transitions for a given temperature was chosen as 5, the number of the outer iterations was chosen as 3000, and the initial temperature was chosen as 100. A random detector locations were generated as shown in Fig. 7, the black marked assemblies only belong to the random detector location strategy. The comparison of the AOKV of three location patterns were shown in Table 1, where r is the standard deviation in Eq. (15).
0.35
0.3
SD/
0.25
0.2
0.15
0.1 0
20
40
60
80
100
120
Uninstrumented assembly number Fig. 4. Standard deviation of interpolation result.
Fig. 6. Comparison of the optimal detector location with the current detector location.
X. Peng et al. / Annals of Nuclear Energy 68 (2014) 118–123
123
(1) The reconstruction results show that the ordinary kriging method provides credible results and performs better than the CECOR method. The RMS reconstruction errors are kept at less than 0.35%, and the maximum RE reconstruction errors are kept at less than 1.02% for the entire operating cycle. (2) The ordinary kriging method for power distribution monitoring is fast enough for on-line monitoring. (3) The ordinary kriging method can be used together with combinatorial optimization algorithm, such as simulated annealing, to determine the optimal neutron detector locations. More testing needs to be done to ensure that the ordinary kriging performs good enough in terms of power distribution reconstruction. And also, there is a potential for expanding the ordinary kriging method to a three-dimensional approach that could produce even better results. References
Fig. 7. Comparison of the random detector location with the current detector location.
Table 1 Comparison of the AOKVs. Locations pattern
AOKV (r2)
Current locations Random locations Optimal locations
0.0319 0.0510 0.0309
From Table 1, we can see that the current detector locations of Unit 1 reactor of DayaBay Nuclear Power Plant has a near optimal AOKV, and the random detector locations cannot give a good AOKV. 4. Conclusion and future work This paper presents the concept of using the ordinary kriging method to reconstruct reactor power distributions and determine the optimal neutron detector locations. The DayaBay Unit 1 PWR is taken as an example to verify this method. The following conclusions are drawn from the study:
Bertrand, B. et al., 2012. Robustness of nuclear core activity reconstruction by data assimilation. J. Power Eng. Syst. 6 (2), 289–301. Boyd, A., Miller, W., 1996. The BEACON on-line core monitoring system: functional upgrades and applications. In: In-Core Instrumentation and Reactor Core Assessment. Mito-shi, Japan. Brus, D. et al., 2013. Optimization of sample patterns for universal kriging of environmental variables. Geoderma 138 (1–2), 86–95. Jang, J. et al., 2004. On-line generation of three-dimensional core power distribution using in-core detector signals to monitor safety limits. J. Korean Nucl. Soc. 36 (6), 528–539. Karlson, C., 1995. Continuing advancements in in-core power distribution measurement methods using SIMULATE-3 and CECOR-3.4. Nucl. Sci. Eng. 121, 57–66. Lee, K., Kim, C., 2003. The least-squares method for three-dimensional core power distribution monitoring in pressurized water reactors. Nucl. Sci. Eng. 143 (3), 268–280. Li, W., et al., 2013. Development of a new flux map processing code for moveable detector system in PWR. In: (M&C 2013). Sun Valley, America. Lockwood, B., Anitescu, M., 2012. Gradient-enhanced universal kriging for uncertainty propagation. Nucl. Sci. Eng. 170, 168–195. Pokhrel, R. et al., 2013. A kriging method of interpolation used to map liquefaction potential over alluvial ground. Eng. Geol. 152 (1), 26–37. Tang, T., et al., 1978. Analytical Design of the CANDU-600 On-line Flux Mapping System. TDAI-152, Atomic Energy of Canada Limited. Tojo, M., et al., 2008. Application of the new core monitoring system ‘GNF-ARGOS’ to the start-up tracking calculation of ABWR. In: International Conference on the Physics of Reactors. Interlaken, Switzerland. Wang, K. et al., 1991. Application of cubic spline synthesis in on-line core axial power distribution monitoring. J. Korean Nucl. Soc. 23 (3), 316–320. Webb, R., Brittingham, J., 2000. Comparison of CECOR algorithm to Lagrange multiplier method to estimate reactor power distributions. Nucl. Technol. 132, 206–213. Zheng, S., 2003. GNPS 1&2 18 Month Fuel Cycles Reload Study Procedure. TFJN/DC/ 2322, Framatome ANP. Zhong, B., et al., 2010. A new method on flux mapping by function expansion. In: International Conference on Nuclear Engineering (ICONE18). Xi’an, China.