- Email: [email protected]
Reactivity estimation during a reactivity-initiated accident using the extended Kalman filter
Annals of Nuclear Energy 85 (2015) 753–762
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Reactivity estimation during a reactivity-initiated accident using the extended Kalman filter R. Busquim e Silva a,⇑, A.L.F. Marques a, J.J. Cruz a, K. Shirvan b, M.S. Kazimi b a b
University of Sao Paulo, Av. Prof. Luciano Gualberto, 380, Sao Paulo 05508-010, Brazil Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02142, USA
a r t i c l e
i n f o
Article history: Received 7 April 2015 Accepted 19 June 2015 Available online 20 July 2015 Keywords: Extended Kalman filter Reactivity-initiated accident Inverse point kinetics Nuclear power plant Thermal hydraulic/neutronic coupling
a b s t r a c t This study implements the extended Kalman filter (EKF) to estimate the nuclear reactor reactivity behavior under a reactivity-initiated accident (RIA). A coupled neutronics/thermal hydraulics code PARCS/RELAP5 (P3D/R5) simulates a control rod assembly ejection (CRE) on a traditional 2272 MWt PWR to generate the reactor power profile. A MATLAB script adds random noise to the simulated reactor power. For comparison, the inverse point kinetics (IPK) deterministic method is also implemented. Three different cases of CRE are simulated and the EKF, IPK and the P3D/R5 reactivity are compared. It was found that the EKF method presents better results compared to the IPK method. Furthermore, under a RIA due to small reactivity insertion and slow power response, the IPK reactivity varies widely from positive to negative, which may add extra difficulty to the task of controlling a supercritical reactor. This feature is also confirmed by a sensitivity analysis for five different noise loads and three distinct noise measurements standard deviations (SD). Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction The core neutron population drives the nuclear reaction rate, and resulting power drives the thermal–hydraulic (TH) behavior in light water reactor (LWR). Therefore, the determination of criticality, i.e. the core multiplication factor, using the measured neutron flux, is among the most significant issues in nuclear reactor studies. However, the current online reactivity meters use the deterministic inverse point kinetics (IPK) method (Bhatt et al., 2012; Shimazu and van Rooijen, 2014), which does not account for the stochastic nature of neutron flux interactions with materials. In addition, the neutron flux determination through signal detectors generates noisy data. State estimation is the method to estimate unknown state variables of dynamic systems based on indirect and/or inaccurate series of measurements. One of the best examples of an optimum estimator for linear systems is the Kalman filter (Gelb et al., 1974; Jazwinski, 1970; Anderson and Moore, 1979): an optimal estimator that minimizes the mean square error of the estimated variable if the noise is of Gaussian distribution (or the best linear estimator if it is not). The EKF is an extension of the Kalman filter (KF) for nonlinear systems, which works on a stream of noisy input ⇑ Corresponding author. E-mail address: [email protected] (R. Busquim e Silva). http://dx.doi.org/10.1016/j.anucene.2015.06.031 0306-4549/Ó 2015 Elsevier Ltd. All rights reserved.
data to generate a statistically optimal estimate of the variable of interest. Previous works, related to the use of KF/EKF in the nuclear field, include the study of Zhe (2010) that presents the use of a robust KF, a suboptimal state estimator that guarantees a limit to the error variance, to estimate state variables such as the neutron density, the precursor density, the average reactor fuel temperature, the coolant temperature inside the reactor and the coolant temperature entering the reactor core. Bhatt et al. (2013) presents the comparison between two estimation methods for predicting the shutdown reactivity in a nuclear power plant (NPP) using offline data. On the other hand, Shimazu and van Rooijen (2014) presents qualitative results of the application of KF and IPK. In addition, Valsalam (2013) suggests the use of KF to implement on-line failure detection in the NPP instrumentation. This work implements an EKF algorithm to estimate the reactivity during a RIA. The EKF algorithm uses simulated data obtained from the coupled code P3D/R5, after the addition of random noise. In the P3D/R5 simulation (Busquim e Silva, 2015), the reactor goes supercritical due to a control rod assembly (CRA) ejection, i.e. a partial or total control rod sudden withdraw from its reactor core position, which may occur because of a control rod drive mechanism (or its housing) failure. A CRA ejection accident is a significant event that may lead to RIA conditions.
754
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
In order to assess the reliability and consistency of the results, the IPK method is also implemented and both algorithms (EKF and IPK) are compared to the P3D/R5 simulation. Although both algorithms converged in the performed tests, the EKF has lower noise content, i.e. lower SD, compared to the IPK - thus, the EKF presented superior performance.
Denote by ^ xðtk =tk Þ the estimated state at time tk just after processing the measurements taken up to time tk (seconds). In the same way, Pðtk =tk Þ denotes the error covariance matrix of the estimated state. The pair ^ xðt k =t k Þ and Pðtk =tk Þ represents the EKF state at time tk . Consider a reference trajectory obtained by integrating:
x_ ðtÞ ¼ FðxðtÞÞ
2. The extended Kalman filter
ð2-7Þ
in the interval [t k ,t kþ1 ] with initial conditions The EKF for a continuous time nonlinear stochastic system, assuming discrete time measurement, is implemented using MATLAB (MATLAB, 2014) by a set of equations that models the estimated plant-sensor system dynamics. The estimated state is calculated according to the EKF block diagram presented in Fig. 2-1, where: FðxðtÞ is the state model; xðtÞ is the state vector; Hðxðtk Þ; tk Þ is the observation model; WðtÞ is the state noise, which is assumed to be an independent Gaussian white process with zero mean and spectral density Q0 (t); yðtk Þ is the measurement vector and Vðt k Þ is the measurement noise, which is also assumed to be independent Gaussian white process with zero mean and covariance Rðt k Þ; and T means the discretization/sampling method according to the time step of the P3D/R5 simulation. The following equations describe the non linear system and the process and measurement noise:
dxðtÞ ¼ FðxðtÞÞ þ WðtÞ dt
ð2-1Þ
yðt k Þ ¼ Hðxðtk ÞÞ þ Vðt k Þ
ð2-2Þ
E½WðtÞ ¼ 0
ð2-3Þ
h i E WðtÞ W T ðsÞ ¼ Q 0 ðtÞ dðt sÞ
ð2-4Þ
E½Vðt k Þ ¼ 0
ð2-5Þ
h i E Vðtk Þ V T ðt m Þ ¼ Rðt k Þ dkm ; where dkm is the Kronecker delta, i.e,
dkm ¼
ð2-6Þ
xðt k Þ ¼ ^xðtk =t k Þ
ð2-8Þ
Defining:
dxðtÞ ¼ xðtÞ xðtÞ
And considering it is sufficiently small, the following linear approximation can be taken as:
_ ¼ dxðtÞ
@F dxðtÞ þ WðtÞ @xxðtÞ
ð2-10Þ
The integration of the previous equation in the interval [tk ,tkþ1 ] winds up:
dxðtkþ1 Þ ¼ Uðtkþ1 ; t k Þ dxðtk Þ þ Wðtkþ1 Þ
ð2-11Þ
where Uðtkþ1 ; tk Þ is the state transition matrix associated to the linearized system (2-10), which can be computed by integrating the equation:
@ @F Uðt; tk Þ ¼ Uðt; tk Þ @t @x xðtÞ
ð2-12Þ
in the interval [t k ,t kþ1 ] with initial conditions
Uðt k ; tk Þ ¼ I
ð2-13Þ
Wðtkþ1 Þ is a white Gaussian stochastic process with zero mean and covariance given by
Qðt kþ1 Þ ¼
Z
tkþ1
Uðt kþ1 ; sÞ Q 0 ðsÞ UT ðt kþ1 ; sÞds
ð2-14Þ
tk
Consider the nominal measurement defined by:
1 if k ¼ m
ðtk Þ ¼ Hðxðtk ÞÞ y
0
and let dyðt k Þ be
otherwise
ð2-9Þ
Fig. 2-1. EKF filter block diagram.
ð2-15Þ
755
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
ðt k Þ dyðtk Þ ¼ yðtk Þ y
ð2-16Þ
b.3) Compute the updated error covariance matrix
Then the linearized measurement equation can be written as
dyðtk Þ ¼ Mðtk Þ dxðt k Þ þ Vðtk Þ
Pðt kþ1 =t kþ1 Þ ¼ ½I Kðt kþ1 Þ Mðtkþ1 Þ Pðtkþ1 =tk Þ
ð2-17Þ
½I Kðt kþ1 Þ Mðt kþ1 ÞT þ Kðtkþ1 Þ Rðt kþ1 Þ
where
@H Mðt k Þ ¼ @x xðtk Þ
K T ðtkþ1 Þ ð2-18Þ
Hence, KF in its original linear form (Gelb et al., 1974; Anderson and Moore, 1979) can be applied to the problem given by Eqs. (2-11) and (2-17) giving thus rise to the EKF, which can be summarized by: i. Take the measurement yðtkþ1 Þ and compute dyðt kþ1 Þ by using Eq. (2-16); ii. Apply the KF to compute the estimate d^ xðt kþ1 =t kþ1 Þ of dxðtkþ1 Þ; and xðt kþ1 Þ to obtain ^ xðtkþ1 =tkþ1 Þ. iii. Add ii to xðtk =tk Þ; Pðtk =tk ÞÞ at time t k In other words, given the filter state (^ and the measurement yðt kþ1 Þ taken at tkþ1 , the EKF computes the filter state at t kþ1 by applying the following steps: a. Prediction Step a.1) Compute the predicted value ^ xðt kþ1 =t k Þ by integrating
x_ ¼ FðxðtÞÞ
ð2-23Þ
ð2-19Þ
in the interval t k to t kþ1 with initial conditions xðt k Þ ¼ ^ xðt k =t k Þ. a.2) Compute the predicted error covariance matrix:
Pðtkþ1 =tk Þ ¼ Uðt kþ1 =t k Þ Pðtk =t k Þ UT ðtkþ1 =tk Þ þ Q ðtkþ1 Þ ð2-20Þ b. Update Step b.1) Compute the EKF gain matrix: h i1 Kðtkþ1 Þ ¼ Pðtkþ1 =tk Þ MT ðtkþ1 Þ Mðtkþ1 Þ Pðtkþ1 =tk Þ MT ðtkþ1 Þ þ Rðtkþ1 Þ
The EKF initial state at time t0 ð^ x0 ; P0 Þ is assumed to be given. Hence
^xðt 0 =t 0 Þ ¼ ^x0
ð2-24Þ
Pðt0 =t 0 Þ ¼ P0
ð2-25Þ
By implementing the algorithm above, which works for a continuous system with discrete measurement, the state is estimated at each time step according to the scheme shown in Fig. 2-2.
3. Simulation strategy We evaluate the EKF compared to the IPK algorithm for reactivity estimation using the power output during RIAs from the P3D/R5 2772 MWt NPP simulations results. The cases under investigation are shown in Table 3-1. The input for both algorithms is the reactor power Pðt k Þ generated by the P3D/R5 coupled code under RIA conditions, after the addition of random noise Vðt k Þ. The benchmark strategy is presented in Fig. 3-1. The EKF reactivity is found by applying the stochastic EKF method, the IPK reactivity is the result of a direct application of the deterministic IPK algorithm, and the P3D/R5 reactivity is the result of a CRE evaluated using the coupled P3D/R5 package code. The reference time for all reactivity calculations is the same. The 2772 MWt NPP is modeled using the R5 and P3D codes. The IPK and the estimation tool are based on the point kinetics reactor
ð2-21Þ
b.2) Process the measurement yðt kþ1 Þ and compute the updated estimated state:
^xðt kþ1 =t kþ1 Þ ¼ ^xðt kþ1 =t k Þ þ Kðt kþ1 Þ ½yðtkþ1 Þ Hð^xðt kþ1 =t k ÞÞ ð2-22Þ
Table 3-1 Cases under investigation. Case
Core conditions
Rod worth $
a b c
CRA ejection under hot full power (HFP)⁄ CRA under hot zero power (HZP) Control Rod Bank (CRB) ejection under HZP
0.2 0.4 1.3
⁄ HFP: core at 100% of its nominal power. HZP: core running at only 0.01% of its nominal power.
Fig. 2-2. Prediction and update cycles block diagram.
Fig. 3-1. Reactivity estimation benchmark strategy.
756
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
The PKRE are simple ordinary differential equations, which relate the neutron density and delayed neutron precursor concentration. These equations are usually solved using a six delayed precursor groups (Aboanber and Hamada, 2002; Quintero-Leyva and Core, 2008; Sathiyasheela, 2009): enough to provide an adequate degree of accuracy. The PKRE are:
6 X dnðtÞ qðtÞ b nðtÞ þ ¼ ki C i ðtÞ dt K i¼1
ð3-1Þ
dC i ðtÞ bi ¼ nðtÞ ki C i ðtÞ 1 6 i 6 6 dt K
ð3-2Þ
where:
Fig. 3-2. TH/NK coupling scheme.
equations (PKRE) and are implemented using the MATLAB simulation environment. Vðtk Þ is also generated using a MATLAB script. 3.1. R5/P3D simulation The P3D/R5 RIA simulation strategy starts first with the construction of P3D and R5 nuclear power plant model. These include the core nodalization in the radial and axial directions, their correlation using a mapping scheme, the secondary reactor system and the generation of the input files, as can be shown in Fig. 3-2. Then, the HZP (core at start up) and HFP (reactor under normal operation) rods to be ejected, and the correspondent reactivity insertion $, are identified. Next, the coupled code is applied to simulate the rod ejection within 0.1 s, which is the worst possible scenario for addition of reactivity in LWRs (Busquim e Silva, 2015). However, it must be noted that the ejection time depends on the extent of the core or control rod drive mechanical failure and on the reactor coolant pressure during the failure. The coupling simulation strategy uses a MATLAB script to generate the mapping file used by the general interface (GI) to correlate neutronic (NK) nodes and TH cells. The 2772 MWt NPP modeling is based on reference (Organization for Economic Co-Operation and Development, 1999). The NPP core has a 17x17 assembly configuration. There are 177 fuel assemblies (FA) and 64 reflector assemblies (total of 241 assemblies: each one is 21.81 cm wide). There are 60 CRAs grouped into 7 banks (full length control rods). Radially, the model has 18 thermal–hydraulic (TH) core channels and one channel that represent the reflectors (channel 19). The core active height is 357.12 cm and there is a 21.81 cm high reflector on top and another one at the bottom of the core, adding up to a total height of 400.74 cm. The core operates under typical PWR operating pressure and temperature before the RIAs. The core is also divided axially in 28 neutronic (NK) regions. The reactor is assumed to be at end of cycle, 650 EFPD, equilibrium cycle (24.6 GWD/MT average core exposures). In addition, 29 assembly types with different U235 enrichment (from 4 to 5w/o) are presented in the core. The NPPs model has two once-through steam generator (OTSG) loops, each one with one hot leg and two cold legs. There are four reactor coolant pumps. The prompt neutron lifetime (kÞ is 0.18e-04 s. 3.2. Point Kinetics Reactor Equations (PKRE) Simple core models may be useful only for limited situations: the complete reactor core analysis demands complex models that anticipate the consequences of transients and postulated accidents that must be mitigated by design solutions (Griggs et al., 1982).
-
qðtÞ is the reactivity at time t; nðtÞ is the time dependent neutron density; C i ðtÞ the effective concentration of delayed neutron of group i; b is the total effective delayed neutron fraction. Thus
b¼
6 X bi ; i¼1
- bi is the effective delayed neutron fraction of group i; - ki is the effective decay constant of group i; - K is the mean neutron generation time, which is defined as
K¼
‘ ¼ ‘ ð1 qÞ; kðtÞ
- ‘ is the mean lifetime of a neutron in the a reactor; and - kðtÞ is the multiplication factor. The delayed neutron precursor initial concentration is given by:
C i ðt ¼ 0Þ ¼
bi Pð0Þ 1 6 i 6 6 ki K
ð3-3Þ
3.3. IPK reactor equation The IPK reactor equation correlates the reactivity qðtÞ as a function of the time-dependent power level PðtÞ. Therefore, it may be used to calculate the reactivity from a specific power profile. The IPK is derived directly from Eqs. (2-1) and (2-2) after simple manipulation (Duderstadt and Hamilton, 1976):
qðtÞ ¼ b þ K
d ½ln PðtÞ b dt
Z 0
1
DðsÞ
Pðt sÞ ds PðtÞ
ð3-4Þ
where DðsÞ is defined as the ‘‘Delayed Neutron Kernel’’:
DðsÞ ¼
6 X ki bi i¼1
b
eki s
ð3-5Þ
3.4. State-space model of the PKRE for the EKF The implementation of the EKF estimation tool depends on the adequate mathematical modeling of the phenomenon, according to the following state-space form:
dxðtÞ ¼ FðxðtÞÞ þ WðtÞ ¼ An xðtÞ þ Bn þ WðtÞ dt
ð3-6Þ
The representation of the reactivity as a state variable, similar to the work on references (Duderstadt and Hamilton, 1976 and Bhatt et al., 2013), uses the following modeling equations:
qðtÞ ¼ a þ w t
ð3-7Þ
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
dqðtÞ ¼w dt
ð3-8Þ
dwðtÞ ¼0 dt
ð3-9Þ
S¼
q0 n0 ‘
ð3-10Þ
where a and w stand for the stochastic reactivity jump and its slope, respectively, and S is the neutron source modeled as a function of the initial reactivity and neutron density. The state is defined as:
x ¼ ½ n C1
C2
C3
C4
C5
C6
T
q w ;
ð3-11Þ
where the superscript T denotes the transpose. After simple manipulation of the previous equations, we build the matrices An and Bn . The system is nonlinear due to the term qn þ S. The only terms of the Jacobian A ¼ @F that are different from ‘ @x the An matrix are the Að1; 1Þ ¼ qb and the Að1; 8Þ ¼ n‘ . ‘
2
b‘ k1
6 b1 6 6‘ 6 b2 6‘ 6 6 b3 6‘ 6 b An ¼ 6 6 ‘4 6 6 b5 6‘ 6b 6 6 6‘ 6 40 0 Bn ¼
qn ‘
k2
k3
k4
k5
k6
0
k1
0
0
0
0
0
0
0
k2
0
0
0
0
0
0
0
k3
0
0
0
0
0
0
0
k4
0
0
0
0
0
0
0
k5
0
0
0
0
0
0
0
0
k6
0
0
0
0
0
0
0
0
0
0
0
0
0
0
þS 0 0 0 0 0 0 0 0
T
0
3
7 07 7 7 07 7 7 07 7 07 7 7 07 7 7 07 7 7 15 0
ð3-12Þ
ð3-13Þ
4. Results and discussion 4.1. EKF tuning The EKF tuning process includes the task of selecting the filter initial conditions and the adequate values for the matrices Q0 and R (Cruz, 1981). The matrix Q0 is symmetric, positive semi-definite, and it models the noise existing in the following state variables: neutron flux, delayed neutron precursor concentrations, reactivity and its slope. The neutron flux is the main source of noise fluctuations. The evaluation of the error covariance matrix relies basically on the evaluation of the Uðtkþ1 =tk Þ Pðt k =t k Þ UT ðtkþ1 =tk Þ term of Eq. (2-20) (Shimazu and van Rooijen, 2014). This evaluation winds up with a 1=‘ dependence of elements (1, 1) and (8, 8), where 1=‘ is a very large number. In addition, element Q0 (9, 9) is direct related to the reactivity slope. Therefore, only the elements directly related (see Eq. (3-11)) to the neutron flux/reactivity, which are Q0 (1, 1), Q0 (8, 8) and Q0 (9, 9), are assumed to be nonzero. This assumption was also proved to be correct by simulations. The initial Q 0 ði; iÞ-entries are chosen based on a representative change in the variables (Gelb et al., 1974) during the time step Dt, i.e. by assuming that the state variable is likely to change by an amount Dy over the interval of interest Dt :
Q 0ii ¼
ðDyi Þ2 Dt
output channel (this work takes a single sensor). The measurement error, due to the flux sensor meter uncertainty, is assumed to be 3% (Shimazu and van Rooijen, 2014). The initial value of the slope is assumed to be w = 0 for the three cases under analysis. The time step is obtained automatically from the P3D/R5 simulation. The matrix Q is updated at each time step and it represents the independent stochastic fluctuations in the state variables. In addition, a sensitivity analysis for five different 3% error noise loads, and for 1%, 2% and 3% uncertainties, are performed. 4.2. Reactivity assessment using the IPK and EKF The power responses (and the reactivity behavior) during the RIA are very different among the three cases under investigation: case a has a fast increase in the reactivity, a power peak of about 25% due to the lack of reactivity control, and a slow power return to steady-state (Fig. 4-1); case b has a very slow power dynamics (Fig. 4-5); and case c has a sharp power peak (Fig. 4-10). These different trends are explained by the reactivity worth and the core initial conditions for each case. Nevertheless, the random noise fluctuates stochastically around the reactor power, as shown in Fig. 4-2 (both time and power scale zooms for better noise visualization). Figs. 4-3 and 4-4 show that for the smallest (0.2$) but fast reactivity insertion, which is followed by a moderate feedback effect, both IPK and EKF algorithms converge. However, the SD is higher for the IPK simulation results due to its deterministic method of solving the point kinetics equations, as it is seen in Table 4-1. Also, the steady state EKF reactivity is closer to the P3D/R5 results compared to the IPK reactivity. On the other hand, a HZP insertion of 0.4$ (case bÞ has a very slow power response together with a power increase from its initial value of more than 480 times within 100 s, as can be seen in Fig. 4-5. For this case, the EKF presents superior results compared to the IPK because: i) The IPK reactivity varies from positive to negative (from 120$ to 78$) within 12 s of simulation, as presented in Figs. 4-7 and 4-8 (reactivity zoom scale), due to the noise in the measurement and the deterministic nature of IPK algorithm – this magnitude depends on the noise load, as can be seen in Fig. 4-15. This may add extra difficulty to the task of controlling the reactor based on the deterministic method. While EKF overprotects the peak reactivity by 20% and matches the P3D/R5 reactivity prediction during most of the simulation time as shown in Fig. 4-6; and
ð4-1Þ
Due to the fast flux response to the insertion of reactivity and the distinct power trend among the three cases under evaluation shown in Table 3-1, these Q0 initial entries are calculated using their variation up to the first time step. The covariance matrix R represents the measurement error due to the sensor meters per
757
Fig. 4-1. (Case a) P3D/R5 power distribution with and without noise.
758
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762 Table 4-1 EKF and IPK standard deviation $ for 3% power noise. Case
a b c
EKF SD
IPK SD
Ratio IPK/EKF SD
SD$
Mean $
SD$
Mean $
0.001 0.074 0.091
0.014 0.005 0.001
0.002 0.449 0.181
0.011 0.001 0.002
1.50 6.10 1.99
Fig. 4-2. (Case a) P3D/R5 power distribution with and without noise zoom.
Fig. 4-5. (Case b) P3D/R5 power distribution with and without noise.
Fig. 4-3. (Case a) P3D/R5 and EKF reactivity.
Fig. 4-6. (Case b) P3D/R5 and EKF reactivity.
Fig. 4-4. (Case a) P3D/R5 and IPK reactivity.
ii) The reactivity’s SD is higher for the IPK method. It is worth noting that the use of IPK is more effective during normal operation, when noise detectors are less sensitive to
reactor parameters and more immune to noise content in the measured signals (Bhatt et al., 2012). The EKF, as can be seen in Fig. 4-6, underpredicts the reactivity during the initial few seconds of simulation. This is due to the w (reactivity slope) initial guess wðt¼0Þ ¼ 0. The simulations for the three cases under investigation (see Table 3-1) assume the same set of initial conditions to compare the results independent of the power and reactivity dynamics. Nevertheless, the initial w guess dictates how long it takes for the EKF reactivity to follow the P3D/R5 simulation results. The w ¼ 0 initial condition proved to be the best guess independent of the case being simulated. Nevertheless, Fig. 4-9 presents simulation results for wðt¼0Þ ¼ 0:01, which leads to a EKF SD of 0.04$ (see Table 4-1 for
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
759
Fig. 4-7. (Case b) P3D/R5 and IPK reactivity.
Fig. 4-10. (Case c) P3D/R5 Power distribution with and without noise.
Fig. 4-8. (Case b) P3D/R5 and IPK reactivity zoom.
Fig. 4-11. (Case c) P3D/R5 and EKF reactivity.
Fig. 4-9. (Case b) P3D/R5 reactivity for w = 0.01.
Fig. 4-12. (Case c) P3D/R5 and IPK reactivity.
SD comparison). It can be seen that the initial delay between the P3D/R5 and the EKF reactivity was significantly reduced. The application of EKF and IPK under a HZP with reactivity insertion of 1.3$ (case c) is shown in Figs. 4-10, 4-11 and 4-12
(due to the power scale, the noise is not seen in Fig. 4-10), where very sharp peaks in reactivity and power are simulated. In this case, the simulation results of both IPK and EKF present similar behavior, although the IPK reactivity SD is higher, as show in
760
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
Fig. 4-13. (Case b) Sensitivity analysis for EKF reactivity.
Fig. 4-16. (Case b) Sensitivity analysis for IPK reactivity zoom.
Table 4-2 EKF and IPK standard deviation $ reactivity for five 3% noise loads. Case
Noise Load #
EKF SD $
IPK SD $
Ratio IPK/EKF SD
a
1 2 3 4 5
0.012 0.012 0.013 0.013 0.013
0.018 0.016 0.019 0.017 0.018
1.50 1.33 1.46 1.31 1.38
b
1 2 3 4 5
0.074 0.065 0.068 0.083 0.079
0.449 0.404 0.419 0.511 0.512
6.10 6.20 6.16 6.18 6.45
c
1 2 3 4 5
0.091 0.021 0.021 0.023 0.018
0.181 0.044 0.044 0.044 0.045
1.99 2.09 2.09 1.91 2.49
Fig. 4-14. (Case b) Sensitivity analysis for EKF reactivity zoom.
P3D/R5 reactivity, and the IPK SD is calculated over the difference between the IPK and the P3D/R5 reactivity. 4.3. Sensitivity analysis for different noise loads
Fig. 4-15. (Case b) Sensitivity analysis for IPK reactivity.
Table 4-1 (the effect of a sudden increase in power, which could probably damage the fuel, has not been covered within this work). To benchmark results against the coupled code, in this study the EKF SD is calculated over the difference between the EKF and the
Results from the section 4.2 suggest that, under the core supercritical conditions, the use of EKF for reactivity estimation provides better results for the three cases under analysis, especially for a very slow power response to a small insertion of reactivity (case b). Therefore, a sensitivity analysis of five distinct carry-over effects of different random noise loads, applied to that case, are presented in Figs. 4-13 and 4-14 (EKF), and Figs. 4-15 and 4-16 (IPK). It can be seen that independent of the noise load applied to the power profile, the case b relative method performance (between the IPK and EKF) is similar to the one discussed in the previous section: the EKF presents superior results. To be able to assess the impact that changes in the noise load will have on both the EKF and IPK methods for all cases under investigation, the five noise loads have the same 3% SD but different initial random seeds, which are automatically generated by MATLAB. As it can be seen in Table 4-2, the simulation results indicate that the EKF method is less sensitive to the differences of noise loads for all cases. Furthermore, it can be seen that for a sharp increase of reactor power (case cÞ, the performance of the IPK is about 2 times worse than the EKF (see the column ‘‘Ratio IPK/EKF SD’’ in Table 4-2), while for a HZP slow reactor power
761
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
Fig. 4-17. (Case b) EKF measurements noise sensitivity analysis.
Fig. 4-20. (Case b) IPK measurement noise sensitivity analysis zoom.
Table 4-3 EKF and IPK reactivity different noise SD $ sensitivity. Case
Noise load (%)
EKF SD $
IPK SD $
Ratio IPK/EKF SD
a
1 2 3
0.012 0.013 0.012
0.014 0.017 0.018
1.17 1.31 1.50
b
1 2 3
0.063 0.068 0.074
0.316 0.388 0.449
5.02 5.70 6.10
c
1 2 3
0.019 0.021 0.091
0.015 0.024 0.181
0.77 1.18 1.99
4.4. Sensitivity analysis for different measurements noise SD
Fig. 4-18. (Case b) EKF measurement noise sensitivity analysis zoom.
In order to evaluate how effective the EKF is to the noise content in the measurement, we simulate the reactivity response to 1%, 2% and 3% noise SD. Figs. 4-17 and 4-18 present the estimated reactivity response against the P3D/R5 results for case b. As expected, the EKF method predicts well the reactivity associated to the P3D/R5 simulation for the 3 distinct noise loads. On the other hand, the IPK simulation results for the same measurements, as shown in Figs. 4-19 and 4-20, indicate that as the noise content increases, the error between the IPK and P3D/R5 reactivity also increases, as expected. Table 4-3 presents the reactivity standard deviation for all cases under investigation, assuming a 1%, 2% and 3% random noise power profile. It can be seen that EKF provides the best results compared to IPK and the IPK standard deviation increases, due to its deterministic calculation method, as the measurement noise contents is higher. The ratio IPK/EKF SD shows the improvements of the EKF compared to the IPK as the noise content increases. Also, it confirms the results from previous sections: the EKF exhibits better performance for a slow insertion of reactivity under HZP (case bÞ. 4.5. Computation running time
Fig. 4-19. (Case b) IPK measurement noise sensitivity analysis zoom.
profile (case bÞ, the IPK performance is about 6 times worse. The IPK method is more effective in the HFP case because the on-line neutron flux meters are less immune to noise content in the input signals under the HFP conditions.
By the application of similar coding techniques and by running the simulations in the same computer, the computation time is presented in Table 4-4. The running time depends on the dynamics of the case being simulated and on the length of the simulation (50 s for case a, 250 s for case b and only 10 s for case cÞ. It can be seen, under the stated assumptions, that the computation time is not a constraint or a key decision point for choosing EKF or IPK
762
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
Table 4-4 Running time (seconds). Case
EKF elapsed time (seconds)
IPK elapsed time (seconds)
Ratio EKF/IPK elapsed time
a b c
3.41 12.66 1.02
3.19 11.54 0.95
1.07 1.10 1.08
techniques in reactivity estimation. Both techniques are very fast on today’s computers. The simulations run on iMac OS X 2.7 GHz i7 Intel 16Gb. It is worth noted that the calculation of the Delayed Neutron Kernel is the time demanding IPK code module. The EKF and IPK algorithms implementation depends on the adequate use of MATLAB nested functions, such as ‘‘random’’, ‘‘load’’ and ‘‘ones’’; on the application of well know basic mathematic techniques, such as ‘‘Taylor series’’ and ‘‘Euler method’’; and on the use of systematic programming techniques. Thus, the same coding approach was used in both IPK and EKF MATLAB scripts. Furthermore, the nowadays instrumentation and control systems computer capacity will not be a key issue for both IPK and EKF implementation.
controlling a reactor that may go supercritical. This was also confirmed by a sensitivity analysis for five different noises loads and for three distinct SD noise measurements. It is worth noting that although the P3D/R5 simulates the reactivity using a NK spatial kinetics method, the use of the PKRE to model the EKF to estimate the reactivity for a given power profile under noise environment gave accurate results. The benefits of using a stochastic tool such as EKF should be weighed against the use of the well known IPK method. The EKF requires special attention to the filter’s tuning task while the application of the IPK method is quite straightforward. In addition, by applying the same coding approach for both algorithms, the computation time is not a constraint or a key decision point for choosing between IPK and EKF for reactivity estimation. As a conclusion, this work indicates that, with the appropriate choice of the Rand Q 0 matrices, based on the sensor error (RÞ and on the initial variation of some of the state variables ðQ 0 Þ, the EKF is a powerful method for accurate estimation of reactivity under various RIA conditions. Acknowledgments The authors gratefully acknowledge the financial support of the Sao Paulo Research Foundation (FAPESP).
5. Conclusions References In this work, we have implemented the EKF method to estimate the reactivity behavior under noisy reactor power measurements. Furthermore, the EKF reactivity is compared to the IPK and to the P3D/R5 reactivity. The P3D/R5 coupled code generates the reactor power input, without noise, for both IPK and EKF algorithms. A MATLAB script adds noise to the P3D/R5 reactor power. In addition, this implementation includes the following features: (a) The reactor goes supercritical upon initiation of RIA. (b) The system is described by continuous time nonlinear stochastic differential equations. (c) The power measurement occurs at discrete times. (d) The EKF initial conditions are the same independent of the case under analysis. (e) The tuning of the filter goes automatically up to the simulation first time step. (f) The covariance matrix of the state noise is updated at each time step. The performances of the EKF and IPK algorithms are assessed by the convergence to the P3D/R5 simulation and by the reactivity’s SD. In the tests performed, it was found that the IPK reactivity had higher noise content compared to the EKF reactivity for all cases under investigation (Table 3-1). Therefore, the EKF method presented more accurate results when compared to the IPK method. Furthermore, under a RIA where small reactivity insertion and slower power response were simulated (case bÞ, the IPK reactivity varied widely from positive to negative and at times order of magnitude different, which may add extra difficulty to the task of
Aboanber, A.E., Hamada, Y.M., 2002. PWS: an efficient code system for solving space independent nuclear reactor dynamics. Ann. Nucl. Energy 29, 2159–2172. Anderson, B.D., Moore, J.B., 1979. Optimal Filtering. Prentice-Hall, Englewood Cliffs. Bhatt, T.U. et al., 2013. Estimation of sub-criticality using extended Kalman filtering technique. Ann. Nucl. Energy 60, 98–105. Bhatt, T.U., Shmmith, SR., Tiwari, A.P., 2012. Estimation of Shutdown Reactivity in Power Reactor: Comparative Study of Inverse Point Kinetics and Kalman Filtering, BARC Newsletter, Issue 324. Busquim e Silva, R., 2015. Application of Advanced Computational Tools for Reactivity and Power Analysis of Nuclear Power Plants under ReactivityInitiated Accidents, University of Sao Paulo, PhD Thesis. Cruz, J.J., 1981. Pilotagem Automática de Embarcações com Emprego de Controle Estocástico (in Portuguese). University of Sao Paulo, Master Thesis. Duderstadt, J., Hamilton, L., 1976. Nuclear Reactor Analysis. John Wiley and Sons. Gelb, A. et al., 1974. Applied Optimal Estimation. The MIT Press. Griggs, D.P., Kazimi, M.S., Henry, A.F., 1982. Advanced Methods Development for LWR Transient Analysis. Final Report, MIT. Jazwinski, A.H., 1970. Stochastic Processes and Filtering Theory. Academic Press. MATLAB, The Language of Technical Computing,, accessed on October, 25th 2014. Organization for Economic Co-Operation and Development, 1999. Pressurized Water Reactor Main Steam Line Break (MSLB) Benchmark: Volume I: Final Specifications. Quintero-Leyva, Core, B., 2008. A numerical algorithm to solve the point kinetics equations. Ann. Nucl. Energy 35, 2136–2138. Sathiyasheela, T., 2009. Power series solution method for solving point kinetics equations with lumped model temperature and feedback. Ann. Nucl. Energy 36, 246–250. Shimazu, Y., van Rooijen, W.R.G., 2014. Qualitative performance comparison of reactivity estimation between the extended kalman filter technique and the inverse point kinetic method. Ann. Nucl. Energy 66, 161–166. Valsalam, S.R., 2013. Applications of Kalman Filtering in Industrial Control, Centre for Development of Advanced Computing, , accessed in September 22th 2013. Zhe, D., 2010. Robust Kalman Filter with Application to State Estimation of a Nuclear Reactor, Institute of Nuclear and New Energy, Intech.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Reactivity estimation during a reactivity-initiated accident using the extended Kalman filter R. Busquim e Silva a,⇑, A.L.F. Marques a, J.J. Cruz a, K. Shirvan b, M.S. Kazimi b a b
University of Sao Paulo, Av. Prof. Luciano Gualberto, 380, Sao Paulo 05508-010, Brazil Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02142, USA
a r t i c l e
i n f o
Article history: Received 7 April 2015 Accepted 19 June 2015 Available online 20 July 2015 Keywords: Extended Kalman filter Reactivity-initiated accident Inverse point kinetics Nuclear power plant Thermal hydraulic/neutronic coupling
a b s t r a c t This study implements the extended Kalman filter (EKF) to estimate the nuclear reactor reactivity behavior under a reactivity-initiated accident (RIA). A coupled neutronics/thermal hydraulics code PARCS/RELAP5 (P3D/R5) simulates a control rod assembly ejection (CRE) on a traditional 2272 MWt PWR to generate the reactor power profile. A MATLAB script adds random noise to the simulated reactor power. For comparison, the inverse point kinetics (IPK) deterministic method is also implemented. Three different cases of CRE are simulated and the EKF, IPK and the P3D/R5 reactivity are compared. It was found that the EKF method presents better results compared to the IPK method. Furthermore, under a RIA due to small reactivity insertion and slow power response, the IPK reactivity varies widely from positive to negative, which may add extra difficulty to the task of controlling a supercritical reactor. This feature is also confirmed by a sensitivity analysis for five different noise loads and three distinct noise measurements standard deviations (SD). Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction The core neutron population drives the nuclear reaction rate, and resulting power drives the thermal–hydraulic (TH) behavior in light water reactor (LWR). Therefore, the determination of criticality, i.e. the core multiplication factor, using the measured neutron flux, is among the most significant issues in nuclear reactor studies. However, the current online reactivity meters use the deterministic inverse point kinetics (IPK) method (Bhatt et al., 2012; Shimazu and van Rooijen, 2014), which does not account for the stochastic nature of neutron flux interactions with materials. In addition, the neutron flux determination through signal detectors generates noisy data. State estimation is the method to estimate unknown state variables of dynamic systems based on indirect and/or inaccurate series of measurements. One of the best examples of an optimum estimator for linear systems is the Kalman filter (Gelb et al., 1974; Jazwinski, 1970; Anderson and Moore, 1979): an optimal estimator that minimizes the mean square error of the estimated variable if the noise is of Gaussian distribution (or the best linear estimator if it is not). The EKF is an extension of the Kalman filter (KF) for nonlinear systems, which works on a stream of noisy input ⇑ Corresponding author. E-mail address: [email protected] (R. Busquim e Silva). http://dx.doi.org/10.1016/j.anucene.2015.06.031 0306-4549/Ó 2015 Elsevier Ltd. All rights reserved.
data to generate a statistically optimal estimate of the variable of interest. Previous works, related to the use of KF/EKF in the nuclear field, include the study of Zhe (2010) that presents the use of a robust KF, a suboptimal state estimator that guarantees a limit to the error variance, to estimate state variables such as the neutron density, the precursor density, the average reactor fuel temperature, the coolant temperature inside the reactor and the coolant temperature entering the reactor core. Bhatt et al. (2013) presents the comparison between two estimation methods for predicting the shutdown reactivity in a nuclear power plant (NPP) using offline data. On the other hand, Shimazu and van Rooijen (2014) presents qualitative results of the application of KF and IPK. In addition, Valsalam (2013) suggests the use of KF to implement on-line failure detection in the NPP instrumentation. This work implements an EKF algorithm to estimate the reactivity during a RIA. The EKF algorithm uses simulated data obtained from the coupled code P3D/R5, after the addition of random noise. In the P3D/R5 simulation (Busquim e Silva, 2015), the reactor goes supercritical due to a control rod assembly (CRA) ejection, i.e. a partial or total control rod sudden withdraw from its reactor core position, which may occur because of a control rod drive mechanism (or its housing) failure. A CRA ejection accident is a significant event that may lead to RIA conditions.
754
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
In order to assess the reliability and consistency of the results, the IPK method is also implemented and both algorithms (EKF and IPK) are compared to the P3D/R5 simulation. Although both algorithms converged in the performed tests, the EKF has lower noise content, i.e. lower SD, compared to the IPK - thus, the EKF presented superior performance.
Denote by ^ xðtk =tk Þ the estimated state at time tk just after processing the measurements taken up to time tk (seconds). In the same way, Pðtk =tk Þ denotes the error covariance matrix of the estimated state. The pair ^ xðt k =t k Þ and Pðtk =tk Þ represents the EKF state at time tk . Consider a reference trajectory obtained by integrating:
x_ ðtÞ ¼ FðxðtÞÞ
2. The extended Kalman filter
ð2-7Þ
in the interval [t k ,t kþ1 ] with initial conditions The EKF for a continuous time nonlinear stochastic system, assuming discrete time measurement, is implemented using MATLAB (MATLAB, 2014) by a set of equations that models the estimated plant-sensor system dynamics. The estimated state is calculated according to the EKF block diagram presented in Fig. 2-1, where: FðxðtÞ is the state model; xðtÞ is the state vector; Hðxðtk Þ; tk Þ is the observation model; WðtÞ is the state noise, which is assumed to be an independent Gaussian white process with zero mean and spectral density Q0 (t); yðtk Þ is the measurement vector and Vðt k Þ is the measurement noise, which is also assumed to be independent Gaussian white process with zero mean and covariance Rðt k Þ; and T means the discretization/sampling method according to the time step of the P3D/R5 simulation. The following equations describe the non linear system and the process and measurement noise:
dxðtÞ ¼ FðxðtÞÞ þ WðtÞ dt
ð2-1Þ
yðt k Þ ¼ Hðxðtk ÞÞ þ Vðt k Þ
ð2-2Þ
E½WðtÞ ¼ 0
ð2-3Þ
h i E WðtÞ W T ðsÞ ¼ Q 0 ðtÞ dðt sÞ
ð2-4Þ
E½Vðt k Þ ¼ 0
ð2-5Þ
h i E Vðtk Þ V T ðt m Þ ¼ Rðt k Þ dkm ; where dkm is the Kronecker delta, i.e,
dkm ¼
ð2-6Þ
xðt k Þ ¼ ^xðtk =t k Þ
ð2-8Þ
Defining:
dxðtÞ ¼ xðtÞ xðtÞ
And considering it is sufficiently small, the following linear approximation can be taken as:
_ ¼ dxðtÞ
@F dxðtÞ þ WðtÞ @xxðtÞ
ð2-10Þ
The integration of the previous equation in the interval [tk ,tkþ1 ] winds up:
dxðtkþ1 Þ ¼ Uðtkþ1 ; t k Þ dxðtk Þ þ Wðtkþ1 Þ
ð2-11Þ
where Uðtkþ1 ; tk Þ is the state transition matrix associated to the linearized system (2-10), which can be computed by integrating the equation:
@ @F Uðt; tk Þ ¼ Uðt; tk Þ @t @x xðtÞ
ð2-12Þ
in the interval [t k ,t kþ1 ] with initial conditions
Uðt k ; tk Þ ¼ I
ð2-13Þ
Wðtkþ1 Þ is a white Gaussian stochastic process with zero mean and covariance given by
Qðt kþ1 Þ ¼
Z
tkþ1
Uðt kþ1 ; sÞ Q 0 ðsÞ UT ðt kþ1 ; sÞds
ð2-14Þ
tk
Consider the nominal measurement defined by:
1 if k ¼ m
ðtk Þ ¼ Hðxðtk ÞÞ y
0
and let dyðt k Þ be
otherwise
ð2-9Þ
Fig. 2-1. EKF filter block diagram.
ð2-15Þ
755
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
ðt k Þ dyðtk Þ ¼ yðtk Þ y
ð2-16Þ
b.3) Compute the updated error covariance matrix
Then the linearized measurement equation can be written as
dyðtk Þ ¼ Mðtk Þ dxðt k Þ þ Vðtk Þ
Pðt kþ1 =t kþ1 Þ ¼ ½I Kðt kþ1 Þ Mðtkþ1 Þ Pðtkþ1 =tk Þ
ð2-17Þ
½I Kðt kþ1 Þ Mðt kþ1 ÞT þ Kðtkþ1 Þ Rðt kþ1 Þ
where
@H Mðt k Þ ¼ @x xðtk Þ
K T ðtkþ1 Þ ð2-18Þ
Hence, KF in its original linear form (Gelb et al., 1974; Anderson and Moore, 1979) can be applied to the problem given by Eqs. (2-11) and (2-17) giving thus rise to the EKF, which can be summarized by: i. Take the measurement yðtkþ1 Þ and compute dyðt kþ1 Þ by using Eq. (2-16); ii. Apply the KF to compute the estimate d^ xðt kþ1 =t kþ1 Þ of dxðtkþ1 Þ; and xðt kþ1 Þ to obtain ^ xðtkþ1 =tkþ1 Þ. iii. Add ii to xðtk =tk Þ; Pðtk =tk ÞÞ at time t k In other words, given the filter state (^ and the measurement yðt kþ1 Þ taken at tkþ1 , the EKF computes the filter state at t kþ1 by applying the following steps: a. Prediction Step a.1) Compute the predicted value ^ xðt kþ1 =t k Þ by integrating
x_ ¼ FðxðtÞÞ
ð2-23Þ
ð2-19Þ
in the interval t k to t kþ1 with initial conditions xðt k Þ ¼ ^ xðt k =t k Þ. a.2) Compute the predicted error covariance matrix:
Pðtkþ1 =tk Þ ¼ Uðt kþ1 =t k Þ Pðtk =t k Þ UT ðtkþ1 =tk Þ þ Q ðtkþ1 Þ ð2-20Þ b. Update Step b.1) Compute the EKF gain matrix: h i1 Kðtkþ1 Þ ¼ Pðtkþ1 =tk Þ MT ðtkþ1 Þ Mðtkþ1 Þ Pðtkþ1 =tk Þ MT ðtkþ1 Þ þ Rðtkþ1 Þ
The EKF initial state at time t0 ð^ x0 ; P0 Þ is assumed to be given. Hence
^xðt 0 =t 0 Þ ¼ ^x0
ð2-24Þ
Pðt0 =t 0 Þ ¼ P0
ð2-25Þ
By implementing the algorithm above, which works for a continuous system with discrete measurement, the state is estimated at each time step according to the scheme shown in Fig. 2-2.
3. Simulation strategy We evaluate the EKF compared to the IPK algorithm for reactivity estimation using the power output during RIAs from the P3D/R5 2772 MWt NPP simulations results. The cases under investigation are shown in Table 3-1. The input for both algorithms is the reactor power Pðt k Þ generated by the P3D/R5 coupled code under RIA conditions, after the addition of random noise Vðt k Þ. The benchmark strategy is presented in Fig. 3-1. The EKF reactivity is found by applying the stochastic EKF method, the IPK reactivity is the result of a direct application of the deterministic IPK algorithm, and the P3D/R5 reactivity is the result of a CRE evaluated using the coupled P3D/R5 package code. The reference time for all reactivity calculations is the same. The 2772 MWt NPP is modeled using the R5 and P3D codes. The IPK and the estimation tool are based on the point kinetics reactor
ð2-21Þ
b.2) Process the measurement yðt kþ1 Þ and compute the updated estimated state:
^xðt kþ1 =t kþ1 Þ ¼ ^xðt kþ1 =t k Þ þ Kðt kþ1 Þ ½yðtkþ1 Þ Hð^xðt kþ1 =t k ÞÞ ð2-22Þ
Table 3-1 Cases under investigation. Case
Core conditions
Rod worth $
a b c
CRA ejection under hot full power (HFP)⁄ CRA under hot zero power (HZP) Control Rod Bank (CRB) ejection under HZP
0.2 0.4 1.3
⁄ HFP: core at 100% of its nominal power. HZP: core running at only 0.01% of its nominal power.
Fig. 2-2. Prediction and update cycles block diagram.
Fig. 3-1. Reactivity estimation benchmark strategy.
756
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
The PKRE are simple ordinary differential equations, which relate the neutron density and delayed neutron precursor concentration. These equations are usually solved using a six delayed precursor groups (Aboanber and Hamada, 2002; Quintero-Leyva and Core, 2008; Sathiyasheela, 2009): enough to provide an adequate degree of accuracy. The PKRE are:
6 X dnðtÞ qðtÞ b nðtÞ þ ¼ ki C i ðtÞ dt K i¼1
ð3-1Þ
dC i ðtÞ bi ¼ nðtÞ ki C i ðtÞ 1 6 i 6 6 dt K
ð3-2Þ
where:
Fig. 3-2. TH/NK coupling scheme.
equations (PKRE) and are implemented using the MATLAB simulation environment. Vðtk Þ is also generated using a MATLAB script. 3.1. R5/P3D simulation The P3D/R5 RIA simulation strategy starts first with the construction of P3D and R5 nuclear power plant model. These include the core nodalization in the radial and axial directions, their correlation using a mapping scheme, the secondary reactor system and the generation of the input files, as can be shown in Fig. 3-2. Then, the HZP (core at start up) and HFP (reactor under normal operation) rods to be ejected, and the correspondent reactivity insertion $, are identified. Next, the coupled code is applied to simulate the rod ejection within 0.1 s, which is the worst possible scenario for addition of reactivity in LWRs (Busquim e Silva, 2015). However, it must be noted that the ejection time depends on the extent of the core or control rod drive mechanical failure and on the reactor coolant pressure during the failure. The coupling simulation strategy uses a MATLAB script to generate the mapping file used by the general interface (GI) to correlate neutronic (NK) nodes and TH cells. The 2772 MWt NPP modeling is based on reference (Organization for Economic Co-Operation and Development, 1999). The NPP core has a 17x17 assembly configuration. There are 177 fuel assemblies (FA) and 64 reflector assemblies (total of 241 assemblies: each one is 21.81 cm wide). There are 60 CRAs grouped into 7 banks (full length control rods). Radially, the model has 18 thermal–hydraulic (TH) core channels and one channel that represent the reflectors (channel 19). The core active height is 357.12 cm and there is a 21.81 cm high reflector on top and another one at the bottom of the core, adding up to a total height of 400.74 cm. The core operates under typical PWR operating pressure and temperature before the RIAs. The core is also divided axially in 28 neutronic (NK) regions. The reactor is assumed to be at end of cycle, 650 EFPD, equilibrium cycle (24.6 GWD/MT average core exposures). In addition, 29 assembly types with different U235 enrichment (from 4 to 5w/o) are presented in the core. The NPPs model has two once-through steam generator (OTSG) loops, each one with one hot leg and two cold legs. There are four reactor coolant pumps. The prompt neutron lifetime (kÞ is 0.18e-04 s. 3.2. Point Kinetics Reactor Equations (PKRE) Simple core models may be useful only for limited situations: the complete reactor core analysis demands complex models that anticipate the consequences of transients and postulated accidents that must be mitigated by design solutions (Griggs et al., 1982).
-
qðtÞ is the reactivity at time t; nðtÞ is the time dependent neutron density; C i ðtÞ the effective concentration of delayed neutron of group i; b is the total effective delayed neutron fraction. Thus
b¼
6 X bi ; i¼1
- bi is the effective delayed neutron fraction of group i; - ki is the effective decay constant of group i; - K is the mean neutron generation time, which is defined as
K¼
‘ ¼ ‘ ð1 qÞ; kðtÞ
- ‘ is the mean lifetime of a neutron in the a reactor; and - kðtÞ is the multiplication factor. The delayed neutron precursor initial concentration is given by:
C i ðt ¼ 0Þ ¼
bi Pð0Þ 1 6 i 6 6 ki K
ð3-3Þ
3.3. IPK reactor equation The IPK reactor equation correlates the reactivity qðtÞ as a function of the time-dependent power level PðtÞ. Therefore, it may be used to calculate the reactivity from a specific power profile. The IPK is derived directly from Eqs. (2-1) and (2-2) after simple manipulation (Duderstadt and Hamilton, 1976):
qðtÞ ¼ b þ K
d ½ln PðtÞ b dt
Z 0
1
DðsÞ
Pðt sÞ ds PðtÞ
ð3-4Þ
where DðsÞ is defined as the ‘‘Delayed Neutron Kernel’’:
DðsÞ ¼
6 X ki bi i¼1
b
eki s
ð3-5Þ
3.4. State-space model of the PKRE for the EKF The implementation of the EKF estimation tool depends on the adequate mathematical modeling of the phenomenon, according to the following state-space form:
dxðtÞ ¼ FðxðtÞÞ þ WðtÞ ¼ An xðtÞ þ Bn þ WðtÞ dt
ð3-6Þ
The representation of the reactivity as a state variable, similar to the work on references (Duderstadt and Hamilton, 1976 and Bhatt et al., 2013), uses the following modeling equations:
qðtÞ ¼ a þ w t
ð3-7Þ
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
dqðtÞ ¼w dt
ð3-8Þ
dwðtÞ ¼0 dt
ð3-9Þ
S¼
q0 n0 ‘
ð3-10Þ
where a and w stand for the stochastic reactivity jump and its slope, respectively, and S is the neutron source modeled as a function of the initial reactivity and neutron density. The state is defined as:
x ¼ ½ n C1
C2
C3
C4
C5
C6
T
q w ;
ð3-11Þ
where the superscript T denotes the transpose. After simple manipulation of the previous equations, we build the matrices An and Bn . The system is nonlinear due to the term qn þ S. The only terms of the Jacobian A ¼ @F that are different from ‘ @x the An matrix are the Að1; 1Þ ¼ qb and the Að1; 8Þ ¼ n‘ . ‘
2
b‘ k1
6 b1 6 6‘ 6 b2 6‘ 6 6 b3 6‘ 6 b An ¼ 6 6 ‘4 6 6 b5 6‘ 6b 6 6 6‘ 6 40 0 Bn ¼
qn ‘
k2
k3
k4
k5
k6
0
k1
0
0
0
0
0
0
0
k2
0
0
0
0
0
0
0
k3
0
0
0
0
0
0
0
k4
0
0
0
0
0
0
0
k5
0
0
0
0
0
0
0
0
k6
0
0
0
0
0
0
0
0
0
0
0
0
0
0
þS 0 0 0 0 0 0 0 0
T
0
3
7 07 7 7 07 7 7 07 7 07 7 7 07 7 7 07 7 7 15 0
ð3-12Þ
ð3-13Þ
4. Results and discussion 4.1. EKF tuning The EKF tuning process includes the task of selecting the filter initial conditions and the adequate values for the matrices Q0 and R (Cruz, 1981). The matrix Q0 is symmetric, positive semi-definite, and it models the noise existing in the following state variables: neutron flux, delayed neutron precursor concentrations, reactivity and its slope. The neutron flux is the main source of noise fluctuations. The evaluation of the error covariance matrix relies basically on the evaluation of the Uðtkþ1 =tk Þ Pðt k =t k Þ UT ðtkþ1 =tk Þ term of Eq. (2-20) (Shimazu and van Rooijen, 2014). This evaluation winds up with a 1=‘ dependence of elements (1, 1) and (8, 8), where 1=‘ is a very large number. In addition, element Q0 (9, 9) is direct related to the reactivity slope. Therefore, only the elements directly related (see Eq. (3-11)) to the neutron flux/reactivity, which are Q0 (1, 1), Q0 (8, 8) and Q0 (9, 9), are assumed to be nonzero. This assumption was also proved to be correct by simulations. The initial Q 0 ði; iÞ-entries are chosen based on a representative change in the variables (Gelb et al., 1974) during the time step Dt, i.e. by assuming that the state variable is likely to change by an amount Dy over the interval of interest Dt :
Q 0ii ¼
ðDyi Þ2 Dt
output channel (this work takes a single sensor). The measurement error, due to the flux sensor meter uncertainty, is assumed to be 3% (Shimazu and van Rooijen, 2014). The initial value of the slope is assumed to be w = 0 for the three cases under analysis. The time step is obtained automatically from the P3D/R5 simulation. The matrix Q is updated at each time step and it represents the independent stochastic fluctuations in the state variables. In addition, a sensitivity analysis for five different 3% error noise loads, and for 1%, 2% and 3% uncertainties, are performed. 4.2. Reactivity assessment using the IPK and EKF The power responses (and the reactivity behavior) during the RIA are very different among the three cases under investigation: case a has a fast increase in the reactivity, a power peak of about 25% due to the lack of reactivity control, and a slow power return to steady-state (Fig. 4-1); case b has a very slow power dynamics (Fig. 4-5); and case c has a sharp power peak (Fig. 4-10). These different trends are explained by the reactivity worth and the core initial conditions for each case. Nevertheless, the random noise fluctuates stochastically around the reactor power, as shown in Fig. 4-2 (both time and power scale zooms for better noise visualization). Figs. 4-3 and 4-4 show that for the smallest (0.2$) but fast reactivity insertion, which is followed by a moderate feedback effect, both IPK and EKF algorithms converge. However, the SD is higher for the IPK simulation results due to its deterministic method of solving the point kinetics equations, as it is seen in Table 4-1. Also, the steady state EKF reactivity is closer to the P3D/R5 results compared to the IPK reactivity. On the other hand, a HZP insertion of 0.4$ (case bÞ has a very slow power response together with a power increase from its initial value of more than 480 times within 100 s, as can be seen in Fig. 4-5. For this case, the EKF presents superior results compared to the IPK because: i) The IPK reactivity varies from positive to negative (from 120$ to 78$) within 12 s of simulation, as presented in Figs. 4-7 and 4-8 (reactivity zoom scale), due to the noise in the measurement and the deterministic nature of IPK algorithm – this magnitude depends on the noise load, as can be seen in Fig. 4-15. This may add extra difficulty to the task of controlling the reactor based on the deterministic method. While EKF overprotects the peak reactivity by 20% and matches the P3D/R5 reactivity prediction during most of the simulation time as shown in Fig. 4-6; and
ð4-1Þ
Due to the fast flux response to the insertion of reactivity and the distinct power trend among the three cases under evaluation shown in Table 3-1, these Q0 initial entries are calculated using their variation up to the first time step. The covariance matrix R represents the measurement error due to the sensor meters per
757
Fig. 4-1. (Case a) P3D/R5 power distribution with and without noise.
758
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762 Table 4-1 EKF and IPK standard deviation $ for 3% power noise. Case
a b c
EKF SD
IPK SD
Ratio IPK/EKF SD
SD$
Mean $
SD$
Mean $
0.001 0.074 0.091
0.014 0.005 0.001
0.002 0.449 0.181
0.011 0.001 0.002
1.50 6.10 1.99
Fig. 4-2. (Case a) P3D/R5 power distribution with and without noise zoom.
Fig. 4-5. (Case b) P3D/R5 power distribution with and without noise.
Fig. 4-3. (Case a) P3D/R5 and EKF reactivity.
Fig. 4-6. (Case b) P3D/R5 and EKF reactivity.
Fig. 4-4. (Case a) P3D/R5 and IPK reactivity.
ii) The reactivity’s SD is higher for the IPK method. It is worth noting that the use of IPK is more effective during normal operation, when noise detectors are less sensitive to
reactor parameters and more immune to noise content in the measured signals (Bhatt et al., 2012). The EKF, as can be seen in Fig. 4-6, underpredicts the reactivity during the initial few seconds of simulation. This is due to the w (reactivity slope) initial guess wðt¼0Þ ¼ 0. The simulations for the three cases under investigation (see Table 3-1) assume the same set of initial conditions to compare the results independent of the power and reactivity dynamics. Nevertheless, the initial w guess dictates how long it takes for the EKF reactivity to follow the P3D/R5 simulation results. The w ¼ 0 initial condition proved to be the best guess independent of the case being simulated. Nevertheless, Fig. 4-9 presents simulation results for wðt¼0Þ ¼ 0:01, which leads to a EKF SD of 0.04$ (see Table 4-1 for
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
759
Fig. 4-7. (Case b) P3D/R5 and IPK reactivity.
Fig. 4-10. (Case c) P3D/R5 Power distribution with and without noise.
Fig. 4-8. (Case b) P3D/R5 and IPK reactivity zoom.
Fig. 4-11. (Case c) P3D/R5 and EKF reactivity.
Fig. 4-9. (Case b) P3D/R5 reactivity for w = 0.01.
Fig. 4-12. (Case c) P3D/R5 and IPK reactivity.
SD comparison). It can be seen that the initial delay between the P3D/R5 and the EKF reactivity was significantly reduced. The application of EKF and IPK under a HZP with reactivity insertion of 1.3$ (case c) is shown in Figs. 4-10, 4-11 and 4-12
(due to the power scale, the noise is not seen in Fig. 4-10), where very sharp peaks in reactivity and power are simulated. In this case, the simulation results of both IPK and EKF present similar behavior, although the IPK reactivity SD is higher, as show in
760
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
Fig. 4-13. (Case b) Sensitivity analysis for EKF reactivity.
Fig. 4-16. (Case b) Sensitivity analysis for IPK reactivity zoom.
Table 4-2 EKF and IPK standard deviation $ reactivity for five 3% noise loads. Case
Noise Load #
EKF SD $
IPK SD $
Ratio IPK/EKF SD
a
1 2 3 4 5
0.012 0.012 0.013 0.013 0.013
0.018 0.016 0.019 0.017 0.018
1.50 1.33 1.46 1.31 1.38
b
1 2 3 4 5
0.074 0.065 0.068 0.083 0.079
0.449 0.404 0.419 0.511 0.512
6.10 6.20 6.16 6.18 6.45
c
1 2 3 4 5
0.091 0.021 0.021 0.023 0.018
0.181 0.044 0.044 0.044 0.045
1.99 2.09 2.09 1.91 2.49
Fig. 4-14. (Case b) Sensitivity analysis for EKF reactivity zoom.
P3D/R5 reactivity, and the IPK SD is calculated over the difference between the IPK and the P3D/R5 reactivity. 4.3. Sensitivity analysis for different noise loads
Fig. 4-15. (Case b) Sensitivity analysis for IPK reactivity.
Table 4-1 (the effect of a sudden increase in power, which could probably damage the fuel, has not been covered within this work). To benchmark results against the coupled code, in this study the EKF SD is calculated over the difference between the EKF and the
Results from the section 4.2 suggest that, under the core supercritical conditions, the use of EKF for reactivity estimation provides better results for the three cases under analysis, especially for a very slow power response to a small insertion of reactivity (case b). Therefore, a sensitivity analysis of five distinct carry-over effects of different random noise loads, applied to that case, are presented in Figs. 4-13 and 4-14 (EKF), and Figs. 4-15 and 4-16 (IPK). It can be seen that independent of the noise load applied to the power profile, the case b relative method performance (between the IPK and EKF) is similar to the one discussed in the previous section: the EKF presents superior results. To be able to assess the impact that changes in the noise load will have on both the EKF and IPK methods for all cases under investigation, the five noise loads have the same 3% SD but different initial random seeds, which are automatically generated by MATLAB. As it can be seen in Table 4-2, the simulation results indicate that the EKF method is less sensitive to the differences of noise loads for all cases. Furthermore, it can be seen that for a sharp increase of reactor power (case cÞ, the performance of the IPK is about 2 times worse than the EKF (see the column ‘‘Ratio IPK/EKF SD’’ in Table 4-2), while for a HZP slow reactor power
761
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
Fig. 4-17. (Case b) EKF measurements noise sensitivity analysis.
Fig. 4-20. (Case b) IPK measurement noise sensitivity analysis zoom.
Table 4-3 EKF and IPK reactivity different noise SD $ sensitivity. Case
Noise load (%)
EKF SD $
IPK SD $
Ratio IPK/EKF SD
a
1 2 3
0.012 0.013 0.012
0.014 0.017 0.018
1.17 1.31 1.50
b
1 2 3
0.063 0.068 0.074
0.316 0.388 0.449
5.02 5.70 6.10
c
1 2 3
0.019 0.021 0.091
0.015 0.024 0.181
0.77 1.18 1.99
4.4. Sensitivity analysis for different measurements noise SD
Fig. 4-18. (Case b) EKF measurement noise sensitivity analysis zoom.
In order to evaluate how effective the EKF is to the noise content in the measurement, we simulate the reactivity response to 1%, 2% and 3% noise SD. Figs. 4-17 and 4-18 present the estimated reactivity response against the P3D/R5 results for case b. As expected, the EKF method predicts well the reactivity associated to the P3D/R5 simulation for the 3 distinct noise loads. On the other hand, the IPK simulation results for the same measurements, as shown in Figs. 4-19 and 4-20, indicate that as the noise content increases, the error between the IPK and P3D/R5 reactivity also increases, as expected. Table 4-3 presents the reactivity standard deviation for all cases under investigation, assuming a 1%, 2% and 3% random noise power profile. It can be seen that EKF provides the best results compared to IPK and the IPK standard deviation increases, due to its deterministic calculation method, as the measurement noise contents is higher. The ratio IPK/EKF SD shows the improvements of the EKF compared to the IPK as the noise content increases. Also, it confirms the results from previous sections: the EKF exhibits better performance for a slow insertion of reactivity under HZP (case bÞ. 4.5. Computation running time
Fig. 4-19. (Case b) IPK measurement noise sensitivity analysis zoom.
profile (case bÞ, the IPK performance is about 6 times worse. The IPK method is more effective in the HFP case because the on-line neutron flux meters are less immune to noise content in the input signals under the HFP conditions.
By the application of similar coding techniques and by running the simulations in the same computer, the computation time is presented in Table 4-4. The running time depends on the dynamics of the case being simulated and on the length of the simulation (50 s for case a, 250 s for case b and only 10 s for case cÞ. It can be seen, under the stated assumptions, that the computation time is not a constraint or a key decision point for choosing EKF or IPK
762
R. Busquim e Silva et al. / Annals of Nuclear Energy 85 (2015) 753–762
Table 4-4 Running time (seconds). Case
EKF elapsed time (seconds)
IPK elapsed time (seconds)
Ratio EKF/IPK elapsed time
a b c
3.41 12.66 1.02
3.19 11.54 0.95
1.07 1.10 1.08
techniques in reactivity estimation. Both techniques are very fast on today’s computers. The simulations run on iMac OS X 2.7 GHz i7 Intel 16Gb. It is worth noted that the calculation of the Delayed Neutron Kernel is the time demanding IPK code module. The EKF and IPK algorithms implementation depends on the adequate use of MATLAB nested functions, such as ‘‘random’’, ‘‘load’’ and ‘‘ones’’; on the application of well know basic mathematic techniques, such as ‘‘Taylor series’’ and ‘‘Euler method’’; and on the use of systematic programming techniques. Thus, the same coding approach was used in both IPK and EKF MATLAB scripts. Furthermore, the nowadays instrumentation and control systems computer capacity will not be a key issue for both IPK and EKF implementation.
controlling a reactor that may go supercritical. This was also confirmed by a sensitivity analysis for five different noises loads and for three distinct SD noise measurements. It is worth noting that although the P3D/R5 simulates the reactivity using a NK spatial kinetics method, the use of the PKRE to model the EKF to estimate the reactivity for a given power profile under noise environment gave accurate results. The benefits of using a stochastic tool such as EKF should be weighed against the use of the well known IPK method. The EKF requires special attention to the filter’s tuning task while the application of the IPK method is quite straightforward. In addition, by applying the same coding approach for both algorithms, the computation time is not a constraint or a key decision point for choosing between IPK and EKF for reactivity estimation. As a conclusion, this work indicates that, with the appropriate choice of the Rand Q 0 matrices, based on the sensor error (RÞ and on the initial variation of some of the state variables ðQ 0 Þ, the EKF is a powerful method for accurate estimation of reactivity under various RIA conditions. Acknowledgments The authors gratefully acknowledge the financial support of the Sao Paulo Research Foundation (FAPESP).
5. Conclusions References In this work, we have implemented the EKF method to estimate the reactivity behavior under noisy reactor power measurements. Furthermore, the EKF reactivity is compared to the IPK and to the P3D/R5 reactivity. The P3D/R5 coupled code generates the reactor power input, without noise, for both IPK and EKF algorithms. A MATLAB script adds noise to the P3D/R5 reactor power. In addition, this implementation includes the following features: (a) The reactor goes supercritical upon initiation of RIA. (b) The system is described by continuous time nonlinear stochastic differential equations. (c) The power measurement occurs at discrete times. (d) The EKF initial conditions are the same independent of the case under analysis. (e) The tuning of the filter goes automatically up to the simulation first time step. (f) The covariance matrix of the state noise is updated at each time step. The performances of the EKF and IPK algorithms are assessed by the convergence to the P3D/R5 simulation and by the reactivity’s SD. In the tests performed, it was found that the IPK reactivity had higher noise content compared to the EKF reactivity for all cases under investigation (Table 3-1). Therefore, the EKF method presented more accurate results when compared to the IPK method. Furthermore, under a RIA where small reactivity insertion and slower power response were simulated (case bÞ, the IPK reactivity varied widely from positive to negative and at times order of magnitude different, which may add extra difficulty to the task of
Aboanber, A.E., Hamada, Y.M., 2002. PWS: an efficient code system for solving space independent nuclear reactor dynamics. Ann. Nucl. Energy 29, 2159–2172. Anderson, B.D., Moore, J.B., 1979. Optimal Filtering. Prentice-Hall, Englewood Cliffs. Bhatt, T.U. et al., 2013. Estimation of sub-criticality using extended Kalman filtering technique. Ann. Nucl. Energy 60, 98–105. Bhatt, T.U., Shmmith, SR., Tiwari, A.P., 2012. Estimation of Shutdown Reactivity in Power Reactor: Comparative Study of Inverse Point Kinetics and Kalman Filtering, BARC Newsletter, Issue 324. Busquim e Silva, R., 2015. Application of Advanced Computational Tools for Reactivity and Power Analysis of Nuclear Power Plants under ReactivityInitiated Accidents, University of Sao Paulo, PhD Thesis. Cruz, J.J., 1981. Pilotagem Automática de Embarcações com Emprego de Controle Estocástico (in Portuguese). University of Sao Paulo, Master Thesis. Duderstadt, J., Hamilton, L., 1976. Nuclear Reactor Analysis. John Wiley and Sons. Gelb, A. et al., 1974. Applied Optimal Estimation. The MIT Press. Griggs, D.P., Kazimi, M.S., Henry, A.F., 1982. Advanced Methods Development for LWR Transient Analysis. Final Report, MIT. Jazwinski, A.H., 1970. Stochastic Processes and Filtering Theory. Academic Press. MATLAB, The Language of Technical Computing,