A Novel Nonlinear BWR Stability Indicator Based on the Sample Entropy

BWRs are thus far the simplest energy systems to transformfission energy into electrical power.However, there are stillmany aspects in their operation that, under certain conditions, may induce BWR unstable behavior.The default indicator to study BWR unstable behavior is the Decay Ratio (DR). However, due to the fact that BWRs show very complex responses under instability and responses that may even be chaotic, the DRmight not be a suitable choice to rely on to accommodate for such intricate behavior. In this work a novel methodology based on the Sample entropy (SampEn) and the noise-assisted multivariate empirical mode decomposition (NA-MEMD) is introduced. Such methodology was developed thinking for a real time-implementation of a stability monitor. The proposed methodology was tested with a set of signals that stem from several nuclear power plants in operation today that have experienced in the past unstable events, each one of a different nature.


Introduction
The issue of BWR instability has been a topic of concern of many scientific and technological works for more than forty years.BWR unstable events are quite rare and may happen during BWR start up or during transients that may interfere with the conventional operating region of the reactor.BWR instability phenomena are strongly linked with the interaction between thermal hydraulic and neutron kinetics processes such interaction under certain conditions and may induce undesirable BWR behavior that could jeopardize BWR operation.Thus, a reliable prediction for the onset of BWR instability is of utter importance for BWR core safety [1].D' Auria [2] presented all BWR unstable events recorded up to 1999, which were due to various incidents (e.g., [3,4]), while others were intentionally induced for experimental purposes (e.g., [5,6]).Unstable oscillations in the neutron flux were recorded during these unstable events through the use of the electronic instrumentation available in the reactor.After the first events of instability were observed, the respective authorities fostered the growth of research projects to study the physical phenomena involved in BWR instability.
The most common instability type for commercial BWRs is the density-wave instability [7][8][9].This type of instability can be described as follows: given a flow perturbation, a wave of voids travels upwards through the channel producing a pressure drop (Δ = Δ 1 + Δ 2 , is approximately constant, where Δ 1 is the pressure drop in the single phase in the liquid region and Δ 2 is the pressure drop in the two-phase flow region) which is delayed with respect to the original perturbation.An increase in flow might induce an increase in pressure drop and a negative feedback that reduces the flow perturbation.The density wave phenomena delays this feedback and, at some frequency, the delay is equivalent to a 180 ∘ phase lag; so at this frequency the pressure drop feedback is positive and if the gain is large enough, the channel flow becomes unstable and oscillates at that frequency highly concentrated around 0.5 Hz, where the associated resonance frequency is related with the fluid transit time along the heated channel, which could differ from 2 seconds, i.e., approximately 1 Hz.
In practice, the most commonly used indicator to study BWR stability due to wave oscillation is the decay ratio (DR), which is computed from the impulse response function obtained from a model fitting a linear second order system, accommodating the core behavior of a BWR.The DR is the current output indicator of most, if not all, of the stability monitors implemented in NPPs today [10].The DR as a feasible BWR stability measure has been widely accepted but it has been observed that a BWR working at an operating point with a small DR can be close to instability [11].Sometimes, the DR jumps from the stable to the far unstable region [12].According to these works, the DR might not be a reliable option, under certain conditions.Besides, the need for linear and stationary signals might be a handicap for DR estimation.Therefore, it is necessary to explore new methodologies and indicators adapted to accommodate at their theoretical framework nonlinear and nonstationary behavior in order to study BWR unstable phenomena with as much realism as possible.
We are trying to convey in this work this idea, for this we explore the Sample Entropy (SampEn) linked to the noiseassisted multivariate empirical mode decomposition (NA-MEMD) to infer whether BWR signals are stable or not.The SampEn [13] is a measure that provides an index of signal complexity or irregularity of a time series.In this way, the SampEn might work as a possible nonlinear BWR stability indicator.SampEn was in principle developed almost exclusively to analyze physiological time series [13] but its utility has expanded to other domains.For instance, SampEn has been tested with daily weather temperature to measure climate complexity [14].It has been used to measure the complexity of the dynamic reconfiguration of the brain [15] to infer its association with normal aging.
To properly estimate the SampEn from real BWR signals, the NA-MEMD [16] was explored.The NA-MEMD is an algorithm that decomposes non-stationary signals that stem from nonlinear sources.This technique also mitigates the mode mixing phenomena [17] of the default EMD [18] technique.The EMD is the root technique that inspired the development of the NA-MEMD.The NA-MEMD produces a local and fully data-driven decomposition of a studied signal into its fast and slow oscillations.At the end of the procedure, the studied signal can be expressed as a sum of amplitude and frequency modulated (AM-FM) functions called intrinsic mode functions (IMFs), or simply called modes, plus a residue (the trend) of the decomposition.The combination of the NA-MEMD and the Hilbert transform is known as the Hilbert-Huang transform (HHT).The methodology we introduce here is based on the HHT and it computes an indicator linked to BWR stability, in this case the aforementioned SampEn.The NA-MEMD decomposes the analyzed BWR signal into IMFs.One or more of these extracted modes can be associated to the instability problem in BWRs.Through HHT it is possible to get the instantaneous frequency (IF) linked to each IMF.By tracking this IF and the SampEn of the IMF linked to instability, the estimation of the SampEnbased stability indicator is accomplished.The methodology proposed in here is a continuation of two previous works [19,20], in which a Shannon Entropy estimator [21] was used in conjunction with other members of the EMD family to study the stability of artificial and real BWR signals.In the present work the Shannon Entropy estimator was replaced by the SampEn practical formula to infer whether BWR signals are stable or not.
This work is organized as follows: in Section 2, a full review of the SampEn and of the NA-MEMD techniques is given.In Section 3, the methodology to compute the instantaneous frequency and the proposed SampEn novel stability indicator is described.In Section 4, the validation of the methodology presented in this paper is performed doing experiments with real signals taken from the Forsmark [5] and Ringhals [6] stability benchmarks and from a typical BWR power plant during stable and unstable operation.Our major observations regarding our novel methodology based on SampEn and the NA-MEMD are thoroughly discussed in Section 5.

