- Email: [email protected]
Decay Ratio estimation in boiling water reactors based on the empirical mode decomposition and the Hilbert–Huang transform
Progress in Nuclear Energy 71 (2014) 122e133
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
Decay Ratio estimation in boiling water reactors based on the empirical mode decomposition and the HilberteHuang transform Alfonso Prieto-Guerrero, Gilberto Espinosa-Paredes* División de Ciencias Básicas e Ingeniería, Ingeniería de Procesos e Hidráulica, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, México, D.F. 09340, Mexico
a r t i c l e i n f o
a b s t r a c t
Article history: Received 16 July 2013 Received in revised form 23 November 2013 Accepted 27 November 2013
In this paper a new method based on the empirical mode decomposition (EMD) to estimate a parameter associated with instability in boiling water reactors (BWR), is explored. This instability parameter is not exactly the classical Decay Ratio (DR), but it will be associated with this. The proposed method allows to decompose the analyzed signal in different levels or intrinsic mode functions (IMF). One or more of these different modes can be associated to the instability problem in BWRs. By tracking the instantaneous frequency (obtained through the HilberteHuang transform) and the autocorrelation function of the IMF associated to the instability of the BWR, the estimation of the proposed instability parameter can be achieved. The methodology was validated with two events reported in the Forsmark stability benchmark. Ó 2013 Elsevier Ltd. All rights reserved.
Keywords: Decay Ratio Empirical mode decomposition Boiling water reactors Intrinsic mode functions HilberteHuang transform
1. Introduction In the last few decades, research has been devoted to the study of power oscillations and the mechanisms that generate them (e.g. Saha and Zuber, 1978; Peng et al., 1984; Lahey and Podoswski, 1989; March-Leuba, 1990; March-Leuba and Blakeman, 1991; MarchLeuba and Rey, 1993; Cai et al., 2009). Several approaches have been taken to address the stability of BWRs, March-Leuba (1986) pioneered the study of reduced-order models for coupled thermal-neutron dynamics, followed by Turso et al. (1997) and Muñoz-Cobo et al. (2004), among others. Some of those works were developed to gain insights about the BWRs dynamics, while others were focused on a more complete description of the heat transfer process (e.g. Uehiro et al., 1996; Guido et al., 1997; Podowski and Pinheiro, 1997). The models presented in the mentioned works were able to describe, in a qualitative way, low-frequency oscillations and even instabilities, but neither sustained oscillations of relatively high frequency nor highly non-stationary behavior could be described accurately by such models (Verdu et al., 2001). When the signal is non-stationary Navarro-Esbri et al. (2003) studied the time dependence of the natural frequency using two different tools: the Short-Time Fourier Transform (STFT) and the Time Dependent Power Spectral Density (TDPSD). Both techniques split
* Corresponding author. E-mail address: [email protected] (G. Espinosa-Paredes). 0149-1970/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.pnucene.2013.11.015
the signal in short segments with a high degree of overlapping. Their results show that the STFT provides better accuracy than the TDPSD despite the lower sensitivity to noisy neutronic signals of the TDPSD method. It is worth mentioning that an autoregressive method (AR) was used by those authors to estimate the Decay Ratio (DR). However one of the most important disadvantages in the methods above discussed is the large number of floating operations (multiplications and additions) needs to be implemented. Recently, the wavelet theory has been used to explore new alternatives for transient instability analysis (Espinosa-Paredes et al., 2005, 2007). It has been shown that stability depends on several variables such as control rod patterns, void fraction, burnup, inlet mass flow, among others. A key point is that in general, BWR signals are non-stationary, therefore traditional methods such as the Fourier transform, might lead to biased stability parameters. Sunde and Pázsit (2007) proposed an original work using the wavelet transform in combination with the autocorrelation function (ACF) to estimate the DR. Prieto-Guerrero and Espinosa-Paredes (2008) propose the application of wavelet ridges to track the instantaneous frequency and determine the DR. Also recently (TorresFernández et al., 2010a,b, 2012), the idea of an instantaneous DR was introduced. In this work, a new method based on empirical mode decomposition (EMD) to estimate a parameter associated to instability in BWRs is introduced. This instability parameter is not exactly the classical DR, but it will be associated with this. The methodology is based on the implementation of the empirical mode decomposition
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
algorithm that allows the decomposition of the analyzed signal in different levels or intrinsic mode functions (IMF). One or more of these different modes can be associated to the instability problem in BWRs. Based on the HilberteHuang transform it is possible to get the instantaneous frequency (IF) associated to each IMF. By tracking this instantaneous frequency and the autocorrelation function of the IMF associated to the instability of the BWR, the estimation of the proposed instability parameter can be achieved, this is the main contribution in this work. The effectiveness of EMD has been demonstrated before for processing biomedical signals (Wu and Huang, 2009), audio signals (Khaldi et al., 2009), mechanical signals (Peng et al., 2005), gravitational waves (Lin et al., 2009), geophysical signals (Huang and Wu, 2008), among others. In the nuclear signals domain, there had been a previous research effort on how to estimate the DR of BWRs using EMD (Montesinos et al., 2003). It is worth mentioning the differences between previous research and the present one. The research in of Montesinos et al. (2003) deals with estimating the DR of BWRs by also using EMD and IMF which, conceptually speaking, have no significantly difference from the method employed in the present paper as mentioned before. However, the proposed method is different in two important points: first, we consider the signal (from BWRs) in short segments of 15 s, tracking the estimated instantaneous frequency and second, the IMF associated to detected IF is processed in order to get a parameter similar to the DR along time. Montesinos et al. (2003) considers the complete signal of the IMF in order to estimate the impulse response of the system based on an autoregressive (AR) model. The classical global DR is then obtained based on this impulse response. In counterpart, in our method the proposed Decay Ratio is calculated directly on the autocorrelation function of the IMF associated to the instability phenomenon. To validate our method, simulated and real neutronic signals were used. The methodology was validated with two cases (4 and 6) reported in the Forsmark stability benchmark (Verdu et al., 2001). The rest of this paper is organized as follows: in Section 2 the basic background to understand our methodology is presented. In Section 3, the methodology to estimate the instantaneous frequency and the proposed DR is discussed, additionally a simulated case is proposed to validate our hypothesis. Then in Section 4, the validation of the methodology presented in this paper is performed doing experiments with real signals taken from the Forsmark stability benchmark. Last, in Section 5, our conclusions are presented.
2.
3. 4.
5.
1 and Emin ¼ fxmin ðt1 Þ; xmin ðt2 Þ; xmin ðt3 Þ; .; xmin ðtL Þg. We also set x(t) ¼ s0(t). 1 The set points of Emax are interpolated to form the upper en1 velope of the signal, b x U ðtÞ. Similarly, the set points of Emin are interpolated to form the minimum envelope, b x L ðtÞ. The average envelope b x A ðtÞ ¼ b x U ðtÞ þ b x U ðtÞ=2, is computed. This average envelope is subtracted from the original signal x(t) resulting in a residue signal: r k ðtÞ ¼ xðtÞ b x A ðtÞ with k indicating the iteration, and k ¼ 1 for the first iteration. The iteration on k is continued until the scalar product hrk(t), rk(t)i ¼ 0 and the number of extreme (maxima and minima) and the number of zero-crossings of rk(t) may differ by no more than one. This sifting process produces the first IMF given by IMFj ðtÞ ¼ rjk ðtÞ with j ¼ 1 obtained at the kth iteration. Following this, the function with s1(t) ¼ s0(t) IMF1(t) is created, and the sifting process is repeated (steps 1e4), resulting in the second IMF, i.e. IMF2(t). Considering this procedure, the other P IMFs are generated until the residue rðtÞ ¼ xðtÞ N j ¼ 1 IMFj ðtÞ is accomplished. The functions IMFj(t), j ¼ 1,2,.,N decompose x(t) and are nearly orthogonal to one another.
The schematic diagram of the EMD algorithm is given in Fig. 1. The sifting process essentially extracts scales of the signal. Since each IMF has only one extreme between any two successive zero crossings, the frequency of the signal can be directly inferred by measuring the distribution of the zero crossings of the signal. Further, the IMF has symmetric envelopes and a zero mean value. Due to these characteristics, the IMF is referred to as being monocomponent. Since the residue is computed by successively subtracting the sifted functions from the original signal, the EMD algorithm is data driven and adaptive, i.e. the basis functions are derived from the signal itself in contrast to the traditional methods where the basis functions are fixed. Furthermore, interpolation is an inexact
2. Preliminaries 2.1. Empirical mode decomposition and intrinsic mode functions The empirical mode decomposition (EMD) algorithm was proposed in Huang et al. (1998) in order to analyze non-stationary signals from non-linear processes. EMD extracts intrinsic oscillatory modes defined by the time scales of oscillation. The components that result from the EMD algorithm are called Intrinsic Mode Functions (IMFs). These obtained IMFs result in a composed AMe FM (Amplitude ModulationeFrequency Modulation) signal. 2.1.1. Sifting process The fundamental step of the empirical mode decomposition (EMD) is the next iterative sifting process: 1. Consider a signal x(t) with M maxima and L minima. The sifting process starts with identifying the extrema of the signal, x(t), given 1 by the sets Emax ¼ fxmax ðt1 Þ; xmax ðt2 Þ; xmax ðt3 Þ; .; xmax ðtM Þg,
123
Fig. 1. Flowchart of the EMD method.
124
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
procedure. The performance of the EMD algorithm is sensitive to the interpolation procedure. The sifting process is defined for continuous signals which means, that the performance of the EMD depends on the sampling rate (Huang et al., 1998). The dependence of the EMD algorithm on these factors makes precludes a general, unique framework for EMD. Thus defining a function space for the EMD algorithm, is an ill posed problem, making it difficult to construct an analytical description of the EMD. However, empirically, the EMD has been shown to be effective in extracting relevant components in a variety of applications involving non-stationary signals. 2.2. HilberteHuang transform Consider the analytic representation of a signal given by: jfðtÞ
Zx ðtÞbxðtÞ þ jb x ðtÞ ¼ AðtÞe
ZN N
xðt sÞ
ps
(1)
ds ;
xðtÞ ¼ Re½Zx ðtÞ ¼ AðtÞcosðfðtÞÞ
(2)
The operation denoted in Eq. (2) is merely the convolution of x(t) with 1/t or it can be thought of as a phase shift of the signal by p/2.
(3)
we can observe that this canonical representation of x(t) corresponds to a signal varying simultaneously in amplitude and phase all the time (exactly like each IMF from the EMD). Based on this representation, Ville (1948) proposed to calculate the instantaneous frequency associated to the instantaneous phase 4(t). From Eq. (1), the instantaneous frequency is defined by:
uinst ðtÞ ¼
where A(t) and 4(t) are the envelope and the instantaneous phase of the signal, respectively. In this equation b x ðtÞ is the Hilbert transform of the signal x(t), which is defined by Delprat et al. (1992):
b x ðtÞ ¼
As a result, x(t) and b x ðtÞ are said to be in quadrature. This complex or analytic signal is completely characterized by their amplitude A(t) and their phase 4(t) with values in the interval [0, 2p), forming the so-called canonical pair. Considering that
dfðtÞ dt
(4)
It is clear that the analytic signal Zx(t) associated with x(t) has the same amplitude and frequency range with x(t), it also comprises the phase information of the original signal x(t), which can be used to compute the instantaneous amplitude and to analyze the frequency performance of x(t). Based on this fact, we can construct the analytic signal corresponding to each IMF using the Hilbert transform, defined in Eq. (2). The combination of the EMD applied to x(t) to generate IMFs, and the Hilbert transform of each IMF is called the HilberteHuang Transform (HHT). The different IMFs resulting from the EMD algorithm are orthogonal to each other. The IMFs, thus represent different time-scales of oscillations, which form a set of bases
Original Signal: Forsmark, Case 4, LPRM 22 80 60 40
0
5
10
Instantaneous Frequency (IF)
15
IF Histrograms 4
IMF 1 0
0
5
10
2
IMF 2
0 -2
0
5
10
0
5
10
2
15
IMF 4
0 -2
0
5
10
IF [Hz]
-5
IF [Hz]
IMF 3
IMF 5 -2
0
5
10
Time [s]
15
IF [Hz]
2
0
5
10
15
0
1
4
0.5
2
0
5
10
15
0
1
10
0.5
5
0
0
5
10
15
0
1
20
0.5
10
0
15
0
2
0
15
5 0
5
0
15
IF [Hz]
-1
IF [Hz]
1
0
5
10
15
0
1
40
0.5
20
0
0
5
10
15
Time [s] Fig. 2. The EMD-HHS algorithm applied to a real signal from a BWR.
0
0
1
2
3
4
5
6
7
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
Time [s]
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
functions. This implies that there is no redundancy in the information contained in the different IMFs. Using this property, a distribution can be constructed from the instantaneous frequencies of each of the IMFs. This distribution is called the HilberteHuang Spectrum (HHS) (Huang et al., 1998). Since the instantaneous frequency of the EMD-HHS approach is not based on the Fourier transform, the time-frequency resolution is not limited by uncertainty. It is worth mentioning that the EMD can be viewed as a reconstruction problem where each IMF contributes to the signal reconstruction. It is also clear that the reconstruction error depends on the levels of the EMD: if the EMD is obtained in all the possible levels (exhaustive problem), the reconstructed signal is practically the original one. In our case, this reconstruction situation is irrelevant because we are only interested in a specific IMF: the one associated to the instability. Fig. 2 plots the first 5 IMFs of the EMD decomposition of real signal obtained from a BWR during 15 s. This BWR signal corresponds to the LPRM number 22 from the case 4 of the Forsmark benchmark (Verdu et al., 2001). This case is presented in detail in Section 4 of this work. We observe that the IMFs have a zero mean, symmetric envelopes (AMeFM signals), and the number of extreme is at most one more than the number of zero crossings. The frequency of an IMF is directly proportional to the number of zero crossings in one time duration of IMFs. Also we can observe perfectly the frequencies tracking instantaneous frequency (IF) associated to each IMF. This IF permits to establish which level of the EMD
125
contains the mode associated to the BWR instability. Based on this, a simple method based on the distribution (histogram) of the IF is implemented. Both IMFs and IFs determine the mode or modes associated to the instability phenomena. In the presented signal, this instability problem can be clearly observed at the IMF 3 of Fig. 2. In this figure, the frame (in red) shows the selected IMF and the corresponding IF and its distribution (histogram). Comparing this IF with those found in other levels of the EMD, it is clear that the level 3 (IMF 3) is the right one to be considered. 3. Decay Ratio estimation using EMDeHHT In this work, the EMDeHHT technique to the BWR signals was applied, following the next methodology: 1. The considered signal (APRM or LPRM) obtained from the BWR is segmented in time windows of 15 s of duration. 2. Each segmented signal (APRM or LPRM) is analyzed (decomposed) using the EMD algorithm, described in Section 2.1, obtaining in this way the corresponding IMFs. It is worth mentioning that the APRM or LPRM signals are not being processed before. For instance, to remove the signal trend, it is not necessary to filter the BWR signal before, due that this information is contained in the last scale computed from the EMD (Fig. 3). It is also worth mentioning that the trend does not represent an IMF. Another important feature of the EMD is that
Original Signal 10 0 -10
0
5
10
15
20
25
30
Instantaneous Frequency (IF)
35
IF Histrograms
0.05
40
IF [Hz]
IMF 1 0 -0.05
0
5
10
15
20
0.5 20 0
25
0
5
10
15
0 0.4
0.45
0.5
0.55
0.6
0.35
0.4
0.45
0.5
0.1
0.2
0.3
0.4
-3
IF [Hz]
5
x 10
IMF 2 0 -5
0
5
10
15
20
10 0
25
20
0.5
0
5
10
15
0 0.3
-3
IF [Hz]
5
x 10
IMF 3 0 -5
0
5
10
15
20
20 0
25
40
0.5
0
5
10
15
0
0
-3
IF [Hz]
2
x 10
IMF 4 0 -2
0
5
10
15
20
25
100
0.5
50 0
0
5
10
15
0
0
0.05
0.1
0
0.05
0.1
0.15
0.2
0.25
0.3
0.15
0.2
0.25
0.3
-3
1
x 10
IF [Hz]
Trend 0 -1
0
5
10
15
Time [s]
20
25
100
0.5
50 0
0
5
10
15
Time [s] Fig. 3. EMD-HHS method applied to a synthetic signal.
0
IF [Hz]
126
3. 4.
5. 6.
7.
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
anytime the analyzed signal should contain an additive and high-frequency noise, it will be isolated in the first scale (IMF 1). Thus completely eliminating it from the signal and the other IMFs (including our IMF of interest) permitting a better analysis. The Hilbert transform for each IMF is computed in order to get the instantaneous frequencies contained in each IMF. When tracking these frequencies, it is possible to get the mode or modes associated to instability processes. In this way, only the mode (IMF) or modes associated to the BWRs instability are considered for further processing. The autocorrelation function (AFC) of the considered IMF is calculated. Based on these ACF values, an instability parameter is estimated. This parameter is determined as the Decay Ratio (DR) but unlike it, is calculated considering two consecutive maxima from the ACF magnitude of the analyzed IMF. This consideration (DR estimation) was inferred based on the analysis of a synthetic signal (see below) and the discussion on the Forsmark stability benchmark (see Section 4). Finally, the mean and variance of the DR are calculated averaging the estimated DR obtained in each window of 15 s along time.
Before taking the challenge of real signals, first we use a synthetic signal to validate our method. For this we considered three different oscillating modes driving a second order system which impulse response is given by:
Table 1 Estimated DR and frequency of the Forsmark stability benchmark, Case 4. Power sensor
IMF
Decay Ratio DR (mean)
DR (standard deviation)
Frequency (Hz)
C4_APRM_1 C4_LPRM_1 C4_LPRM_2 C4_LPRM_3 C4_LPRM_4 C4_LPRM_5 C4_LPRM_6 C4_LPRM_7 C4_LPRM_8 C4_LPRM_9 C4_LPRM_10 C4_LPRM_11 C4_LPRM_12 C4_LPRM_13 C4_LPRM_14 C4_LPRM_15 C4_LPRM_16 C4_LPRM_17 C4_LPRM_18 C4_LPRM_19 C4_LPRM_20 C4_LPRM_21 C4_LPRM_22
2 2 2 2 2 3 3 3 2 3 2 3 3 3 3 3 2 3 2 3 3 3 3
0.7143 0.7843 0.8162 0.8413 0.7272 0.8004 0.8068 0.7593 0.6660 0.7635 0.5217 0.7652 0.7001 0.7519 0.7410 0.6985 0.7837 0.7828 0.7046 0.6610 0.6654 0.7880 0.7014
0.1761 0.1244 0.1132 0.1057 0.1842 0.1220 0.0907 0.1559 0.1436 0.1289 0.1847 0.1415 0.1347 0.1694 0.1802 0.1877 0.1145 0.1192 0.1503 0.1219 0.1689 0.1596 0.1174
0.4718 0.5031 0.4923 0.5024 0.4905 0.5041 0.5299 0.5091 0.4951 0.5174 0.5217 0.4793 0.4873 0.5031 0.5150 0.5178 0.5103 0.4819 0.4824 0.4864 0.4873 0.5270 0.5017
Fig. 4. Decay Ratio estimation of the IMF 1 from the simulation considering a synthetic signal.
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 hðtÞ ¼ A0 ezun t cos 1 z un t
(5)
where un is the natural resonant angular frequency, A0 is the magnitude of this impulse response, and z is the damping parameter. Here the oscillation term is represented by the cosine function, and the damping effect is modeled using an exponentially decaying factor. The Decay Ratio (DR) is calculated by:
0
1
2pz C B DR ¼ exp@ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiA 2 1z
(6)
The input signal applied to this system is established as:
xðtÞ ¼ 5 cosð2pf0 tÞ þ cosð2p1:2f0 tÞ þ cosð2p0:35f0 tÞ
(7)
with f0 ¼ 0.5 Hz. It is clearly that the signal in Eq. (7) will be in resonance with the system, which the impulse response is defined in Eq. (5), at the frequency u0 ¼ 2pf0. The idea behind the creation of this synthetic signal is to prove that the proposed method is well adapted to estimate the DR in a second order system (like in a BWR) considering that the output is composed by different modes
127
attenuated or reinforced by the system along time. Indeed, a real BWR signal from the LPRM or APRM, the density wave phenomenon is presented like a mode around 0.5 Hz growing in amplitude along time. Harmonics of this mode can be added reinforcing its presence. For this simulation a sampling frequency of 12.5 Hz was considered. The input signal given in Eq. (7) was introduced in the system defined by Eq. (5). To obtain the corresponding IMFs and IFs plotted in Fig. 3 the proposed EMDeHHT method was applied to the output signal of this system. Two important remarks can be pointed out: We can observe that the interesting mode can be tracked along time. In this case, corresponding to the first IMF. The IF associated to this IMF is clearly showed in the middle (at the top) of Fig. 3. Values of this IF are confirmed by the histogram showed in the third column of Fig. 3. It is worth mentioning that the ACF obtained from the analyzed IMF vanishes quickly in time after 15 s. Based on the fact that the proposed figure of merit is estimated from this autocorrelation function, this means that there is not relevant information in the LPRM signals (directly on the IMF) after this time window. This is a very important point to consider in order to apply this
Fig. 5. IMFs and IFs of a specific LPRM for Forsmark case 4.
128
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
method on BWRs real signals. It is also important to mention that this result is also comparable to our simulation case. Following with the proposed methodology, we computed the ACF of the first IMF of the output signal. We obtain the ACF magnitude and we detect the maxima contained in this function. Fig. 4 presented three different cases corresponding to DR values of 0.5, 0.8 and 1. In Eq. (5) the respective damping ratios (i.e. z ¼ 0.1097, 0.0361, 0.0152) are considered for simulation proposals. Maxima of the ACF magnitude are calculated and represented by the red marks (*) in Fig. 4. Again two important remarks can be pointed out: Considering the classical DR (ratio of third and first peaks), in the case of DR equal to 0.5 and 1, values are slightly (less than 10%) underestimated. In the case of the DR equal to 0.8, the estimated DR is practically the real value considered in the simulation. We decided to use as figure of merit also the DR but measured over the two first peaks of the ACF magnitude of the considered IMF. This decision was based on the analysis of several real signals from BWRs, where in certain cases the third peak sometimes is not well calculated (sometimes even lost when the IMF is quite small in amplitude) from the EMDeHHT algorithm. In the case of this simulation, we can observe that for DR equal to 0.5 and 0.8 is slightly (less than 10%) overestimated, for a DR equal to 1 is practically the real considered value. This decision
about the DR and its impact in real signals from BRWs will be discussed in the Section 4 of this work.
4. Experiments and results The methodology described in the previous section, was applied to BWRs real signals. Among those, we include two complete cases from the Forsmark stability benchmark corresponding to cases 4 and 6. It is worth mentioning that all cases from the Forsmark stability benchmark were studied. However, for reasons of space only the most challenging cases are presented. The potential of this proposed method to predict an instability event is to implement a real-time stability monitor in a BWR. In our implementation, a modified version of the EMD code provided in Rilling and Goncalves (2008) was used. 4.1. Cases of the Forsmark benchmark The EMDeHHT method was applied to neutronic signals of Forsmark stability benchmark (Verdu et al., 2001). The benchmark is based on data from several measurements performed in the Swedish BWR reactor Forsmark 1 and 2, in the period 1989 to 1997. We discuss the results obtained specially from two cases: 4 and 6. Case 4. This case contains a mixture between a global oscillation mode and a regional (half core) oscillation. This event corresponds
Fig. 6. IF and DR estimation for Forsmark case 4, APRM 1 signal.
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
Forsmark Signal. Case 4. APRM 1 150
100
50
0
0
0.1
0.2
0.3
0.4
0.5
0.6
Instantaneous Frequency (IF)
0.7
0.8
[Hz]
Fig. 7. Histogram of the IF for Forsmark case 4, APRM 1 signal.
to a situation where the neutronic power reactor suffers abnormal and apparently unstable oscillations. The C4_APRM and C4_LPRM_x signals correspond to average power range monitor (APRM) and local power range monitor (LPRM) registers respectively, during the instability event. Results
129
of the proposed method are shown in Table 1. These results agree with those obtained by different methodologies developed by participants and included in the final report of the Forsmark stability benchmark (Verdu et al., 2001). In Fig. 5, a specific LPRM (7) is considered for the EMD-HSS algorithm. In this figure, the IMFs, IFs associated to these IMFs, and the histograms from the instantaneous frequencies are plotted. The analyzed segment of 15 s is illustrated in the top of the figure and marked by the box (in green) in the complete record from the LPRM. It is clear that the first IMF contains only the information that we will considers like noise which contents is distributed along all frequencies as shown by the histogram of the corresponding IF and, applying this algorithm we can consider an automatic denoising of the analyzed signal. IMFs 4 and 5 showed that the frequencies contained in these modes are the tendency towards low part of the spectrum and they are not considered in the analysis. However IMFs 2 and 3 results show the interesting modes to consider for further processing. Histograms of the corresponding IMF 2 and 3 show both modes containing information about stability (frequency around 0.5 Hz), nevertheless the IMF 3 clearly shows the frequency around 0.5 Hz more persistently. Because that, this is the chosen mode to further processing to get the DR parameter. In Table 1 it is also included the number of IMF which is processing in order to get the values of the DR and IF. So having decided which mode (IMF) to be processed, we obtain the ACF of this mode and we get the proposed DR considering only the two consecutives peaks from the magnitude of the ACF. Fig. 6 shows the application of this algorithm to a specific signal, in this
Fig. 8. IF and DR estimation for Forsmark case 4, LPRM 1 signal.
130
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
Decay Ratio. Forsmark Case 4. 1Level 1 Decay Ratio. Forsmark Case 4. Level 1.2 lprm1-pos23 lprm17-pos9 lprm11-pos11 lprm21-pos31
1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0
50
100
150
200
250
300
350
Time [s]
Time [s] Fig. 9. DR estimation along time for Forsmark case 4.
case the APRM 1. In the top of this figure the complete original is plotted. In this figure each analyzed segment of the signal is marked by a green box and plotted immediately below. In the middle of Fig. 6 the corresponding instantaneous frequency (IF) for the chosen IMFs are shown. Two last parts (below) of this figure show the obtained ACF magnitude and the corresponding estimated DR from this ACF. Estimated DR is plotted along time and computed its corresponding mean and standard deviation. It is worth mentioning that this method permits two important things along time: tracking the frequency associated to an instability and corresponding Decay Ratio. In Fig. 6 we can observe that the obtained DR corresponds to those reported in the benchmark by the different applied methodologies (between 0.459 and 0.9, but mostly between 0.7 and 0.85). In our method only two segments of 15 s are under a DR of 0.5 (see Fig. 6), the most part of results are in the
interval going from 0.7 to 0.85 (see Fig. 6), in perfect agree with the benchmark. In Fig. 7 the distribution of the IF from the analyzed IMFs is plotted. We can observe a Gaussian distribution of this IF for this APRM. This distribution is centered in the frequency 0.48 Hz completely in agree with those presented in the benchmark (between 0.48 and 0.51). Fig. 8 corresponds to the signal from the LPRM 1. We can observe how this method perfectly tracks the IF from the IMF. The corresponding DR along time varies from 0.7 to 0.9, only three analyzed segments give a DR lower than 0.6. Most part of estimated DR is between 0.8 and 0.9 in completely agree with results reported in the benchmark. It is clear that the observed mode (IMF 2) indicates a core more unstable in this part of the reactor than the corresponding global measure (APRM 1). This Gaussian distribution is
Table 2 Estimated DR and frequency of the Forsmark stability benchmark, Case 6.1.
Table 3 Estimated DR and frequency of the Forsmark stability benchmark, Case 6.2.
Power sensor
IMF
Decay Ratio DR (mean)
DR (standard deviation)
Frequency (Hz)
Power sensor
IMF
Decay Ratio DR (mean)
DR (standard deviation)
Frequency (Hz)
C6_APRM_1 C6_LPRM_1 C6_LPRM_2 C6_LPRM_3 C6_LPRM_4 C6_LPRM_5 C6_LPRM_6 C6_LPRM_7 C6_LPRM_8 C6_LPRM_9 C6_LPRM_10 C6_LPRM_11 C6_LPRM_12 C6_LPRM_13 C6_LPRM_14 C6_LPRM_15 C6_LPRM_16 C6_LPRM_17 C6_LPRM_18
3 3 3 3 3 3 3 3 3 4 4 4 3 4 3 3 3 4 3
0.6722 0.6183 0.6492 0.5981 0.6288 0.5797 0.6906 0.6148 0.6852 0.6264 0.6466 0.7217 0.5687 0.6295 0.6087 0.5508 0.6242 0.6233 0.6323
0.1078 0.1323 0.1285 0.1555 0.1015 0.1522 0.1254 0.1327 0.1608 0.1573 0.1385 0.1069 0.1592 0.1484 0.1178 0.1399 0.1379 0.1981 0.1416
0.4932 0.5662 0.5016 0.5093 0.5036 0.5079 0.5127 0.4718 0.4639 0.4450 0.4430 0.4278 0.4934 0.4297 0.5499 0.4940 0.5207 0.4739 0.4943
C6_APRM_2 C6_LPRM_1 C6_LPRM_2 C6_LPRM_3 C6_LPRM_4 C6_LPRM_5 C6_LPRM_6 C6_LPRM_7 C6_LPRM_8 C6_LPRM_9 C6_LPRM_10 C6_LPRM_11 C6_LPRM_12 C6_LPRM_13 C6_LPRM_14 C6_LPRM_15 C6_LPRM_16 C6_LPRM_17 C6_LPRM_18
3 3 3 3 3 3 2 3 2 3 3 3 3 3 3 3 3 4 3
0.8682 0.5937 0.6790 0.8529 0.6596 0.8882 0.8679 0.8422 0.8609 0.5445 0.6104 0.6709 0.7714 0.7059 0.7754 0.8425 0.7942 0.6842 0.6100
0.0589 0.1230 0.1776 0.0862 0.2000 0.0544 0.1320 0.1099 0.1089 0.1309 0.1158 0.1493 0.1125 0.1316 0.1281 0.1238 0.1187 0.1312 0.1727
0.4986 0.5494 0.4810 0.5226 0.5020 0.5065 0.5024 0.5065 0.5162 0.4757 0.4734 0.5051 0.5240 0.5318 0.5377 0.5223 0.5214 0.4459 0.5798
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
131
centered at 0.495 Hz and it corresponds to reported values (between 0.48 and 0.5 Hz) in the benchmark. All values (DRs and frequencies) corresponding to the reported LPRMs in the Forsmark benchmark are included in Table 1. To end this case, we pointed out the discussion on the benchmark conclusions, especially about these two remarks (Verdu et al., 2001):
performed close to each other, both in time and in operating conditions. For test 1 (6.1): power at 64.4% and core flow of 4 416 kg/s, for test 2 (6.2): power at 63.3% and core flow of 4 298 kg/s. The C6_APRM and C6_LPRM_x signals correspond to an average power range monitor (APRM) and local power range monitor (LPRM) registers respectively, during the instability event considered.
There is a half of the reactor (locations 23 and 9) where the DR is high and the other half (locations 31 and 11) where the DR is lower. The upper part of the reactor seems to be more stable than the lower part. Spectral analysis of the signals indicates that there is a phase shift between the LPRM at radial locations 23 and 11, and locations 23 and 31, but the out-of-phase oscillation is not totally developed.
Case 6.1. Results of the proposed method are shown in Table 2. These results are completely agreed with those reported in the Forsmark stability benchmark (Verdu et al., 2001). However, we focused the discussion on the second test (case 6.2).
In Fig. 9, the estimated DRs along time for the LPRMs 1, 7, 11 and 21 (positions 23, 9, 11 and 31 respectively) at Level 1 are plotted. In our opinion (in agree with remarks above) the out-of-phase oscillation is not completely developed. In this figure, we can observe that only in certain time regions we have this regional oscillation but it is not presented in the complete analyzed signals. Case 6. This test shows local (channel) oscillations. The data contain APRM and LPRM signals from two tests that were
Case 6.2. Results of the proposed method are shown in Table 3. Like the other shown cases, results are in complete agree with those presented in the benchmark. In Fig. 10 a signal segment from LPRM 15 is plotted. Like in Case 4, we can observe that the IMF 3 contains the frequency associated to instability event (around 0.5 Hz). Almost in all APRM and LPRMs for this case, the IF is tracked on the IMF 3. Fig. 11 shows an analyzed segment from the APRM 2. In the middle of the figure, the tracking of the IF associated to the IMF 3 along time is shown. We can clearly observe estimated values of the IF around 0.5 Hz. In the bottom of Fig. 11, estimated DR is plotted along time. For this APRM, we can see almost all DR values varying
Fig. 10. IMFs and IFs of a specific LPRM for Forsmark case 6.2.
132
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
Fig. 11. IF and DR estimation for Forsmark case 6.2, APRM 2 signal.
Decay Decay Ratio. Forsmark Level 1 Ratio. Forsmark CaseCase 6.2. Level6.2. 1 1
0.9
0.8
0.7
0.6
0.5
lprm1-pos23 0.4
lprm17-pos29 lprm7-pos6 lprm15-pos24
0.3
0.2 0
50
100
150
200
250
Time [s] Time [s]
Fig. 12. DR estimation along time for Forsmark case 6.2.
300
350
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
between 0.8 and 0.95. This indicates an almost unstable event. This situation corresponds to those reported in the benchmark. One important conclusion from this case in the benchmark indicates that one half of the reactor is more stable than the other one. To identify this situation we plotted simultaneously in Fig. 12 the estimated DR obtained from the LPRMs 1, 17, 7 and 15 (positions in the reactor 23, 29, 6 and 24 respectively) at Level 1. In this figure we can observe, practically all the time, how a half of the reactor is unstable (LPRM 7 and 15, in solid lines) presenting oscillations varying almost all between 0.8 and 0.95. The other half of the reactor is stable, with estimated DRs (marked with dotted lines) varying almost all between 0.4 and 0.7. This is a very important result that allows visualizing the potential of the proposed method, capable of determining in real time, not only possible unstable events but also possible out-of-phase oscillations in the reactor core.
5. Conclusions In this work, we propose a novel approach to stability analysis of boiling water reactors (BWRs) based on empirical mode decomposition (EMD) and the HilberteHuang transform (HHT). The method consists of obtaining the different intrinsic modes functions (IMF) associated to the EMD. Applying the HHT it is possible to get from each IMF, the associated instantaneous frequency (IF) along time. This instantaneous frequency allows to decide which IMF will be associated to a possible instability event. In this manner, it is possible to track this instability situation. An estimated Decay Ratio (DR) is obtained time along, considering the magnitude of the autocorrelation function (ACF) of followed IMF. Some important conclusions can be delivered. First, the method is relatively simple to implement and it does not represent a high computational complexity. Second, results obtained from the implementation of this method on the Forsmark stability benchmark, are completely in concordance with those reported by the different methodologies developed by the participants of this benchmark. For Cases 4 and 6, included and discussed in this work, we can conclude that the proposed method clearly contributes on the fact to detect events of instability along time inclusive detecting possible cases of out-of-phase oscillations. Conclusions reported in the benchmark (Verdu et al., 2001) were corroborated by the proposed method in this work. This method, with the proposed estimation for the DR, also allows predicting the instability event through time before the oscillation going to a critical state. References Cai, J., Ishiwatari, Y., Ikejiri, S., Oka, Y., 2009. Thermal and stability considerations for a supercritical water-cooled fast reactor with downward-flow channels during power-raising phase of plant startup. Nucl. Eng. Design 239 (4), 665e679. Delprat, N., Escudie, B., Guillemain, P., Kronland-Martinet, R., Tchamitchian, P., Torresani, B., 1992. Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies. IEEE Trans. Info. Theory 38, 644e664. Espinosa-Paredes, G., Prieto-Guerrero, A., Núñez-Carrera, A., Amador-García, R., 2005. Wavelet-based method for instability analysis in boiling water reactors. Nucl. Technol. 151, 250e260.
133
Espinosa-Paredes, G., Nuñez-Carrera, A., Prieto-Guerrero, A., Ceceñas, M., 2007. Wavelet approach for analysis of neutronic power using data of Ringhals stability benchmark. Nucl. Eng. Design 237, 1009e1015. Guido, G., Converti, J., Clause, A., 1997. Density wave oscillations in parallel channels, an analytical approach. Nucl. Eng. Design 125, 121e139. Huang, N.E., Shen, Z., Long, S.R., Wu, M.C., Shih, H.H., Zheng, Q., Yen, N.C., Tung, C.C., Liu, H.H., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Math. Phys. Eng. Sci. 454, 679e699. Huang, N.E., Wu, Z., 2008. A review on HilberteHuang transform: method and its applications to geophysical studies. Rev. Geophys. 46, 1e23. Khaldi, K., Boudraa, A.O., Turki, M., Chonavel, T., Samaali, I., 2009. Audio Encoding Based on the Empirical Mode Decomposition. EUSIPCO, Glasgow. Lahey, R.T., Podoswski, M.Z., 1989. On the analysis of various instabilities in twophase flow. Multiph. Sci. Technol. 4, 183e371. Lin, S.L., Tung, P.C., Huang, N.E., 2009. Data analysis using a combination of independent component analysis and empirical mode decomposition. Phys. Rev. E. Stat. Nonlin. Soft. Matter. Phys. 79 (6 Pt 2), 066705. March-Leuba, J., 1986. A reduced-order model of boiling water reactor linear dynamics. Nucl. Technol. 75, 15e22. March-Leuba, J., 1990. LAPUR Benchmark against In-phase and Out-of-phase Stability Test. NUREG/CR-5605, ORNL/TM-116210. March-Leuba, J., Blakeman, E.D., 1991. A mechanism for out-of-phase power instabilities in boiling water reactors. Nucl. Sci. Eng. 107, 173e179. March-Leuba, J., Rey, J.M., 1993. Coupled thermal-hydraulic neutronic stabilities, in boiling water reactor a review of the state of the art. Nucl. Eng. Design 145, 97e111. Montesinos, M.E., Muñoz-Cobo, J.L., Pérez, C., 2003. HilberteHuang analysis of BWR neutron detector signals: application to DR calculation and to corrupted signal analysis. Ann. Nucl. Energy 30, 715e727. Muñoz-Cobo, J.L., Chiva, S., Sekhri, A., 2004. A reduced order model of BWR dynamics with subcooled boiling and modal kinetics: application to out of phase oscillations. Ann. Nucl. Energy 31, 1135e1162. Navarro-Esbri, J., Ginestar, D., Verdu, G., 2003. Time dependence of linear stability parameters of a BWR. Prog. Nucl. Energy 43, 187e194. Peng, S.J., Podowski, M.Z., Lahey, R.T., Becker, M., 1984. NUFREQ-NP: a computer code for the instability analysis of boiling water nuclear reactors. Nucl. Sci. Eng. 88, 404e411. Peng, Z.K., Tseb, P.W., Chu, F.L., 2005. An improved HilberteHuang transform and its application in vibration signal analysis. J. Sound Vib. 286, 187e205. Podowski, M.Z., Pinheiro, M., 1997. Modeling and numerical simulation of oscillatory two phase flow with application to boiling water nuclear reactor. Nucl. Eng. Design 177, 179e188. Prieto-Guerrero, A., Espinosa-Paredes, G., 2008. Decay Ratio estimation of BWR signals based on wavelet ridges. Nucl. Sci. Eng. 160, 302e317. Rilling, G., Goncalves, P., 2008. EMD Toolbox for Matlab. http://perso.ens-lyon.fr/ patrick.flandrin/emd.html. Saha, P., Zuber, N., 1978. An analytical study of the thermally induced two-phase flow instabilities including the effects of thermal non-equilibrium. Int. J. Heat Mass Transf. 21, 415e426. Sunde, C., Pázsit, I., 2007. Wavelet techniques for the determination of the Decay Ratio in boiling water reactors. Kerntechnik 72, 7e19. Torrres-Fernández, J.E., Espinosa-Paredes, G., Prieto-Guerrero, A., 2012. Decay Ratio estimation in BWRs using the exponential time-scale representation. Ann. Nucl. Energy 49, 143e153. Torres-Fernández, J.E., Prieto-Guerrero, A., Espinosa-Paredes, G., 2010a. Decay Ratio estimation using the instantaneous autocorrelation function. Nucl. Eng. Design 240, 3559e3570. Torres-Fernández, J.E., Prieto-Guerrero, A., Espinosa-Paredes, G., 2010b. Decay Ratio estimation based on time-frequency representations. Ann. Nucl. Energy 37, 93e102. Turso, J.A., March-Leuba, J., Edwards, M., 1997. A modal-based reduced-order model of BWR out-of-phase instabilities. Ann. Nucl. Energy 12, 921e934. Uehiro, M., Rao, Y.F., Fukuda, K., 1996. Linear stability analysis on instabilities of inphase and out-of-phase modes in boiling water reactors. J. Nucl. Sci. Technol. 33, 628e635. Verdu, G., Ginestar, D., Muñoz-Cobo, J.L., Navarro-Esbri, J., Palomo, M.J., Lansaker, P., Conde, J.M., Recio, M., Sartori, E., 2001. Forsmark 1&2 Stability Benchmark. Time Series Analysis Methods for Oscillations during BWR Operation. Final Report. NEA/NSC/DOC, p. 2. Ville, J., 1948. Theorie et application de la notion de signal analytique. Cables Trans. A (1), 61e74. Wu, M.C., Huang, N.E., 2009. Biomedical Data Processing Using HHT: A Review. Chapter 16. In: Advanced Biosignal Processing. Springer-Verlag, pp. 335e352.
Contents lists available at ScienceDirect
Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene
Decay Ratio estimation in boiling water reactors based on the empirical mode decomposition and the HilberteHuang transform Alfonso Prieto-Guerrero, Gilberto Espinosa-Paredes* División de Ciencias Básicas e Ingeniería, Ingeniería de Procesos e Hidráulica, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186, Col. Vicentina, México, D.F. 09340, Mexico
a r t i c l e i n f o
a b s t r a c t
Article history: Received 16 July 2013 Received in revised form 23 November 2013 Accepted 27 November 2013
In this paper a new method based on the empirical mode decomposition (EMD) to estimate a parameter associated with instability in boiling water reactors (BWR), is explored. This instability parameter is not exactly the classical Decay Ratio (DR), but it will be associated with this. The proposed method allows to decompose the analyzed signal in different levels or intrinsic mode functions (IMF). One or more of these different modes can be associated to the instability problem in BWRs. By tracking the instantaneous frequency (obtained through the HilberteHuang transform) and the autocorrelation function of the IMF associated to the instability of the BWR, the estimation of the proposed instability parameter can be achieved. The methodology was validated with two events reported in the Forsmark stability benchmark. Ó 2013 Elsevier Ltd. All rights reserved.
Keywords: Decay Ratio Empirical mode decomposition Boiling water reactors Intrinsic mode functions HilberteHuang transform
1. Introduction In the last few decades, research has been devoted to the study of power oscillations and the mechanisms that generate them (e.g. Saha and Zuber, 1978; Peng et al., 1984; Lahey and Podoswski, 1989; March-Leuba, 1990; March-Leuba and Blakeman, 1991; MarchLeuba and Rey, 1993; Cai et al., 2009). Several approaches have been taken to address the stability of BWRs, March-Leuba (1986) pioneered the study of reduced-order models for coupled thermal-neutron dynamics, followed by Turso et al. (1997) and Muñoz-Cobo et al. (2004), among others. Some of those works were developed to gain insights about the BWRs dynamics, while others were focused on a more complete description of the heat transfer process (e.g. Uehiro et al., 1996; Guido et al., 1997; Podowski and Pinheiro, 1997). The models presented in the mentioned works were able to describe, in a qualitative way, low-frequency oscillations and even instabilities, but neither sustained oscillations of relatively high frequency nor highly non-stationary behavior could be described accurately by such models (Verdu et al., 2001). When the signal is non-stationary Navarro-Esbri et al. (2003) studied the time dependence of the natural frequency using two different tools: the Short-Time Fourier Transform (STFT) and the Time Dependent Power Spectral Density (TDPSD). Both techniques split
* Corresponding author. E-mail address: [email protected] (G. Espinosa-Paredes). 0149-1970/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.pnucene.2013.11.015
the signal in short segments with a high degree of overlapping. Their results show that the STFT provides better accuracy than the TDPSD despite the lower sensitivity to noisy neutronic signals of the TDPSD method. It is worth mentioning that an autoregressive method (AR) was used by those authors to estimate the Decay Ratio (DR). However one of the most important disadvantages in the methods above discussed is the large number of floating operations (multiplications and additions) needs to be implemented. Recently, the wavelet theory has been used to explore new alternatives for transient instability analysis (Espinosa-Paredes et al., 2005, 2007). It has been shown that stability depends on several variables such as control rod patterns, void fraction, burnup, inlet mass flow, among others. A key point is that in general, BWR signals are non-stationary, therefore traditional methods such as the Fourier transform, might lead to biased stability parameters. Sunde and Pázsit (2007) proposed an original work using the wavelet transform in combination with the autocorrelation function (ACF) to estimate the DR. Prieto-Guerrero and Espinosa-Paredes (2008) propose the application of wavelet ridges to track the instantaneous frequency and determine the DR. Also recently (TorresFernández et al., 2010a,b, 2012), the idea of an instantaneous DR was introduced. In this work, a new method based on empirical mode decomposition (EMD) to estimate a parameter associated to instability in BWRs is introduced. This instability parameter is not exactly the classical DR, but it will be associated with this. The methodology is based on the implementation of the empirical mode decomposition
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
algorithm that allows the decomposition of the analyzed signal in different levels or intrinsic mode functions (IMF). One or more of these different modes can be associated to the instability problem in BWRs. Based on the HilberteHuang transform it is possible to get the instantaneous frequency (IF) associated to each IMF. By tracking this instantaneous frequency and the autocorrelation function of the IMF associated to the instability of the BWR, the estimation of the proposed instability parameter can be achieved, this is the main contribution in this work. The effectiveness of EMD has been demonstrated before for processing biomedical signals (Wu and Huang, 2009), audio signals (Khaldi et al., 2009), mechanical signals (Peng et al., 2005), gravitational waves (Lin et al., 2009), geophysical signals (Huang and Wu, 2008), among others. In the nuclear signals domain, there had been a previous research effort on how to estimate the DR of BWRs using EMD (Montesinos et al., 2003). It is worth mentioning the differences between previous research and the present one. The research in of Montesinos et al. (2003) deals with estimating the DR of BWRs by also using EMD and IMF which, conceptually speaking, have no significantly difference from the method employed in the present paper as mentioned before. However, the proposed method is different in two important points: first, we consider the signal (from BWRs) in short segments of 15 s, tracking the estimated instantaneous frequency and second, the IMF associated to detected IF is processed in order to get a parameter similar to the DR along time. Montesinos et al. (2003) considers the complete signal of the IMF in order to estimate the impulse response of the system based on an autoregressive (AR) model. The classical global DR is then obtained based on this impulse response. In counterpart, in our method the proposed Decay Ratio is calculated directly on the autocorrelation function of the IMF associated to the instability phenomenon. To validate our method, simulated and real neutronic signals were used. The methodology was validated with two cases (4 and 6) reported in the Forsmark stability benchmark (Verdu et al., 2001). The rest of this paper is organized as follows: in Section 2 the basic background to understand our methodology is presented. In Section 3, the methodology to estimate the instantaneous frequency and the proposed DR is discussed, additionally a simulated case is proposed to validate our hypothesis. Then in Section 4, the validation of the methodology presented in this paper is performed doing experiments with real signals taken from the Forsmark stability benchmark. Last, in Section 5, our conclusions are presented.
2.
3. 4.
5.
1 and Emin ¼ fxmin ðt1 Þ; xmin ðt2 Þ; xmin ðt3 Þ; .; xmin ðtL Þg. We also set x(t) ¼ s0(t). 1 The set points of Emax are interpolated to form the upper en1 velope of the signal, b x U ðtÞ. Similarly, the set points of Emin are interpolated to form the minimum envelope, b x L ðtÞ. The average envelope b x A ðtÞ ¼ b x U ðtÞ þ b x U ðtÞ=2, is computed. This average envelope is subtracted from the original signal x(t) resulting in a residue signal: r k ðtÞ ¼ xðtÞ b x A ðtÞ with k indicating the iteration, and k ¼ 1 for the first iteration. The iteration on k is continued until the scalar product hrk(t), rk(t)i ¼ 0 and the number of extreme (maxima and minima) and the number of zero-crossings of rk(t) may differ by no more than one. This sifting process produces the first IMF given by IMFj ðtÞ ¼ rjk ðtÞ with j ¼ 1 obtained at the kth iteration. Following this, the function with s1(t) ¼ s0(t) IMF1(t) is created, and the sifting process is repeated (steps 1e4), resulting in the second IMF, i.e. IMF2(t). Considering this procedure, the other P IMFs are generated until the residue rðtÞ ¼ xðtÞ N j ¼ 1 IMFj ðtÞ is accomplished. The functions IMFj(t), j ¼ 1,2,.,N decompose x(t) and are nearly orthogonal to one another.
The schematic diagram of the EMD algorithm is given in Fig. 1. The sifting process essentially extracts scales of the signal. Since each IMF has only one extreme between any two successive zero crossings, the frequency of the signal can be directly inferred by measuring the distribution of the zero crossings of the signal. Further, the IMF has symmetric envelopes and a zero mean value. Due to these characteristics, the IMF is referred to as being monocomponent. Since the residue is computed by successively subtracting the sifted functions from the original signal, the EMD algorithm is data driven and adaptive, i.e. the basis functions are derived from the signal itself in contrast to the traditional methods where the basis functions are fixed. Furthermore, interpolation is an inexact
2. Preliminaries 2.1. Empirical mode decomposition and intrinsic mode functions The empirical mode decomposition (EMD) algorithm was proposed in Huang et al. (1998) in order to analyze non-stationary signals from non-linear processes. EMD extracts intrinsic oscillatory modes defined by the time scales of oscillation. The components that result from the EMD algorithm are called Intrinsic Mode Functions (IMFs). These obtained IMFs result in a composed AMe FM (Amplitude ModulationeFrequency Modulation) signal. 2.1.1. Sifting process The fundamental step of the empirical mode decomposition (EMD) is the next iterative sifting process: 1. Consider a signal x(t) with M maxima and L minima. The sifting process starts with identifying the extrema of the signal, x(t), given 1 by the sets Emax ¼ fxmax ðt1 Þ; xmax ðt2 Þ; xmax ðt3 Þ; .; xmax ðtM Þg,
123
Fig. 1. Flowchart of the EMD method.
124
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
procedure. The performance of the EMD algorithm is sensitive to the interpolation procedure. The sifting process is defined for continuous signals which means, that the performance of the EMD depends on the sampling rate (Huang et al., 1998). The dependence of the EMD algorithm on these factors makes precludes a general, unique framework for EMD. Thus defining a function space for the EMD algorithm, is an ill posed problem, making it difficult to construct an analytical description of the EMD. However, empirically, the EMD has been shown to be effective in extracting relevant components in a variety of applications involving non-stationary signals. 2.2. HilberteHuang transform Consider the analytic representation of a signal given by: jfðtÞ
Zx ðtÞbxðtÞ þ jb x ðtÞ ¼ AðtÞe
ZN N
xðt sÞ
ps
(1)
ds ;
xðtÞ ¼ Re½Zx ðtÞ ¼ AðtÞcosðfðtÞÞ
(2)
The operation denoted in Eq. (2) is merely the convolution of x(t) with 1/t or it can be thought of as a phase shift of the signal by p/2.
(3)
we can observe that this canonical representation of x(t) corresponds to a signal varying simultaneously in amplitude and phase all the time (exactly like each IMF from the EMD). Based on this representation, Ville (1948) proposed to calculate the instantaneous frequency associated to the instantaneous phase 4(t). From Eq. (1), the instantaneous frequency is defined by:
uinst ðtÞ ¼
where A(t) and 4(t) are the envelope and the instantaneous phase of the signal, respectively. In this equation b x ðtÞ is the Hilbert transform of the signal x(t), which is defined by Delprat et al. (1992):
b x ðtÞ ¼
As a result, x(t) and b x ðtÞ are said to be in quadrature. This complex or analytic signal is completely characterized by their amplitude A(t) and their phase 4(t) with values in the interval [0, 2p), forming the so-called canonical pair. Considering that
dfðtÞ dt
(4)
It is clear that the analytic signal Zx(t) associated with x(t) has the same amplitude and frequency range with x(t), it also comprises the phase information of the original signal x(t), which can be used to compute the instantaneous amplitude and to analyze the frequency performance of x(t). Based on this fact, we can construct the analytic signal corresponding to each IMF using the Hilbert transform, defined in Eq. (2). The combination of the EMD applied to x(t) to generate IMFs, and the Hilbert transform of each IMF is called the HilberteHuang Transform (HHT). The different IMFs resulting from the EMD algorithm are orthogonal to each other. The IMFs, thus represent different time-scales of oscillations, which form a set of bases
Original Signal: Forsmark, Case 4, LPRM 22 80 60 40
0
5
10
Instantaneous Frequency (IF)
15
IF Histrograms 4
IMF 1 0
0
5
10
2
IMF 2
0 -2
0
5
10
0
5
10
2
15
IMF 4
0 -2
0
5
10
IF [Hz]
-5
IF [Hz]
IMF 3
IMF 5 -2
0
5
10
Time [s]
15
IF [Hz]
2
0
5
10
15
0
1
4
0.5
2
0
5
10
15
0
1
10
0.5
5
0
0
5
10
15
0
1
20
0.5
10
0
15
0
2
0
15
5 0
5
0
15
IF [Hz]
-1
IF [Hz]
1
0
5
10
15
0
1
40
0.5
20
0
0
5
10
15
Time [s] Fig. 2. The EMD-HHS algorithm applied to a real signal from a BWR.
0
0
1
2
3
4
5
6
7
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
0
0.2
0.4
0.6
0.8
1
Time [s]
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
functions. This implies that there is no redundancy in the information contained in the different IMFs. Using this property, a distribution can be constructed from the instantaneous frequencies of each of the IMFs. This distribution is called the HilberteHuang Spectrum (HHS) (Huang et al., 1998). Since the instantaneous frequency of the EMD-HHS approach is not based on the Fourier transform, the time-frequency resolution is not limited by uncertainty. It is worth mentioning that the EMD can be viewed as a reconstruction problem where each IMF contributes to the signal reconstruction. It is also clear that the reconstruction error depends on the levels of the EMD: if the EMD is obtained in all the possible levels (exhaustive problem), the reconstructed signal is practically the original one. In our case, this reconstruction situation is irrelevant because we are only interested in a specific IMF: the one associated to the instability. Fig. 2 plots the first 5 IMFs of the EMD decomposition of real signal obtained from a BWR during 15 s. This BWR signal corresponds to the LPRM number 22 from the case 4 of the Forsmark benchmark (Verdu et al., 2001). This case is presented in detail in Section 4 of this work. We observe that the IMFs have a zero mean, symmetric envelopes (AMeFM signals), and the number of extreme is at most one more than the number of zero crossings. The frequency of an IMF is directly proportional to the number of zero crossings in one time duration of IMFs. Also we can observe perfectly the frequencies tracking instantaneous frequency (IF) associated to each IMF. This IF permits to establish which level of the EMD
125
contains the mode associated to the BWR instability. Based on this, a simple method based on the distribution (histogram) of the IF is implemented. Both IMFs and IFs determine the mode or modes associated to the instability phenomena. In the presented signal, this instability problem can be clearly observed at the IMF 3 of Fig. 2. In this figure, the frame (in red) shows the selected IMF and the corresponding IF and its distribution (histogram). Comparing this IF with those found in other levels of the EMD, it is clear that the level 3 (IMF 3) is the right one to be considered. 3. Decay Ratio estimation using EMDeHHT In this work, the EMDeHHT technique to the BWR signals was applied, following the next methodology: 1. The considered signal (APRM or LPRM) obtained from the BWR is segmented in time windows of 15 s of duration. 2. Each segmented signal (APRM or LPRM) is analyzed (decomposed) using the EMD algorithm, described in Section 2.1, obtaining in this way the corresponding IMFs. It is worth mentioning that the APRM or LPRM signals are not being processed before. For instance, to remove the signal trend, it is not necessary to filter the BWR signal before, due that this information is contained in the last scale computed from the EMD (Fig. 3). It is also worth mentioning that the trend does not represent an IMF. Another important feature of the EMD is that
Original Signal 10 0 -10
0
5
10
15
20
25
30
Instantaneous Frequency (IF)
35
IF Histrograms
0.05
40
IF [Hz]
IMF 1 0 -0.05
0
5
10
15
20
0.5 20 0
25
0
5
10
15
0 0.4
0.45
0.5
0.55
0.6
0.35
0.4
0.45
0.5
0.1
0.2
0.3
0.4
-3
IF [Hz]
5
x 10
IMF 2 0 -5
0
5
10
15
20
10 0
25
20
0.5
0
5
10
15
0 0.3
-3
IF [Hz]
5
x 10
IMF 3 0 -5
0
5
10
15
20
20 0
25
40
0.5
0
5
10
15
0
0
-3
IF [Hz]
2
x 10
IMF 4 0 -2
0
5
10
15
20
25
100
0.5
50 0
0
5
10
15
0
0
0.05
0.1
0
0.05
0.1
0.15
0.2
0.25
0.3
0.15
0.2
0.25
0.3
-3
1
x 10
IF [Hz]
Trend 0 -1
0
5
10
15
Time [s]
20
25
100
0.5
50 0
0
5
10
15
Time [s] Fig. 3. EMD-HHS method applied to a synthetic signal.
0
IF [Hz]
126
3. 4.
5. 6.
7.
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
anytime the analyzed signal should contain an additive and high-frequency noise, it will be isolated in the first scale (IMF 1). Thus completely eliminating it from the signal and the other IMFs (including our IMF of interest) permitting a better analysis. The Hilbert transform for each IMF is computed in order to get the instantaneous frequencies contained in each IMF. When tracking these frequencies, it is possible to get the mode or modes associated to instability processes. In this way, only the mode (IMF) or modes associated to the BWRs instability are considered for further processing. The autocorrelation function (AFC) of the considered IMF is calculated. Based on these ACF values, an instability parameter is estimated. This parameter is determined as the Decay Ratio (DR) but unlike it, is calculated considering two consecutive maxima from the ACF magnitude of the analyzed IMF. This consideration (DR estimation) was inferred based on the analysis of a synthetic signal (see below) and the discussion on the Forsmark stability benchmark (see Section 4). Finally, the mean and variance of the DR are calculated averaging the estimated DR obtained in each window of 15 s along time.
Before taking the challenge of real signals, first we use a synthetic signal to validate our method. For this we considered three different oscillating modes driving a second order system which impulse response is given by:
Table 1 Estimated DR and frequency of the Forsmark stability benchmark, Case 4. Power sensor
IMF
Decay Ratio DR (mean)
DR (standard deviation)
Frequency (Hz)
C4_APRM_1 C4_LPRM_1 C4_LPRM_2 C4_LPRM_3 C4_LPRM_4 C4_LPRM_5 C4_LPRM_6 C4_LPRM_7 C4_LPRM_8 C4_LPRM_9 C4_LPRM_10 C4_LPRM_11 C4_LPRM_12 C4_LPRM_13 C4_LPRM_14 C4_LPRM_15 C4_LPRM_16 C4_LPRM_17 C4_LPRM_18 C4_LPRM_19 C4_LPRM_20 C4_LPRM_21 C4_LPRM_22
2 2 2 2 2 3 3 3 2 3 2 3 3 3 3 3 2 3 2 3 3 3 3
0.7143 0.7843 0.8162 0.8413 0.7272 0.8004 0.8068 0.7593 0.6660 0.7635 0.5217 0.7652 0.7001 0.7519 0.7410 0.6985 0.7837 0.7828 0.7046 0.6610 0.6654 0.7880 0.7014
0.1761 0.1244 0.1132 0.1057 0.1842 0.1220 0.0907 0.1559 0.1436 0.1289 0.1847 0.1415 0.1347 0.1694 0.1802 0.1877 0.1145 0.1192 0.1503 0.1219 0.1689 0.1596 0.1174
0.4718 0.5031 0.4923 0.5024 0.4905 0.5041 0.5299 0.5091 0.4951 0.5174 0.5217 0.4793 0.4873 0.5031 0.5150 0.5178 0.5103 0.4819 0.4824 0.4864 0.4873 0.5270 0.5017
Fig. 4. Decay Ratio estimation of the IMF 1 from the simulation considering a synthetic signal.
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 hðtÞ ¼ A0 ezun t cos 1 z un t
(5)
where un is the natural resonant angular frequency, A0 is the magnitude of this impulse response, and z is the damping parameter. Here the oscillation term is represented by the cosine function, and the damping effect is modeled using an exponentially decaying factor. The Decay Ratio (DR) is calculated by:
0
1
2pz C B DR ¼ exp@ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiA 2 1z
(6)
The input signal applied to this system is established as:
xðtÞ ¼ 5 cosð2pf0 tÞ þ cosð2p1:2f0 tÞ þ cosð2p0:35f0 tÞ
(7)
with f0 ¼ 0.5 Hz. It is clearly that the signal in Eq. (7) will be in resonance with the system, which the impulse response is defined in Eq. (5), at the frequency u0 ¼ 2pf0. The idea behind the creation of this synthetic signal is to prove that the proposed method is well adapted to estimate the DR in a second order system (like in a BWR) considering that the output is composed by different modes
127
attenuated or reinforced by the system along time. Indeed, a real BWR signal from the LPRM or APRM, the density wave phenomenon is presented like a mode around 0.5 Hz growing in amplitude along time. Harmonics of this mode can be added reinforcing its presence. For this simulation a sampling frequency of 12.5 Hz was considered. The input signal given in Eq. (7) was introduced in the system defined by Eq. (5). To obtain the corresponding IMFs and IFs plotted in Fig. 3 the proposed EMDeHHT method was applied to the output signal of this system. Two important remarks can be pointed out: We can observe that the interesting mode can be tracked along time. In this case, corresponding to the first IMF. The IF associated to this IMF is clearly showed in the middle (at the top) of Fig. 3. Values of this IF are confirmed by the histogram showed in the third column of Fig. 3. It is worth mentioning that the ACF obtained from the analyzed IMF vanishes quickly in time after 15 s. Based on the fact that the proposed figure of merit is estimated from this autocorrelation function, this means that there is not relevant information in the LPRM signals (directly on the IMF) after this time window. This is a very important point to consider in order to apply this
Fig. 5. IMFs and IFs of a specific LPRM for Forsmark case 4.
128
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
method on BWRs real signals. It is also important to mention that this result is also comparable to our simulation case. Following with the proposed methodology, we computed the ACF of the first IMF of the output signal. We obtain the ACF magnitude and we detect the maxima contained in this function. Fig. 4 presented three different cases corresponding to DR values of 0.5, 0.8 and 1. In Eq. (5) the respective damping ratios (i.e. z ¼ 0.1097, 0.0361, 0.0152) are considered for simulation proposals. Maxima of the ACF magnitude are calculated and represented by the red marks (*) in Fig. 4. Again two important remarks can be pointed out: Considering the classical DR (ratio of third and first peaks), in the case of DR equal to 0.5 and 1, values are slightly (less than 10%) underestimated. In the case of the DR equal to 0.8, the estimated DR is practically the real value considered in the simulation. We decided to use as figure of merit also the DR but measured over the two first peaks of the ACF magnitude of the considered IMF. This decision was based on the analysis of several real signals from BWRs, where in certain cases the third peak sometimes is not well calculated (sometimes even lost when the IMF is quite small in amplitude) from the EMDeHHT algorithm. In the case of this simulation, we can observe that for DR equal to 0.5 and 0.8 is slightly (less than 10%) overestimated, for a DR equal to 1 is practically the real considered value. This decision
about the DR and its impact in real signals from BRWs will be discussed in the Section 4 of this work.
4. Experiments and results The methodology described in the previous section, was applied to BWRs real signals. Among those, we include two complete cases from the Forsmark stability benchmark corresponding to cases 4 and 6. It is worth mentioning that all cases from the Forsmark stability benchmark were studied. However, for reasons of space only the most challenging cases are presented. The potential of this proposed method to predict an instability event is to implement a real-time stability monitor in a BWR. In our implementation, a modified version of the EMD code provided in Rilling and Goncalves (2008) was used. 4.1. Cases of the Forsmark benchmark The EMDeHHT method was applied to neutronic signals of Forsmark stability benchmark (Verdu et al., 2001). The benchmark is based on data from several measurements performed in the Swedish BWR reactor Forsmark 1 and 2, in the period 1989 to 1997. We discuss the results obtained specially from two cases: 4 and 6. Case 4. This case contains a mixture between a global oscillation mode and a regional (half core) oscillation. This event corresponds
Fig. 6. IF and DR estimation for Forsmark case 4, APRM 1 signal.
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
Forsmark Signal. Case 4. APRM 1 150
100
50
0
0
0.1
0.2
0.3
0.4
0.5
0.6
Instantaneous Frequency (IF)
0.7
0.8
[Hz]
Fig. 7. Histogram of the IF for Forsmark case 4, APRM 1 signal.
to a situation where the neutronic power reactor suffers abnormal and apparently unstable oscillations. The C4_APRM and C4_LPRM_x signals correspond to average power range monitor (APRM) and local power range monitor (LPRM) registers respectively, during the instability event. Results
129
of the proposed method are shown in Table 1. These results agree with those obtained by different methodologies developed by participants and included in the final report of the Forsmark stability benchmark (Verdu et al., 2001). In Fig. 5, a specific LPRM (7) is considered for the EMD-HSS algorithm. In this figure, the IMFs, IFs associated to these IMFs, and the histograms from the instantaneous frequencies are plotted. The analyzed segment of 15 s is illustrated in the top of the figure and marked by the box (in green) in the complete record from the LPRM. It is clear that the first IMF contains only the information that we will considers like noise which contents is distributed along all frequencies as shown by the histogram of the corresponding IF and, applying this algorithm we can consider an automatic denoising of the analyzed signal. IMFs 4 and 5 showed that the frequencies contained in these modes are the tendency towards low part of the spectrum and they are not considered in the analysis. However IMFs 2 and 3 results show the interesting modes to consider for further processing. Histograms of the corresponding IMF 2 and 3 show both modes containing information about stability (frequency around 0.5 Hz), nevertheless the IMF 3 clearly shows the frequency around 0.5 Hz more persistently. Because that, this is the chosen mode to further processing to get the DR parameter. In Table 1 it is also included the number of IMF which is processing in order to get the values of the DR and IF. So having decided which mode (IMF) to be processed, we obtain the ACF of this mode and we get the proposed DR considering only the two consecutives peaks from the magnitude of the ACF. Fig. 6 shows the application of this algorithm to a specific signal, in this
Fig. 8. IF and DR estimation for Forsmark case 4, LPRM 1 signal.
130
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
Decay Ratio. Forsmark Case 4. 1Level 1 Decay Ratio. Forsmark Case 4. Level 1.2 lprm1-pos23 lprm17-pos9 lprm11-pos11 lprm21-pos31
1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0
50
100
150
200
250
300
350
Time [s]
Time [s] Fig. 9. DR estimation along time for Forsmark case 4.
case the APRM 1. In the top of this figure the complete original is plotted. In this figure each analyzed segment of the signal is marked by a green box and plotted immediately below. In the middle of Fig. 6 the corresponding instantaneous frequency (IF) for the chosen IMFs are shown. Two last parts (below) of this figure show the obtained ACF magnitude and the corresponding estimated DR from this ACF. Estimated DR is plotted along time and computed its corresponding mean and standard deviation. It is worth mentioning that this method permits two important things along time: tracking the frequency associated to an instability and corresponding Decay Ratio. In Fig. 6 we can observe that the obtained DR corresponds to those reported in the benchmark by the different applied methodologies (between 0.459 and 0.9, but mostly between 0.7 and 0.85). In our method only two segments of 15 s are under a DR of 0.5 (see Fig. 6), the most part of results are in the
interval going from 0.7 to 0.85 (see Fig. 6), in perfect agree with the benchmark. In Fig. 7 the distribution of the IF from the analyzed IMFs is plotted. We can observe a Gaussian distribution of this IF for this APRM. This distribution is centered in the frequency 0.48 Hz completely in agree with those presented in the benchmark (between 0.48 and 0.51). Fig. 8 corresponds to the signal from the LPRM 1. We can observe how this method perfectly tracks the IF from the IMF. The corresponding DR along time varies from 0.7 to 0.9, only three analyzed segments give a DR lower than 0.6. Most part of estimated DR is between 0.8 and 0.9 in completely agree with results reported in the benchmark. It is clear that the observed mode (IMF 2) indicates a core more unstable in this part of the reactor than the corresponding global measure (APRM 1). This Gaussian distribution is
Table 2 Estimated DR and frequency of the Forsmark stability benchmark, Case 6.1.
Table 3 Estimated DR and frequency of the Forsmark stability benchmark, Case 6.2.
Power sensor
IMF
Decay Ratio DR (mean)
DR (standard deviation)
Frequency (Hz)
Power sensor
IMF
Decay Ratio DR (mean)
DR (standard deviation)
Frequency (Hz)
C6_APRM_1 C6_LPRM_1 C6_LPRM_2 C6_LPRM_3 C6_LPRM_4 C6_LPRM_5 C6_LPRM_6 C6_LPRM_7 C6_LPRM_8 C6_LPRM_9 C6_LPRM_10 C6_LPRM_11 C6_LPRM_12 C6_LPRM_13 C6_LPRM_14 C6_LPRM_15 C6_LPRM_16 C6_LPRM_17 C6_LPRM_18
3 3 3 3 3 3 3 3 3 4 4 4 3 4 3 3 3 4 3
0.6722 0.6183 0.6492 0.5981 0.6288 0.5797 0.6906 0.6148 0.6852 0.6264 0.6466 0.7217 0.5687 0.6295 0.6087 0.5508 0.6242 0.6233 0.6323
0.1078 0.1323 0.1285 0.1555 0.1015 0.1522 0.1254 0.1327 0.1608 0.1573 0.1385 0.1069 0.1592 0.1484 0.1178 0.1399 0.1379 0.1981 0.1416
0.4932 0.5662 0.5016 0.5093 0.5036 0.5079 0.5127 0.4718 0.4639 0.4450 0.4430 0.4278 0.4934 0.4297 0.5499 0.4940 0.5207 0.4739 0.4943
C6_APRM_2 C6_LPRM_1 C6_LPRM_2 C6_LPRM_3 C6_LPRM_4 C6_LPRM_5 C6_LPRM_6 C6_LPRM_7 C6_LPRM_8 C6_LPRM_9 C6_LPRM_10 C6_LPRM_11 C6_LPRM_12 C6_LPRM_13 C6_LPRM_14 C6_LPRM_15 C6_LPRM_16 C6_LPRM_17 C6_LPRM_18
3 3 3 3 3 3 2 3 2 3 3 3 3 3 3 3 3 4 3
0.8682 0.5937 0.6790 0.8529 0.6596 0.8882 0.8679 0.8422 0.8609 0.5445 0.6104 0.6709 0.7714 0.7059 0.7754 0.8425 0.7942 0.6842 0.6100
0.0589 0.1230 0.1776 0.0862 0.2000 0.0544 0.1320 0.1099 0.1089 0.1309 0.1158 0.1493 0.1125 0.1316 0.1281 0.1238 0.1187 0.1312 0.1727
0.4986 0.5494 0.4810 0.5226 0.5020 0.5065 0.5024 0.5065 0.5162 0.4757 0.4734 0.5051 0.5240 0.5318 0.5377 0.5223 0.5214 0.4459 0.5798
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
131
centered at 0.495 Hz and it corresponds to reported values (between 0.48 and 0.5 Hz) in the benchmark. All values (DRs and frequencies) corresponding to the reported LPRMs in the Forsmark benchmark are included in Table 1. To end this case, we pointed out the discussion on the benchmark conclusions, especially about these two remarks (Verdu et al., 2001):
performed close to each other, both in time and in operating conditions. For test 1 (6.1): power at 64.4% and core flow of 4 416 kg/s, for test 2 (6.2): power at 63.3% and core flow of 4 298 kg/s. The C6_APRM and C6_LPRM_x signals correspond to an average power range monitor (APRM) and local power range monitor (LPRM) registers respectively, during the instability event considered.
There is a half of the reactor (locations 23 and 9) where the DR is high and the other half (locations 31 and 11) where the DR is lower. The upper part of the reactor seems to be more stable than the lower part. Spectral analysis of the signals indicates that there is a phase shift between the LPRM at radial locations 23 and 11, and locations 23 and 31, but the out-of-phase oscillation is not totally developed.
Case 6.1. Results of the proposed method are shown in Table 2. These results are completely agreed with those reported in the Forsmark stability benchmark (Verdu et al., 2001). However, we focused the discussion on the second test (case 6.2).
In Fig. 9, the estimated DRs along time for the LPRMs 1, 7, 11 and 21 (positions 23, 9, 11 and 31 respectively) at Level 1 are plotted. In our opinion (in agree with remarks above) the out-of-phase oscillation is not completely developed. In this figure, we can observe that only in certain time regions we have this regional oscillation but it is not presented in the complete analyzed signals. Case 6. This test shows local (channel) oscillations. The data contain APRM and LPRM signals from two tests that were
Case 6.2. Results of the proposed method are shown in Table 3. Like the other shown cases, results are in complete agree with those presented in the benchmark. In Fig. 10 a signal segment from LPRM 15 is plotted. Like in Case 4, we can observe that the IMF 3 contains the frequency associated to instability event (around 0.5 Hz). Almost in all APRM and LPRMs for this case, the IF is tracked on the IMF 3. Fig. 11 shows an analyzed segment from the APRM 2. In the middle of the figure, the tracking of the IF associated to the IMF 3 along time is shown. We can clearly observe estimated values of the IF around 0.5 Hz. In the bottom of Fig. 11, estimated DR is plotted along time. For this APRM, we can see almost all DR values varying
Fig. 10. IMFs and IFs of a specific LPRM for Forsmark case 6.2.
132
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
Fig. 11. IF and DR estimation for Forsmark case 6.2, APRM 2 signal.
Decay Decay Ratio. Forsmark Level 1 Ratio. Forsmark CaseCase 6.2. Level6.2. 1 1
0.9
0.8
0.7
0.6
0.5
lprm1-pos23 0.4
lprm17-pos29 lprm7-pos6 lprm15-pos24
0.3
0.2 0
50
100
150
200
250
Time [s] Time [s]
Fig. 12. DR estimation along time for Forsmark case 6.2.
300
350
A. Prieto-Guerrero, G. Espinosa-Paredes / Progress in Nuclear Energy 71 (2014) 122e133
between 0.8 and 0.95. This indicates an almost unstable event. This situation corresponds to those reported in the benchmark. One important conclusion from this case in the benchmark indicates that one half of the reactor is more stable than the other one. To identify this situation we plotted simultaneously in Fig. 12 the estimated DR obtained from the LPRMs 1, 17, 7 and 15 (positions in the reactor 23, 29, 6 and 24 respectively) at Level 1. In this figure we can observe, practically all the time, how a half of the reactor is unstable (LPRM 7 and 15, in solid lines) presenting oscillations varying almost all between 0.8 and 0.95. The other half of the reactor is stable, with estimated DRs (marked with dotted lines) varying almost all between 0.4 and 0.7. This is a very important result that allows visualizing the potential of the proposed method, capable of determining in real time, not only possible unstable events but also possible out-of-phase oscillations in the reactor core.
5. Conclusions In this work, we propose a novel approach to stability analysis of boiling water reactors (BWRs) based on empirical mode decomposition (EMD) and the HilberteHuang transform (HHT). The method consists of obtaining the different intrinsic modes functions (IMF) associated to the EMD. Applying the HHT it is possible to get from each IMF, the associated instantaneous frequency (IF) along time. This instantaneous frequency allows to decide which IMF will be associated to a possible instability event. In this manner, it is possible to track this instability situation. An estimated Decay Ratio (DR) is obtained time along, considering the magnitude of the autocorrelation function (ACF) of followed IMF. Some important conclusions can be delivered. First, the method is relatively simple to implement and it does not represent a high computational complexity. Second, results obtained from the implementation of this method on the Forsmark stability benchmark, are completely in concordance with those reported by the different methodologies developed by the participants of this benchmark. For Cases 4 and 6, included and discussed in this work, we can conclude that the proposed method clearly contributes on the fact to detect events of instability along time inclusive detecting possible cases of out-of-phase oscillations. Conclusions reported in the benchmark (Verdu et al., 2001) were corroborated by the proposed method in this work. This method, with the proposed estimation for the DR, also allows predicting the instability event through time before the oscillation going to a critical state. References Cai, J., Ishiwatari, Y., Ikejiri, S., Oka, Y., 2009. Thermal and stability considerations for a supercritical water-cooled fast reactor with downward-flow channels during power-raising phase of plant startup. Nucl. Eng. Design 239 (4), 665e679. Delprat, N., Escudie, B., Guillemain, P., Kronland-Martinet, R., Tchamitchian, P., Torresani, B., 1992. Asymptotic wavelet and Gabor analysis: extraction of instantaneous frequencies. IEEE Trans. Info. Theory 38, 644e664. Espinosa-Paredes, G., Prieto-Guerrero, A., Núñez-Carrera, A., Amador-García, R., 2005. Wavelet-based method for instability analysis in boiling water reactors. Nucl. Technol. 151, 250e260.
133
Espinosa-Paredes, G., Nuñez-Carrera, A., Prieto-Guerrero, A., Ceceñas, M., 2007. Wavelet approach for analysis of neutronic power using data of Ringhals stability benchmark. Nucl. Eng. Design 237, 1009e1015. Guido, G., Converti, J., Clause, A., 1997. Density wave oscillations in parallel channels, an analytical approach. Nucl. Eng. Design 125, 121e139. Huang, N.E., Shen, Z., Long, S.R., Wu, M.C., Shih, H.H., Zheng, Q., Yen, N.C., Tung, C.C., Liu, H.H., 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Math. Phys. Eng. Sci. 454, 679e699. Huang, N.E., Wu, Z., 2008. A review on HilberteHuang transform: method and its applications to geophysical studies. Rev. Geophys. 46, 1e23. Khaldi, K., Boudraa, A.O., Turki, M., Chonavel, T., Samaali, I., 2009. Audio Encoding Based on the Empirical Mode Decomposition. EUSIPCO, Glasgow. Lahey, R.T., Podoswski, M.Z., 1989. On the analysis of various instabilities in twophase flow. Multiph. Sci. Technol. 4, 183e371. Lin, S.L., Tung, P.C., Huang, N.E., 2009. Data analysis using a combination of independent component analysis and empirical mode decomposition. Phys. Rev. E. Stat. Nonlin. Soft. Matter. Phys. 79 (6 Pt 2), 066705. March-Leuba, J., 1986. A reduced-order model of boiling water reactor linear dynamics. Nucl. Technol. 75, 15e22. March-Leuba, J., 1990. LAPUR Benchmark against In-phase and Out-of-phase Stability Test. NUREG/CR-5605, ORNL/TM-116210. March-Leuba, J., Blakeman, E.D., 1991. A mechanism for out-of-phase power instabilities in boiling water reactors. Nucl. Sci. Eng. 107, 173e179. March-Leuba, J., Rey, J.M., 1993. Coupled thermal-hydraulic neutronic stabilities, in boiling water reactor a review of the state of the art. Nucl. Eng. Design 145, 97e111. Montesinos, M.E., Muñoz-Cobo, J.L., Pérez, C., 2003. HilberteHuang analysis of BWR neutron detector signals: application to DR calculation and to corrupted signal analysis. Ann. Nucl. Energy 30, 715e727. Muñoz-Cobo, J.L., Chiva, S., Sekhri, A., 2004. A reduced order model of BWR dynamics with subcooled boiling and modal kinetics: application to out of phase oscillations. Ann. Nucl. Energy 31, 1135e1162. Navarro-Esbri, J., Ginestar, D., Verdu, G., 2003. Time dependence of linear stability parameters of a BWR. Prog. Nucl. Energy 43, 187e194. Peng, S.J., Podowski, M.Z., Lahey, R.T., Becker, M., 1984. NUFREQ-NP: a computer code for the instability analysis of boiling water nuclear reactors. Nucl. Sci. Eng. 88, 404e411. Peng, Z.K., Tseb, P.W., Chu, F.L., 2005. An improved HilberteHuang transform and its application in vibration signal analysis. J. Sound Vib. 286, 187e205. Podowski, M.Z., Pinheiro, M., 1997. Modeling and numerical simulation of oscillatory two phase flow with application to boiling water nuclear reactor. Nucl. Eng. Design 177, 179e188. Prieto-Guerrero, A., Espinosa-Paredes, G., 2008. Decay Ratio estimation of BWR signals based on wavelet ridges. Nucl. Sci. Eng. 160, 302e317. Rilling, G., Goncalves, P., 2008. EMD Toolbox for Matlab. http://perso.ens-lyon.fr/ patrick.flandrin/emd.html. Saha, P., Zuber, N., 1978. An analytical study of the thermally induced two-phase flow instabilities including the effects of thermal non-equilibrium. Int. J. Heat Mass Transf. 21, 415e426. Sunde, C., Pázsit, I., 2007. Wavelet techniques for the determination of the Decay Ratio in boiling water reactors. Kerntechnik 72, 7e19. Torrres-Fernández, J.E., Espinosa-Paredes, G., Prieto-Guerrero, A., 2012. Decay Ratio estimation in BWRs using the exponential time-scale representation. Ann. Nucl. Energy 49, 143e153. Torres-Fernández, J.E., Prieto-Guerrero, A., Espinosa-Paredes, G., 2010a. Decay Ratio estimation using the instantaneous autocorrelation function. Nucl. Eng. Design 240, 3559e3570. Torres-Fernández, J.E., Prieto-Guerrero, A., Espinosa-Paredes, G., 2010b. Decay Ratio estimation based on time-frequency representations. Ann. Nucl. Energy 37, 93e102. Turso, J.A., March-Leuba, J., Edwards, M., 1997. A modal-based reduced-order model of BWR out-of-phase instabilities. Ann. Nucl. Energy 12, 921e934. Uehiro, M., Rao, Y.F., Fukuda, K., 1996. Linear stability analysis on instabilities of inphase and out-of-phase modes in boiling water reactors. J. Nucl. Sci. Technol. 33, 628e635. Verdu, G., Ginestar, D., Muñoz-Cobo, J.L., Navarro-Esbri, J., Palomo, M.J., Lansaker, P., Conde, J.M., Recio, M., Sartori, E., 2001. Forsmark 1&2 Stability Benchmark. Time Series Analysis Methods for Oscillations during BWR Operation. Final Report. NEA/NSC/DOC, p. 2. Ville, J., 1948. Theorie et application de la notion de signal analytique. Cables Trans. A (1), 61e74. Wu, M.C., Huang, N.E., 2009. Biomedical Data Processing Using HHT: A Review. Chapter 16. In: Advanced Biosignal Processing. Springer-Verlag, pp. 335e352.