Preliminaries
2.1.The Sample Entropy.The procedure for computing Sam-pEn was introduced by Richman and Moorman [13].In this work, we grant a brief summary of their findings.SampEn is a measure of complexity.Let  = [x 1 , x 2 , ..., x  ] be a time series of length .To compute the complexity of this time series via SampEn, follow the next steps: (1) Build a vector V  with  consecutive data points taken from where  is the length of sequences to be compared, also called the embedding dimension.
(2) For each  define where  varies in the interval (1 ≤ i ≤ N − ).In here,  is the tolerance for accepting matches,  =  × std(), where  is a scaling parameter and std(x) is the standard deviation of .Θ(⋅) is the Heaviside function: and ‖ ⋅ ‖ 1 is the Chebyshev distance, defined as (3)    represents the proportion of V  (i ̸ = j) whose distances to V  are less than .Now, for each  we also define where  +1  represents the proportion corresponding to the dimension of  + 1.    and  +1  have the same mold, but embedding vectors in both cases are defined in different spaces.
(4) Average across all embedding vectors, to obtain (6) and (5) The SampEn is computed as SampEn is the negative natural logarithm of the conditional probability that two sequences similar for  points remain similar at the next point, where self matches are not included in calculating the probability.Thus, a lower value of SampEn indicates more self-similarity (i.e., high order) of the studied time series whereas a higher value of SampEn points to higher complexity of the time series.The calculation of SampEn requires a priori determination of two unknown parameters,  and  (the length of data  is up to the user).The suggested values of  are located in the range of 0.1-0.2times the standard deviation of the studied signal .In this work, we looked for the  value that grants the global SampEn maximum because this maximum value leads to the correct interpretation of signals complexity [22].So, this  of interest is fixed in our simulations at  = 0.1.The value of  can be computed via the estimation of false nearest neighbor [23].In our particular BWR stability discipline,  that was found through false nearest neighbors for our signals was most of the time equal to 2. Therefore,  is fixed at  = 2 for all of our computations.SampEn is, theoretically speaking, a fraction in the interval 0 < SampEn < ∞.However, the next two formulas can be used to find the lower bound and the upper bound of SampEn for fixed values of  and .
The lower bound is computed as The upper bound is computed as For all of our simulations m, r and  are fixed at  = 2,  = 0.1 × std() and  = 300 data points.Thus the lower bound is practically ≈ 0 whereas the upper bound is close to 11.However, the upper bound value was never attained for any simulation.

The Noise-Assisted Multivariate Empirical Mode Decomposition.
The multivariate empirical mode decomposition [16], commonly referred to as the MEMD, is a nonlinear filter to make the classic EMD [18] suitable for processing of multichannel signals.The behavior of the MEMD was analyzed in the presence of white Gaussian noise [20,24] and it was found that the MEMD in essence acts as a dyadic filter bank on each channel of the multivariate input signal; such MEMD property is illustrated in Figure 1 and its steps are given below.
The MEMD has the property of aligning the IMFs or modes from different channels across the same frequency range and the NA-MEMD technique was developed to aid resolve the mode mixing problem observed in the current MEMD, two key issues for real applications.The NA-MEMD technique, which makes use of the quasi-dyadic filter bank properties of MEMD on white noise (Figure 1), is capable of reducing the mode mixing problem for classes of signals where the quasidyadic filter bank structure proves useful (such property applies for BWR signals).The NA-MEMD adds  extra channels of white noise to the input studied signal before being decomposed by the MEMD.The addition of these  extra channels helps to establish a uniformly distributed reference scale which, in turn, results in corresponding modes exhibiting a quasi-dyadic filter bank structure.The extracted IMFs corresponding to the  extra channels of white noise are then discarded yielding a set of IMFs linked with only the original input signal of interest.The NA-MEMD procedure is summarized in the following next steps [20]: Step 1. Create an uncorrelated Gaussian white noise timeseries (q-channels) of the same length as that of the input.
Step 2. Add the noise channels (q-channels) created in Step 1 to the input multivariate (L-channels) signal, obtaining an (L+q)-channel signal.
Step 4. From the resulting (L+q)-variate IMFs, discard the q channels corresponding to the noise, giving a set of L-channel IMFs corresponding to the original signal.However, it should be mentioned that the NA-MEMD for mitigating the mode mixing becomes more useful for signals in which the dyadic filter bank decomposition is relevant.This is the case for the studied BWR signals.

BWR Stability Methodology Based on the SampEn and the NA-MEMD
The methodology to detect stable or unstable states (this methodology works with raw data.So, there is no need to preprocess the analyzed data) in a core of a BWR, based on the Sample entropy and the NA-MEMD, is presented in the next steps: (1) The studied raw signal (APRM or LPRM) obtained from the BWR is segmented in windows of 60 s of duration.
(2) Each segment (of 60 s of time span) is considered as just an independent channel (L=1) added with q independent channels of white Gaussian noise and decomposed through the NA-MEMD (for all of our computer simulations q=2).
(3) After decomposition, the  channels corresponding to the noise are discarded, giving a set of IMFs corresponding to each decomposed studied segment.
(4) The Hilbert transform of each IMF (i.e., the HHT) is computed to get the instantaneous frequency of each IMF in the studied segment.
(5) When tracking these frequencies, it is possible to get the mode (or modes) linked to the density wave oscillation.In this way, only the mode (or modes) associated to BWR instability is considered for further processing whereas the remaining extracted modes are ignored.
(6) The SampEn of the tracked mode of interest (IMF) is computed for  = 300 (the length of the analyzed segment of 60 s),  = 2 and  = 0.1 × std(segment) to attain the global maximum of the SampEn that leads to the correct interpretation of the studied segment complexity [22].(7) The mean, median, and the variance of the SampEn values are computed and averaged along all the studied segments of 60 s of the original signal. A

Results and Discussions
Now we validate the methodology given in Section 3 with the next data sets, where f s is the sampling frequency of each set: (1) Data set from a typical BWR: Hz.

Stable Signals.
The aim of this test is to observe the demeanor of the BWR in steady-state conditions.The results obtained are presented in Figures 2-8. Figure 2 shows the plot of one APRM stable signal of a typical BWR. Figure 3 shows in red dotted line, a studied segment of 60 s of duration of the APRM signal with its power spectral density (PSD) estimate (bottom of figure).Figure 4 shows the first four extracted IMFs of the studied signal, which happens to be noise that is decomposed in a quasidyadic fashion by the NA-MEMD.Figure 5 shows the power spectral density (PSDs) of some IMFs.
Figure 6 shows the tracked IMF (IMF 3) of interest whose instantaneous frequency (IF) oscillates around 0.5 Hz.The computed SampEn of this IMF is 0.9351.Such high SampEn value indicates that the studied signal is BBNCN, so the studied segment is stable and far from an unstable behavior.In this stable case, a mode exist whose IF is close to 0.5 Hz.However, the IMF 3 time series looks distant from a cyclic function we look for when an unstable event triggers in the core.All of the IFs were computed through the Hilbert-Huang transform of the extracted IMFs. Figure 7 shows the computed SampEn along time for all of the studied segments   of 60 s; most of the SampEn estimates have a value slightly greater than 1 (i.e., the studied signal is stable).Table 1 shows the mean(SampEn), median(SampEn) and std(SampEn) of the time series shown in Figure 7. estimates for the studied signals are higher than 1.Thus, the mean SampEn values point to BBNCN.Through SampEn estimation, the proposed stability methodology is able to detect stable signals and classify them in the stable category if their SampEn is greater than 1.

Unstable Signals.
The studied signal in this subsection stems from an instability event.This instability case is registered by 96 LPRM signals.For reasons of space, only the analysis of one signal is presented.Figure 9 shows the plot of the studied signal (from LPRM 2); this signal looks like a cyclic function.Figure 10 shows a studied segment of 60 s (red dotted line) and its PSD estimate (bottom of figure), the segment looks like a noisy cyclic function.Figure 11 shows the first 4 extracted IMFs of the studied segment;   the first 3 IMFs are linked with acquisition noise, whereas IMF 4 looks like a cyclic function (this is in fact the type of waveform commonly associated with density wave unstable events).Figure 12 shows the PSD estimates of the first 6 IMFs, the energetic content of IMF 4 is highly concentrated around 0.5 Hz, and the energetic content of the other IMFs is meaningless next to the energy of IMF 4. Figure 13 shows the time domain series of IMF 4 (which happens to be the mode associated to BWR instability) and its associated IF that oscillates around 0.5 Hz.The IF was computed via HHT.Now, this IF (IF 4) carries physical meaning [7].The SampEn estimate of IMF 4 (IMF or mode linked to instability, see IF in Figure 13 and PSD estimate in Figure 12) is 0.6181; this estimate is a clear indication of signal regularity and low complexity.The SampEn points to a simple function; in this case a cyclic function that in our context is linked with a density wave type of instability.Thus, the studied segment is unstable.Figure 14 shows the SampEn estimates along time for all of the studied segments of the signal of interest; all estimates are smaller than 1.
Table 2 shows the mean(SampEn), median(SampEn), and std(SampEn) of the time series shown in Figure 14.Now, the tracked IF (IF 4) along time carries physical meaning.Thus, the mean(IF) and std(IF) values of all of the studied segments are now provided in this table.
Figure 15 shows a plot of the computed mean SampEn estimates of the 96 unstable signals (each SampEn estimate is a red point).A pattern appears, all of the mean SampEn estimates for the studied signals are smaller than 0.8.Thus, the mean SampEn values point to visually regular signals such as cyclic functions.Through SampEn estimation, we are able to build a binary classifier that separates BBNCN signals (linked to stability) from cyclic ones (linked to density wave instability).One core feature of SampEn is that it accommodates for complex nonlinear and nonstationary data whereas conventional Decay Ratio estimates must assume beforehand that the studied BWR signals behave linearly (an assumption that is false in real life).

Comparison between Stable and Unstable States.
To complete the analysis of the stable/unstable signals through SampEn, Figure 16  There is a distance of ≈ 0.6 between the two states (stable from unstable).So, by fixing a threshold value around 0.9, it is possible to differentiate one state from the other properly.Any SampEn estimate of a BWR signal segment higher than 0.9 is stable whereas any lower SampEn value from this threshold is unstable.difficult case to study by the complexity of the phenomena.For reasons of space, only the analysis of one case will be presented in detail.This case contains a mixture between a global oscillation and a regional oscillation.This particular event corresponds to a situation where the neutronic power reactor suffers unusual unstable problems, i.e., presenting a mix of oscillation modes.The C4 APRM and C4 LPRM x signals correspond to average power monitors (APRM) and local range monitors LPRM.The entire case 4 consists of a total of 23 signals, 22 LPRMs and 1 APRM.Nonetheless, only the analysis of the APRM is presented in this work.Figure 17 shows the studied APRM signal of interest.Figure 18 shows a studied 60 s seconds segment (the dotted red line) that is decomposed through NA-MEMD and its PSD estimate (bottom of figure).Figure 19 shows the first 5 extracted IMFs of the studied segment (as a reminder, the IMFs of the 2 noise extra channels are discarded).Figure 20 shows the PSD estimates of the first 6 IMFs; we highlight that there is a bit of mixing between IMF 4 and IMF 5. Nonetheless IMF 5 is slightly closer to 0.5 Hz than the IMF 4 and has more energy.Figure 21 shows a plot of the IMF 5 of interest that is linked to BWR density wave instability and its associated instantaneous frequency (IF) that oscillates around to 0.5 Hz (we highlight that IMF 5 is slightly mixed with IMF 4 but this is a result due to the complexity of the Forsmark Case 4. Nevertheless, IMF 5 has more energy than IMF 4).The SampEn of this IMF 5 is 0.5443.Thus this segment is unstable (SampEn points to high regularity and low complexity).Figure 22 shows the SampEn estimates of the IMFs or modes of interest along time for all of the studied segments of 60 s.All of the SampEn estimates are smaller than 0.7.So, the studied APRM signal is unstable.
Table 3 shows the mean(SampEn), median(SampEn), and std(SampEn) of the time series shown in Figure 22 (SampEn estimates of IMF 4).The mean(IF) and std(IF) are also given in this table.Figure 23 shows a plot of the computed mean SampEn estimates of the 23 unstable signals (each SampEn estimate is a red point).A pattern appears; all of the mean SampEn estimates for the studied signals are As stated before, only the analysis of one signal (LPRM 1) is detailed in this work.Figure 24 shows a plot of the studied LPRM 1 signal (of level 4). Figure 25 shows a plot of the studied segment (the red dotted line) that is decomposed through NA-MEMD and its PSD estimate (bottom of figure).Thus, the studied segment is unstable).Figure 29 shows the computations of the SampEn along time for all of the studied segments of 60 s.All of the estimates are smaller than 0.8 and located between 0.4-0.6.So, the studied signal is unstable and highly regular like a cyclic function.Table 4 shows the mean(SampEn), median(SampEn), and std(SampEn) of the time series shown in Figure 29.The mean(IF) and std(IF) of all of the studied segments of the targeted LPRM are also given in this table.Figure 30 shows the computed mean(SampEn) values of the 36 LPRMs of level 2 of the Ringhals reactor.It is observed, that all of the SampEn estimates are smaller than 0.8 (The percentage of classification of the studied signals in the unstable category is of 100%).The entire floor 2 is unstable and the SampEn estimates are located in the range 0.5-0.7. Figure 31 shows the computed mean(SampEn) values of the 36 LPRMs of

Conclusions
In this work the Sample Entropy, a measure that provides an index of signal complexity or irregularity of a time series, was In the future, the SampEn might be tested in conjunction with other suitable nonlinear indicators (possibly fractal dimension estimators [25]) to build a robust BWR stability monitor based on 2 nonlinear measurements, where one of them is going to be the SampEn.In general, the SampEn is a suitable and easy to implement BWR nonlinear stability indicator for real time series analysis.

Figure 12 :Figure 13 :
Figure 12: PSD of the first 6 extracted IMFs in the analyzed segment.The PSD of the IMF 4 is highly concentrated around 0.5 Hz.

Figure 14 :
Figure 14: SampEn estimates of the IMFs or modes of interest along time.

FrequencyFigure 20 :Figure 21 :
Figure 20: PSD estimate of the first 6 IMFs extracted in the analyzed segment.

Figure 26 shows the first extracted 5 Figure 22 :Figure 23 :
Figure26shows the first extracted 5 IMFs (by the way, IMF 4 looks like a cyclic function) of the segment.Figure27shows a plot of the PSD estimates of the first 6 IMFs.The PSD estimate of IMF 4 is highly concentrated around 0.5 Hz.Figure28shows the tracked IMF 4 of interest linked to BWR instability, this IMF 4 looks like a cyclic function and its associated IF is almost a line centered around 0.5 Hz.The computed SampEn of this IMF 4 is 0.4544 (the studied IMF 4 is unstable.
[8,9]his high value might indicate that the studied signal could be merely broad band non-coherent noise (BBNCN).This BBNCN, from the BWR point of view is a sign of stable behavior.On the other hand, a low SampEn value (<0.8) points to high predictability and regularity of the studied BWR signal.Such low value might indicate that the studied signal could be a cyclic one.Such waveform is a sign of BWR unstable behavior, due to the presence of a density wave that manifests in the LPRM or APRM signals as a cyclic waveform with a natural frequency that oscillates very close to 0.5 Hz.For forced convection phenomena, in the circulation natural convection the effect of the chimney changes this value[8,9].
high SampEn estimate (>1) points out high irregularity and high unpredictability of the studied BWR signal (LPRM or APRM)
shows a comparison of the mean SampEn Science and Technology of Nuclear Installations of the studied sets (stable set of 208 signals plus a set of 96 unstable signals).The computed SampEn mean values of the stable signals fluctuate very close to 1.2 whereas the SampEn mean values of the unstable signals oscillate very close to 0.6.So the SampEn detects with success stable BWR signals from unstable ones via complexity analysis (of course with the aid of the NA-MEMD to denoise these signals and isolate the mode linked to instability). values
Figure 29: SampEn estimates of the IMFs or modes (IMFs 4) of interest along time.Figure 30: Mean(SampEn) values for the 36 LPRMs of level 2 of the Ringhals BWR unit.used to estimate complexity of a BWR time series for stability analysis.The SampEn was tested in conjunction with the noise-assisted multivariate empirical mode decomposition, a nonlinear filter that decomposes nonstationary data from non-linear sources.The NA-MEMD is used in this application to extract correctly the associated mode to the density wave oscillation, before to compute the SampEn.The proposed SampEn + NA-MEMD rolling window methodology allowed categorize stable BWR signals from unstable ones.The computed SampEn values of the stable signals oscillated around 1.2 (a high SampEN value is linked to a broad band noncoherent noise (BBNCN) associated with stable behavior).The SampEn values of unstable signalsFigure 31: Mean(SampEn) values for the 36 LPRMs of level 4 of the Ringhals BWR unit.fluctuateclose to 0.6 (they might be lower than that, but overall such SampEn estimates are below the threshold SampEn value of 0.8).A low value of SampEn points in this case to a simple cyclic function that is linked to density wave instability.100% of the typical BWR studied cases were classified with success into stable or unstable categories.100% of the Forsmakr Case 4 studied LPRMs confirm BWR instability through our proposed SampEn + NA-MEMD methodology.100% of the studied LPRMs of level 2 and 4 of the Ringhals Case 9 cycle 14 confirm instability.According to these experiments, the SampEn paired with the NA-MEMD is a suitable candidate to be used as a BWR stability indicator.
A studied time series  = [x 1 , x 2 , ..., x  ] of length  Θ:TheHeavisidefunction ‖ ⋅ ‖ 1 : The Chebyshev distance between V  and V +1: An average across all embedding vectors  +1  V  :A subvector of the studied time series, with  consecutive data points taken from  :