Enrico Zio
email: enrico.zio@polimi.it

Francesco Di Maio

A Data-Driven Fuzzy Approach for Predicting the Remaining Useful Life

Keywords: Prognostics, Residual Useful Life, Recovery Time, Emergency Accident Management, Operator Support System, Nuclear Power Plant, Lead-Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS), Pointwise Fuzzy Similarity Analysis

HAL is

Introduction

Existing Nuclear Power Plants (NPPs) strive to improve safety, maintain availability and reduce operation and maintenance costs. Moreover, power uprates and life extensions increase the needs of techniques for diagnosing and prognosing the NPPs health, because of the more probable occurrence of component degradation and failure as load is increased, or changed, and age advances.

In generality, prognostics is an important and challenging task in Reliability, Availability, Maintainability and Safety (RAMS). The primary goal of a prognostic system is to indicate whether the structure, system or component (SSC) of interest can perform its function throughout its lifetime with reasonable assurance and, in case it cannot, to estimate the Remaining Useful Life (RUL), i.e.

the lifetime remaining before it can no longer perform its function [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition based maintenance[END_REF].

Approaches to prognostics can be categorized broadly into model-based and data-driven [START_REF] Chiang | Fault detection and diagnosis in industrial systems[END_REF]. Model-based prognostics attempts to incorporate physical models of the system into the estimation of the RUL. However, uncertainty due to the assumptions and simplifications of the adopted models may pose significant limitations on this approach. On the contrary, data-driven techniques utilize monitored operational data related to system health. They can be beneficial when understanding of first principles of system operation is not straightforward or when the system is so complex that developing an accurate model is prohibitively expensive.

Data-driven approaches can often be deployed quickly and cheaply, and still provide wide coverage of system behavior. An added value of data-driven techniques is their ability to transform highdimensional noisy data into lower dimensional information useful for decision-making [START_REF] Dragomir | [END_REF].

In practice, when the cost/benefit ratio of using physics-based degradation evolution models is not favorable but sufficient data (possibly simulated) are available for constructing a map of the damage space, data-driven techniques may be more apt to the prognostics task. Furthermore, recent advances in sensor technology and refined simulation capabilities enable us to continuously monitor the health of operating components and manage the related large amount of reference data.

Data-driven approaches can be divided into two categories: statistical techniques (regression methods, ARMA models, etc.) and Artificial Intelligence (AI) techniques (neural networks, fuzzy systems, etc.).

The most natural data-driven technique for RUL estimation is to fit a curve of the available data of the component degradation evolution using regression models and extrapolate the curve to the criteria indicating failure. However, the available component history sometimes may be short and incomplete, so that extrapolation may lead to large errors. The same problem arises when employing ARMA models, although the method can handle also situation in which more run-tofailure data are unavailable or insufficient [START_REF] Yan | A Prognostic Algorithm for Machine performance Assessment and its Application[END_REF].

With respect to AI techniques, the most commonly used prediction methods are based on Neural Networks [Peel et al., 2008;[START_REF] Barlett | Nuclear Power Plant Status Diagnostics Using an Artificial Neural Network[END_REF][START_REF] Santosh | Diagnostic System for Identification of Accident Scenarios in Nuclear Power Plants using Artificial Neural Networks[END_REF]. For prognostic tasks, the most promising methods are Recurrent Neural Networks (RNNs) [START_REF] Campolucci | On-Line Learning Algorithms of Locally Recurrent Neural Networks[END_REF] and Neuro-Fuzzy (NF) systems [START_REF] Wang | Prognosis of Machine Health Condition using Neuro-Fuzzy Systems[END_REF]. In spite of the recognized power of neural network modeling techniques, limitations still exist which prevent their wide use in safety critical applications, e.g. in nuclear technology. Skepticism against their use in such application relates to the lack of a systematic approach for selecting the structure and parameters of the model and to their black-box character which limits intuition with respect to the understanding of their performance [START_REF] Wang | A similarity based prognostic approach for Remaining Useful life estimation of Engineered Systems[END_REF]].

An opportunity for increased transparency and openness of data-driven model is offered by fuzzy logic methods, which are increasingly proposed in modern diagnostic technologies. Based on the principles of Zadeh's fuzzy set theory, fuzzy logic provides a formal mathematical framework for dealing with the vagueness of everyday reasoning [Zadeh, 1965]. As opposed to binary reasoning based on ordinary set theory, within the fuzzy logic framework measurement uncertainty and estimation imprecision can be properly accommodated [START_REF] Yuan | Data driven identification of key variables[END_REF][START_REF] Zio | Identification of nuclear transients via optimized fuzzy clustering[END_REF].

The goal of this work is to develop a prognostics tool to be embedded in an operator support system for emergency accident management. In this context, the remaining useful life represents the recovery time available to the operator for deciding on the corrective action to avoid system failure.

In emergency situations, poor decisions may be taken because of the short time available for sorting out the relevant information and lack of expert knowledge [START_REF] Glasstone | Nuclear reactor engineering[END_REF]; on the contrary, timely and correct decisions can prevent an incident from developing into a severe accident or mitigate its undesired consequences.

The paper presents the computational framework devised for integrating on-line monitored data for system RUL prediction. A set of failure trajectory patterns (hereafter called reference trajectory patterns) is collected in a reference library and a fuzzy based data-driven similarity analysis is performed for predicting the remaining life of a newly developing failure trajectory (hereafter called test trajectory pattern). A novelty of the proposal lies in the reliance on a fuzzy definition of trajectory pattern similarity to capture and integrate the ambiguous information carried by the measured signal. More specifically, the pattern matching process is based on a fuzzy distance evaluation algorithm between the reference and test trajectory patterns [Angstenberger, 2001].

An application is presented with reference to dynamic failure scenarios of the Lead Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS) with digital Instrumentation and Control (I&C) [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF].

The paper is structured as follows. Section 2 provides a detailed description of the algorithm for the fuzzy data similarity comparison underpinning the prognostic approach to RUL estimation.

Section 3 presents the mechanistic model of the LBE-XADS. In Section 4, the results of the application of the approach to LBE-XADS accident scenarios are presented. Finally, some conclusions are drawn in Section 5.

Methodology

It is assumed that N trajectories of evolution in time of values of relevant signals (reference trajectory patterns) are available from measurements collected in dynamic failure scenarios of the system of interest. These trajectories last all the way to system failure, i.e., to the time when anyone of the signals reaches the threshold value beyond which the system loses its functionality.

A failure scenario is developing in the system, whose signals trajectory in time (test trajectory pattern) is monitored.

The idea underpinning the RUL estimation is to evaluate the similarity between the test trajectory pattern of the developing failure scenario and the N reference trajectory patterns stored in the database and use the RULs of these latter to estimate the RUL of the former, accounting for how similar they are [Angstenberger, 2001].

Figure 1 shows a schematic sketch of the computational framework; for ease of illustration, a single signal   ft is considered (a generalization to the case of multidimensional trajectories is straightforward).

Figure 1 The flowchart of the fuzzy based data-driven approach

In what follows, the steps of the procedure are reported in details:

-Step 1: fault detection. The signal   ft is continuously monitored throughout the time horizon of observation T , starting from (discrete) time 1 t  ; at each time t , its value is recorded and appended to the vector of the values collected at the time steps preceding t .

For reasons which will become clear in the following, the database constituting the reference trajectory pattern library is organized in a reference matrix

            || 11 ˆ| | | | i ff ii ii fi i t t i t t ff RUL t MTTF t t t RUL t i t t i t t        (1)
where i f t is the system failure time along the i-th trajectory (i.e., the time when the signal value exceeds the threshold beyond which the system loses its functionality),  

|| i f i t t
 is the cardinality of the set of reference trajectories whose failure time is larger than t and   i RUL t is their remaining useful life starting from t . At the following time steps, the algorithm continues to update the estimate of

    RUL t MTTF t 
, until a fault is detected upon a deviation of the signal from its nominal value, outside its range of allowed variability; at this time, the RUL estimation algorithm is put to work to match the similarity of the developing signal trajectory evolution to those in the reference library, and combine their failure times to provide an estimate of the RUL.

From the above explanation, it is clear that fault detection plays a relevant role in the RUL estimation procedure and structures it so that the RUL estimations are made up of two parts:

MTTF-based part, first and RUL matching-based part, afterwards. The dividing instance is the fault detection: when a fault is detected, from then on the test signal trajectory is compared with those of reference to estimate the RUL; before the detection, the RUL is estimated as the MTTF. In spite of the importance of fault detection, in this work this is not treated specifically and the focus is kept on the RUL estimation task after detection. Future works will address the task of fault detection within the proposed approach and with respect to two particular developments: i) the possibility of skipping the fault detection process and proceeding from the beginning with the matching of the evolving test trajectories with the reference trajectories stored in the database: it is foreseeable that before the fault, the test trajectory will match well with the non-faulty segments of the reference trajectories and the estimated RUL will approximate the MTTF; thus, apart from increased computational and storage memory burdens the results are expected to not degrade significantly with respect to those obtained in this work; ii) the possibility of including in the computation of the RUL the accuracy of the fault detection method in identifying a deviation of the test trajectory from the nominal pattern evolution at the early stage as well as account for the reliability of fault detection.

    1, 2,..., f t f t n t n t      is normalized in [0.2,0.8]. The pointwise difference   
between the n values of pattern   ft and those of the reference trajectory segment   , r i j is computed:

      ,, i j f t r i j   , 1, 2,..., iN  , 1, 2,..., jk  (2)
The matrix   Nk   contains the difference measures  

, ij  between all n-long segments of the reference trajectories and the test trajectory pattern of the monitored signal.

-Step 3: trajectory pointwise similarity and distance score computation. Classically, similarity is decided crisply: depending on whether the distance or similarity measure between two objects exceeds a specified threshold, the objects are classified into distinct categories of 'similar' or 'non-similar' [Angstenberger, 2001]. Such binary classification is restrictive when the situation is not so clear-cut and imprecision and uncertainty in similarity classification exist [START_REF] Zimmermann | Quantifying Vagueness in Decision Models[END_REF]. In practice, there are numerous cases in which the similarity measure should allow for a gradual transition between 'similar' and

r (i,1) r (i,15) r (i,5) … … …
'non-similar' [START_REF] Binaghi | Fuzzy Reasoning Approach to Similarity Evaluation in Image Analysis[END_REF][START_REF] Joentgen | Dynamic Fuzzy Data Analysis based on Similarity between Functions[END_REF]. This can be achieved by resorting to a fuzzy logic modeling paradigm in which the pointwise difference of two trajectories is judged for similarity with respect to an "approximately zero" fuzzy set (FS) specified by a function which maps the elements  

, ij

 of the difference matrix   Nk   into their values   , ij
 of membership to the condition of "approximately zero". The distance score   , dij between two trajectory segments is then computed as:

    , 1 , d i j i j   , 1, 2,..., iN  , 1, 2,..., jk  (3)
Common membership functions can be used for the definition of the FS, e.g. triangular, trapezoidal, and bell-shaped [START_REF] Dubois | Weighted Fuzzy Pattern Matching[END_REF]. In the exploration of these membership functions defined on a fixed support of the fuzzy set "approximately zero" and centered in the value zero, no particular sensitivity of the RUL estimation results has been recorded with respect to the particular shape of the function; in the application illustrated in the work, the following bell-shaped function has turned out to give more robust results due to its gradual smoothness:

      2 2 ln , , ij i j e           (4)
The arbitrary parameters  and  can be set by the analyst to shape the desired interpretation of similarity into the fuzzy set: the larger the value of the ratio  

  i RUL t :     ˆfi ii i t t RUL t w RUL t    , 1, 2,..., iN  (5)
The ideas behind the weighting of the individual   i RUL t is that: i) all failure trajectories in the reference library bring useful information for determining the RUL of the trajectory currently developing; ii) those segments of the reference trajectories which are most similar to the most recent segment of length n of the currently developing failure trajectory should be more informative in the extrapolation of the occurring trajectory to failure.

To assign the weight i w , the minimum distance * i d along the i th row of the matrix of Eq. ( 3) is first identified:

  * 1,..., min , i j k d d i j   , 1, 2,..., iN  (6)
The weight i w is then computed, resorting to the arbitrarily chosen decreasing monotone function:

  * 1 * 1 i d ii w d e         , 1, 2,..., iN  (7)
The value of β in ( 7) is arbitrarily taken equal to that used in Eq. ( 4), for reducing the burden of the parametric sensitivity analysis performed for identifying the parameters most affecting the results.

Then, the weight w i is normalized by:

1

N i i e e w w w    (8)
As shown in Figure 4, the smaller the minimum distance the larger the weight given to the ith trajectory. This allows a conservative RUL estimation, biased towards "pessimistic" predictions of the RUL, because in the case that more than one segment along the i-th reference trajectory is closest to the developing test trajectory, the latest one is taken, i.e. the one closest to failure.

Then, the estimate   RUL t of the remaining useful life along the developing trajectory is simply computed as in Eq. ( 5) with weights i w evaluated by Eq. ( 8).

n-long developing test trajectory

i f t j r (i,1) r (i,15) r (i,5) … … …   i RUL t
Latest in-life reference segment with minimum distance M j t j

The LBE-XADS

The Lead-Bismuth Eutectic eXperimental Accelerator Driven System (LBE-XADS) is a sub-critical, fast reactor in which the fission process for providing thermal power  

Pt is sustained by an external neutron source through spallation reaction by a proton beam  

Qt accelerated by a synchrotron on a lead-bismuth eutectic target [START_REF] Bowman | Nuclear energy generation and waste transmutation using an accelerator-driven intense thermal neutron source[END_REF][START_REF] Carminati | An Energy Amplifier for Cleaner and Inexhaustible Nuclear Energy Production Driven by a Particle Beam Accelerator[END_REF][START_REF] Rubbia | Conceptual Design of a Fast Neutron Operated High Power Energy Amplifier[END_REF]Van Tuyle et al., 1993;[START_REF] Venneri | Accelerator-driven Transmutation of Waste (ATW) -A New Method for Reducing the Long-term Radioactivity of Commercial Nuclear Waste[END_REF]. A simplified scheme of the plant is sketched in Figure 7. The primary cooling system is of pool-type with Lead-Bismuth Eutectic A dedicated, dynamic simulation model has been implemented in SIMULINK for providing a simplified, lumped and zero-dimensional description of the coupled neutronic and thermohydraulic evolution of the system [START_REF] Cammi | Modelling and control strategy of the Italian LBE-XADS[END_REF]. The model allows the simulation of the system controlled dynamics as well as of the free dynamics when the control module is deactivated and the air cooler flow is kept constant.

Both feedforward and feedback digital control schemes have been adopted for the operation

Results

Application of the RUL(t) estimation procedure

The application of the procedure for computing the estimate

 

RUL t (Section 2) is hereafter illustrated; the goal is to demonstrate its feasibility for fast and reliable data-driven prognostics which can serve as operator support in emergency accident management. A database of reference trajectories for N = 6400 failure scenarios (differing in faulty components, time of faults occurrence and faults magnitude) is organized in the reference matrix   Nk R  , where N = 6400 and

T k n  =30,
respectively. The elements   , r i j of the reference matrix will be compared for similarity with a test trajectory pattern containing the values of the temperature signal of the latest 100 time steps of the trajectory. For each of the test trajectories the procedural steps are performed as follows:

Step 1: fault detection.

The signal   ft monitored starting from time 1 t  [s] to the mission time T = 3000 [s],
with time step of 1 [s], is the LBE primary coolant temperature. This means that the actual safety parameter driving the RUL, i.e., the diathermic oil temperature, is followed only indirectly through the LBE temperature values, thus challenging further the robustness of the approach with respect to non-observable safety parameters. Obviously, other choices for   ft may also be valid to drive the estimation, e.g. the diathermic oil temperature directly, another indirect parameter like the air flow or any combination of the observable parameters. At each time step t , its value is appended and stored in the vector containing the 1 n  =99 values of the signal collected at the previous times. The Mean Time to Failure   MTTF t is calculated resorting to Eq. ( 1) and

 

RUL t is set equal to   MTTF t for each time step, until any component failure is detected; the fault detection activates the on-line fuzzy based data-driven RUL estimation.

Step 2: trajectory pointwise difference computation. , implying strong sharpness in the FS and thus in the similarity requirement (Figure 4). The distance scores

 

, dij are then computed by Eq. ( 3), 1, 2,...,6400 i  , 1, 2,...,30 j  .

Step 4: weight definition.

The minimum distances * i d are evaluated (Eq. ( 6)), and the relative normalized weights i w calculated through Eqs. ( 7) and ( 8), 1, 2,...,6400 i  .

Step 5: RUL i (t) and RUL(t) estimation.

For each reference trajectory in the library, an estimate   ˆi RUL t for the test trajectory is computed (Eq. ( 9), 1, 2,...,6400 i 

); then, the   i RUL t are aggregated in the weighted sum (Eq.

(

Performance evaluation of the RUL(t) estimation procedure

The performance of the RUL estimation procedure has been tested extensively on a batch of P=1280 test trajectories, different from the reference trajectories. Figure 18 Figure 19 shows the empirical probability density function of the mean relative error evaluated along 3000 [s] of evolution of the 1280 test trajectories. The distribution is highly skewed towards small values, with mean and median values equal to 0.09 and 0.05, respectively. This proves that the procedure most frequently makes small relative estimation errors. 4) and ( 7), and including the case in which the FS "approximately zero" is not introduced in the step 3 of the procedure (i.e., β=∞). It is seen that:

-the accuracy in the estimation of the RUL improves over time: as the RUL decreases, the performance index approximate its mean value 0.05.

-the performance degrades as the parameter  influencing the width of the fuzzy set "approximately zero" increases: the crisper the similarity definition, the worse the estimation performance.

For the application of the present work, the value of 0.01

Figure 3

 3 Figure 2 A reference signal trajectory on a time horizon T=3000 [s]

Figure 4

 4 Figure 4 Bell-shaped membership functions with different values of parameter β

Figure 5

 5 Figure 5 Weight vs. distance: the smaller the distance, the larger the weight value

  Figure 6 The remaining time before failure

(

  LBE) liquid metal coolant leaving the top of the core, at full power nominal conditions, at temperature , CP LB  equal to 400 °C and then re-entering the core from the bottom through the downcomer at temperature , PC LB  equal to 300 °C. The average in-core temperature of the LBE , secondary cooling system is a flow of an organic diathermic oil at 290-320 °C, at full power conditions. Cooling of the diathermic oil in each loop is obtained through an air flow   a t  provided by three air coolers connected in series.

Figure 7

 7 Figure 7 LBE-XADS simplified schematics. A = Accelerator; C = core; P = primary heat exchanger; S = secondary heat exchanger

Figure 12

 12 Figure 8 Transient trajectory belonging to the lowtemperature failure mode. System fails shortly

  zero" FS. The bell-shaped function of Eq. (4) is taken with parameters values 0

Figure 17

 17 Figure13RUL estimation for the trajectory of Figure8belonging to the low-temperature failure mode

Figure 19

 19 Figure 18 Predicted and actual remaining life for 25 test accident scenarios

Figure 20

 20 Figure 20 Relative Error evaluated each 50 [s] starting from 600 [s] before failure, for 1280 trajectories, with different pointwise similarity fuzziness

   

	Nk R  , where r i j is the j-th segment of length n of the values of the i-th reference T k n  , whose generic element   , trajectory, 1, 2,..., iN  , 1, 2,..., jk  (Figures 2-3), normalized in the range [0.2,0.8]. As long as no signal deviation is detected, the system is qualified as working in nominal conditions and the estimate   RUL t of the remaining useful life made at the generic time t is taken equal to the system Mean Time to Failure   MTTF t , obtained from the remaining
	useful life	 

i RUL t of all the failure trajectories in the reference library:

  MTTF t value; this results in a conservative trend of initial anticipation of the RUL associated to trajectories whose failure actually occurs late in life. The largest RUL estimation errors occur for those trajectories in which the component failure is of low-magnitude, whose effect only slowly drives the system to failure and the prediction away from

					shows the	  RUL t
	predictions (continuous line with dark bullets) compared to the actual remaining lives for 25 failure
	trajectories (graphically appended in sequence), after the fault has been detected (i.e., starting from
	step 2 of the algorithm). After fault detection, the initial predictions for each test trajectory tend to
	be similar, regardless of the eventual length of the test trajectory life duration because the initial
	deviation from the nominal evolution is little sensitive and only slightly moves the	  RUL t
	prediction away from the	 
					  MTTF t
	towards the true	  RUL t .
		To globally quantify the performance of the procedure, the mean relative error (RE) at time
	t, between the estimate	  RUL t and its true value	  RUL t , is introduced:
					  RE t	1  P p   P	  1 pp   p RUL t RUL t   RUL t 	(10)
	where	  RUL t is the actual remaining useful life at time t of test pattern p, and p	  RUL t its ˆp
	estimate,	1, 2,..., pP 	.


Enrico Zio
email: enrico.zio@polimi.it

Francesco Di Maio

Jiejuan Tong

SAFETY MARGINS CONFIDENCE ESTIMATION FOR A PASSIVE RESIDUAL HEAT REMOVAL SYSTEM

For licensing purposes, safety cases of Nuclear Power Plants (NPPs) must be presented at the Regulatory Authority with the necessary confidence on the outcomes of the models used to analyze the plant safety behavior. In the present work, we consider the problem of providing a quantitative indication of the confidence in the safety margin estimation by a model with uncertain inputs giving in output the maximum outlet water temperature of the Residual Heat Removal system (RHRs) in accident scenarios of the High Temperature Reactor-Pebble Modular (HTR-PM). The quantitative evaluation is carried out by means of a computational procedure of literature, based on Order Statistics (OS). The procedure is analyzed with respect to some of its key parameters defining the sample size and the number of uncertain inputs.

NOTATION AND LIST OF ACRONYMS

Function that maps the input vector x into the output vector y L j Lower threshold for the j th output parameter U j

Upper threshold for the j th output parameter For the Beyond Design Basis Accidents (BDBAs), the conservative approach is being challenged by a more realistic, Best-Estimate (BE) analysis, which sets forth the calculation of safety margins with realistic models and BE assumptions to account for the consequences related to the failures in some protective barriers.

On one side, the reduction in the conservatism of the analyses leads to more efficient plant design and operation. On the other side, the relaxation of the conservatisms entails that sensitivity and uncertainty analyses be carried out to properly quantify the safety margins while capturing the associated uncertainty for confidence evaluation. This requires a revision in probabilistic terms of the concept of safety margins [START_REF] Gavrilas | A Generalized Framework for Assessment of Safety Margins in Nuclear Power Plants[END_REF] and repeated model runs for the associated sensitivity and uncertainty analyses.

This is even more so when i) analyzing requests for changes to the licensing bases, within a risk-informed decision-making philosophy and ii) checking design solutions of new

NPPs whose safety analysis relies on newly developed models and codes, because the combination of the uncertainties in the analysis could reduce significantly and in an unexpected way the safety margins [START_REF] Usnrc | An Approach for using probabilistic risk assessment in riskinformed decisions on plant-specific changes to the licensing basis[END_REF][START_REF] Martorell | Safety Margins Estimation Method Considering Uncertainties within Risk-Informed Decision Making Framework[END_REF]. This situation may increase the risk of accidents, leading to the dissatisfaction of some of the Basic Safety Functions (BSFs) (i.e., reactivity control, residual heat removal, primary pressure control and containment release) that have to be carried out by the implemented safety systems to avoid major dangerous consequences.

In general, uncertainty affecting the plant behavior can be considered of two types: that due to inherent variability in the system behavior and that due to lack of knowledge and information on the system. The former type of uncertainty is often referred to as objective, aleatory, stochastic whereas the latter is called subjective, epistemic, state-ofknowledge [Apostolakis, 1990;Helton, 2004].

The distinction between aleatory and epistemic uncertainty plays a particularly important role in the risk assessment framework applied to complex engineered systems. In the context of risk analysis, the aleatory uncertainty is related to the occurrence of the events which define the various possible accident scenarios whereas epistemic uncertainty arises from a lack of knowledge of fixed but poorly known parameter values entering the evaluation of the probabilities and consequences of the accident scenarios.

The present work addresses the epistemic uncertainty affecting the evaluation of the safety margins. Under the probabilistic viewpoint here undertaken to represent uncertainties, the code for safety margin evaluation needs to be repeatedly run with different values of the thermal-hydraulic parameters, sampled from predefined probability distributions; the outcomes of these runs are then statistically analyzed to estimate with a specified confidence a given percentile of the distribution of the safety parameter used to calculate the safety margin [START_REF] Guba | Statistical Aspects of Best Estimate Method-1[END_REF][START_REF] Nutt | Evaluation of Nuclear Safety from the outputs of Computer Codes in the Presence of Uncertainties[END_REF]. The confidence intervals of the estimated safety margins is also computed: this additional information provides a realistic refinement of the estimates that is beneficial to power plant owners; on the other hand, from the viewpoint of the regulatory body it increases the robustness of the safety case by allowing verification of the fact that uncertainty in the safety margin estimates does not lead to exceedance of the regulatory safety thresholds.

In general terms, because of the large computing times required to run the codes, the procedure can be computationally quite expensive. Thus, tin this work the statistical analyses of the model evaluations for obtaining confidence intervals for safety parameters estimates rely on the use of Order Statistics (OS), along a non-parametric approach initially explored by [Wilks, 1941;Wilks, 1942]; this brings the advantage that the number of code calculations needed for safety margins evaluation is independent of the number of uncertain input parameters and provides a given confidence on the reliability of the calculated point-estimate obtained with a limited number of code runs.

Figure 1 shows a schematic sketch of the non-parametric procedure here adopted [Secchi et al., 2008]; for ease of illustration, a single safety parameter y is considered. By this procedure for safety margin calculation, the analyst can produce results with the level of confidence against uncertainty required for presenting a robust safety case to the Regulatory Authority.

The approach is applied to a case study regarding a set of accident scenarios related to the Residual Heat Removal system (RHRs) of the High Temperature Reactor-Pebble Modular (HTR-PM) [START_REF] Zhengy | Characteristics of the 250MW Pebble-Bed Modular High Temperature Gas-Cooled Reactor in Depressurized Loss of Coolant Accidents[END_REF]. The application is of safety significance because the RHRs is a passive safety system which operates without external input energy [IAEA, 1991] and is thus expected to contribute significantly to the improvement of plant safety [START_REF] Mathews | Functional reliability analysis of safety grade decay heat removal system of Indian 500MWe PFBR[END_REF].

However, the uncertainties involved in the modeling of the behavior of passive systems are usually larger than in active systems, due to lack of data on some underlying phenomena and scarce or null operating experience over the wide range of conditions encountered during operation [START_REF] Pagani | The impact of uncertainties on the performance of passive systems[END_REF].

This situation may actually increase the risk of accidents leading to the dissatisfaction of some of the BSFs (i.e., reactivity control, residual heat removal, primary pressure control and containment release) that the safety systems are designed for. In fact, deviations in the natural forces and in the conditions of the underlying physical principles from the model expected ones can impair the function of the system itself [START_REF] Marquès | Methodology for the reliability evaluation of a passive system and its integration into a probabilistic safety assessment[END_REF]Burgazzi, 2007].

The paper organization is as follows. In Section 2, the basic principles underpinning the BE nuclear safety analysis in the presence of uncertainties are provided together with an illustration of the method for uncertainty estimation here employed. In Section 3, the main characteristics of the High Temperature Reactor-Pebble Modular (HTR-PM) are briefly introduced, the Residual Heat Removal system (RHRs) accident scenarios are described and the simulations performed to analyze the system response to the accident scenario are presented. In Section 4, the results of the application of the proposed framework for the estimation of the safety margin of the maximum outlet water temperature reached during the accidents described in Section 3 are provided. Finally, some conclusions are drawn in Section 5.

Batch number g:

1 2 … G Code evaluations: ( ) 1 1 ŷ , ( ) 1 2 ŷ , …, ( ) 1 N ŷ ( ) 2 1 ŷ , ( ) 2 2 ŷ , …, ( ) 2 N ŷ … ( ) 1 G ŷ , ( ) 2 G ŷ , …, ( )
1 ŷ γ γ th percentile estimate: ( ) 2 ŷ γ … γ th percentile estimate: ( ) G ŷ γ Empirical γ th percentile distribution ( ) ( ) ( ) { } 1 2 G ˆˆˆŶ y , y ,..., y γ γ γ = Order Statistics 1) δ th percentile: ŷ γ δ
2) Confidence Interval:

( ) ( ) r s ˆŷ , y γ γ       Figure 1
Flowchart of the non-parametric procedure for percentile and confidence interval estimation [Secchi et al., 2008] 2 SAFETY MARGIN ESTIMATION

The mathematical model Code Safety Parameter Evaluation, yi Step 1: Code calculations

Step 2: Code batch-calculations

Step 3: OS batch-percentile estimation

Step 4: OS percentile estimation

Step 5: Confidence interval calculation

i x = (x1 x2 xp)i … Input parameters
A quantitative model for safety analysis of a nuclear power plant may be viewed as composed of three main elements: an input vector

{ } 1 2 , ,..., p x x x x =
, a BE simulator code and an output vector

{ } 1 2
, ,..., l y y y y =

. The elements of the input vector x are all the model parameters and input variables needed to calculate one realization of the output variables y describing the system response. The simulation code can be regarded as a black box which implements the complex, multidimensional, nonlinear function that maps the input vector x into the output vector y [START_REF] Guba | Statistical Aspects of Best Estimate Method-1[END_REF]]:

( )

y T x = (1)
For fixed values of x , the output values y is deterministically computed.

Formulation of safety margins in presence of uncertainty

The defense-in-depth principle of nuclear safety is founded on the implementation of protective barriers between the radioactive products and the environment. Each barrier is a physical device whose integrity is measured with reference to given characteristic safety parameters. When a barrier is subjected to abnormal conditions, some of the related safety parameters may exceed their safety envelope, which results in the failure of the barrier.

With reference to a generic accident scenario k and a characteristic safety variable ( ) k y j to be limited from above by an upper threshold limit j U , the safety margin M(y j ,k) can be defined as:

( )

, j j j M y k U y = - (2)
The dual definition of the safety margin for a safety parameter y j to be limited from below by a threshold value j L is straightforward.

Estimation of safety margins

The safety margins are calculated by running the thermal-hydraulic codes used for safety analysis. In presence of uncertainty, a large number of runs of the code may be required to adequately represent the full distribution of the safety parameter values in output. The problem is that the computer runs of the complex models of plant dynamics used for safety analysis are computationally very expensive. To overcome this hurdle one may resort either to simplified analytical/numerical models, such as those based on lumped effective parameters [START_REF] Marseguerra | Using Genetic Algorithms for Calibrating Simplified Models of Nuclear Reactor Dynamics[END_REF], or to empirical models, e.g. artificial neural networks and fuzzy logic systems [START_REF] Marseguerra | Variance Decomposition-Based Sensitivity Analysis via Neural Networks[END_REF], suitably set up so as to best fit to the data available from the plant.

Another possibility to reduce the computational burden, which may be used in combination with the former ones, is to only compute some percentiles of the output distribution, estimated with a limited number of runs. In this case, the confidence in the estimates becomes crucial for decision making and must, thus, be quantified [Wilks, 1941;Wilks, 1942;[START_REF] Guba | Statistical Aspects of Best Estimate Method-1[END_REF][START_REF] Nutt | Evaluation of Nuclear Safety from the outputs of Computer Codes in the Presence of Uncertainties[END_REF].

This latter approach is undertaken in this work and the problem of confidence building and quantification is addressed. A sample of a small number N of input parameter values is drawn by the Monte Carlo method from their probability distributions. The sample of N input vectors thereby obtained is input to the code which is correspondingly run N times, thus producing a random sample of N output vectors. These can be used to estimate a given percentile of the safety margin probability distribution. To obtain the desired confidence in the safety margin percentile, the number N of code runs is defined on the basis of the Order Statistics (OS) methodology, in its nonparametric formulation which applies independently from the type of probability distribution of the output data under study (in this case unknown) [Wilks, 1941;Wilks, 1942]. As we shall see, this amounts to ordering the elements of the random sample by increasing value, the element in the r th place being the statistic of order r, and using the order statistics for estimating the percentiles of the distribution (Section 2.4) with the desired confidence. Following this methodology, the number of runs required can be kept low because only statistical intervals are estimated and not the full probability distributions of the output.

Estimation of percentiles using Order Statistics

For ease of illustration, let us refer the discussion to a one dimensional output y, e.g. the Pellet Cladding Temperature (PCT) or the cooling water outlet temperature.

The N runs of the code, each one with a different input vector i x , produce N output , it could be possible to give a sufficiently accurate estimate of the full distribution of the output y and draw probabilistic conclusions on where it lies with respect to the threshold value U (or L)

vectors i y , i=1,2,…,N. Let { } 1 2 N Y y ,
defining the safety limit for the integrity of the protective barriers. This would provide a more realistic assessment than the verification that a single run conservative estimate of y is within the safety envelope (e.g., PCT is less than 1200°C or the cooling water outlet temperature is less than 95°C).

Given the computational costs associated with the estimation of the full distribution, one is forced to focus on verifying that with some level of confidence β, a certain percentage γ of the calculated values of y that would be obtained from running the code falls within the safety envelope. The thresholds defining the safety envelope and the values of the confidence and percentage are set by the regulations, considering the risk associated with exceeding the specified range.

With reference to the safety parameter y to be limited from above by U, the approach aims at showing that the m th member y m of the N sorted outputs has a certain probability β of exceeding the unknown true γ th percentile y γ [Wald, 1943;[START_REF] Nutt | Evaluation of Nuclear Safety from the outputs of Computer Codes in the Presence of Uncertainties[END_REF]]. Then, one has a level of confidence β that the actual value of y γ is less than the value obtained for y m : if y m meets the criterion of being less than the safety threshold U, then the unknown y γ will do so, too. In other words, the probabilities γ and β are defined as

{ } P y y γ γ = < and { } P y y γ γ β = < .
Once γ and β are fixed, the OS method for calculating the ( ) β γ -percentile estimate follows the lines of [START_REF] Guba | Statistical Aspects of Best Estimate Method-1[END_REF][START_REF] Nutt | Evaluation of Nuclear Safety from the outputs of Computer Codes in the Presence of Uncertainties[END_REF] and consists in: Note that higher values of m in step (i) imply higher values of the sample size N but generate less conservative estimates of the γ-percentile; in any case, the sample size N,

i.e. the number of BE code runs, can be kept low because only intervals related to the γ th percentile are estimated and not the full probability distribution generating the data.

For the application of interest in the present work, the confidence intervals for the γ th percentile are computed using the γ th percentile estimates obtained from G batches of newly simulated accident transients. To verify that the estimated confidence interval does not exceed the safety thresholds, two-sided confidence intervals are computed.

The computational approach for percentile point-and confidence intervalsestimation

The computational approach of Figure 1 for the evaluation of a point estimate of the γpercentile and of a confidence interval associated to it, to be used for the safety margin evaluation, is taken from [Secchi et al., 2008]. In what follows, the steps of the procedure are repeated for reader's convenience (reference distributions and notations are shown in be the values of the order statistics;

ii) set r and s to positive integers satisfying the inequality 0 < r < (N+1)/2 < s ≤ N;

iii) then, the random interval

[ ] [ ] ˆ, r s y y γ γ   
 covers the median of the distribution of the estimator ŷ γ with probability ( ) ( )

1/ 2, 1, 1/ 2, 1, I N s s I N r r α = -+ - -+ (3)
where ( ) I c, j,k is the Regularized Incomplete Beta Function for non-singular cases [Kendall et al., 1979;Pál et al., 2002].

Hence, by fixing α we may find suitable r and s, for instance in a symmetric position with respect to (N+1)/2, such that At a first glance, the HTR-PM design has the following key technical features [Zhengy et al., 2008]:

• Characteristic coated fuel particles are used, which consist of uranium dioxide (UO 2 ) fuel kernel coated by tri-isotropic (TRISO) ceramics such as pyrolytic y = safety parameter carbon and silicon carbide (SiC), in order to retain fission products in the particle under a fuel cladding temperature of 1600ºC in accident cases.

• A one-zone core design is implemented, consisting of approximately 420,000 spherical fuel elements in a pebble-bed with a diameter of 3m and an average height of 11m.

• Ceramic materials of graphite and carbon bricks, which are high-temperature resistant, surround the active reactor core.

• Decay heat in the fuel elements is dissipated by means of heat conduction and radiation to the outside of the reactor pressure vessel, and then taken away to the ultimate heat sink by water cooling panels on the surface of the primary concrete cell. Therefore, no coolant flow through the reactor core is necessary for decay heat removal in case of loss of coolant flow or loss of pressure accidents.

Maximum accident fuel temperature shall be limited to 1600ºC.

• Spherical fuel elements are charged and discharged in a so-called "multi-pass" mode, which means that before the fuel elements reach the discharge burn-up, they go through the reactor core several times.

• Several of HTR-PM modular reactors can be built at one site to satisfy the power capacity demand of the utility. Some auxiliary systems and facilities can be shared among the modules.

The Passive RHRs

The enhanced safety of the HTR-PM is mainly due to the implementation of passive safety systems [START_REF] Zhao | Assessing Physical Process Failure Probability of Passive System in an Advanced Reactor by Using M.C. Method[END_REF].

Figure 3 sketches the equipment layout of one of the 3 loops of the RHR system implemented in the HTR-PM. The water cooled wall gets the heat from the reactor vessel by thermal radiation; then, the pipe transfers the water to the air-cooled heat exchanger located in the air-cooled tower; the cool air takes the heat away from the heat exchanger to the environment.

The RHRs blockage scenarios

The outlet water temperature , The accident scenarios considered in the present study are:

-Scenario A: 2/3 RHR loops out of service (i.e., simultaneously failed and/or under maintenance: this constitutes a BDBA)

-Scenario B: 1/3 RHR loops out of service (i.e., failed or under maintenance: this is a scenario included in the set of DBA)

-Scenario C: 0/3 RHR loops out of service (i.e., nominal condition)

Scenarios A and B lead to a temporary decrease in the cooling capability of the RHR and to a corresponding increase of , w out T which may exceed c T [START_REF] Zhao | Assessing Physical Process Failure Probability of Passive System in an Advanced Reactor by Using M.C. Method[END_REF]. For this reason, these scenarios are considered safety relevant and a careful analysis of it must be performed.

Simulations of the accident scenarios

A simplified zero-dimensional description of the thermo-hydraulic behavior of the RHRs has been implemented in MATLAB and used to simulate accident blockage transients.

The model allows the computation of the maximum outlet water temperature reached during an accident scenario.

The simulation code models the following phases of the process:

1. The residual heat radiates from the reactor vessel and other thermal sources to the water in the water-cooled wall;

2. Because of the difference in temperature, natural convection will initiate through water, in the water-cooled wall and pipes connected with the air-cooled heat exchanger; then, heat will transfer to the water side of the heat exchanger;

3. The heat will transfer by thermal conduction from the water side to the air side of the heat exchanger, due to the difference of temperature;

4. As the air-cooled heat exchanger is located in the air-cooled tower, natural convection of air will set up and take heat to the final heat trap-atmosphere.

The model fed with the nominal "best estimate" values of the input parameters is assumed to be "sufficiently best estimate". The RHRs accident complete blockage transients are generated by sampling the involved 37 input parameters from probability distributions defined on the basis of previous experience and/or information obtained by skilled experts (Table 1).

RESULTS

The non-parametric procedure for percentile estimation introduced in Section 2.5 is hereafter illustrated with reference to the estimation of the safety margin of the maximum outlet water temperature , w out T reached during the accident scenarios A, B and C of complete/partial blockage of the passive RHRs of the HTR-PM described in Section 3.

Order Statistics has been applied to a sample of maximum outlet water temperature values obtained by simulation, for estimating the γ th percentile with γ=0.95. Then, G estimates of the γ th percentile have been collected by the OS applied to G different batch samples. Finally, confidence intervals for the real γ th percentile have been evaluated.

The procedural steps described in Section 2.5 have been performed as follows:

Step 1: BE code calculations.

To demonstrate the feasibility of the procedure, we take m=1 and β=γ=0.95; this leads to the smallest sample size N=59 for the OS ( ) β γ -percentile estimates.

Step 2: Code batch-calculations.

A number of G=50 batches of N output values have been computed.

Step 3: OS batch-percentile estimation.

For each of the G=50 batches, the ( ) β γ -percentile estimates have been computed and collected in the sample Step 4: OS percentile estimation.

The median of the sample Ŷ and its safety margin are evaluated for the accident scenarios A, B and C. The results are provided in the second column of Tables 2, 3 and 4, respectively.

Step 5: Confidence interval calculation.

Because of the limitation on the sample size used in the estimation, the safety acceptance criteria cannot be based solely on the best estimate results. Hence, the uncertainty of the estimated safety margin must be properly informed, e.g. by computing its confidence interval. With reference to the accident scenarios A, B and C, the confidence interval of level α=0.95, with r=1 and s=49, are provided in the second column of Tables 2, 3 and4 It can be seen that:

-scenario A has a negative safety margin: the occurrence of this accident scenario has to be avoided;

-from scenario C to scenario B the safety margin shrinks;

-scenarios B and C can be classified as safe, because with positive safety margins.

This means that, in absence of any other component failure in the system, although one of the loops is failed the system can continue to produce energy; on the other hand, in case of occurrence of scenario B, a maintenance action has to be adopted quickly to reactivate the failed RHR loop; -as shown in Figure 4, in general terms, the probability distributions of the safety

Improvements of the estimation accuracy by means of higher values of m, N and G

It is known that using higher values of m, N and G would allow increasing the reliability on the estimated confidence interval [START_REF] Zio | Using Order Statistics and the Bootstrap Technique for quantifying Thermal-Hydraulic Code Uncertainties in the Estimation of Probabilistic Safety Margins[END_REF]. Indeed, conservatism is reduced by taking higher values of m, N and G [START_REF] Nutt | Evaluation of Nuclear Safety from the outputs of Computer Codes in the Presence of Uncertainties[END_REF] and, by comparison of the results, the analyst would feel reassured that the estimates obtained have a low probability of differing significantly from the true values (usually unknown), and that the estimated maximum outlet water temperature value satisfies the safety threshold limit U.

In this view, the same procedure detailed in Section 2.5 has been repeated increasing the number of samples G. For the accident scenarios A, B and C, the results are provided T [Yu et al, 2010a;Yu et al., 2010b]. The results are provided in Tables 8910. By comparison to Tables 234, respectively, it can be seen that an overall qualitative agreement exists between the safety margins evaluated resorting to the complete analysis considering all 37 input parameters and those obtained with only the 3 most relevant input parameters: this demonstrates the efficacy of sensitivity analysis and allows for an even faster safety margin evaluation, freeing the procedure from the numerous time-consuming samplings of all the input parameters values and related code calculations.

CONCLUSIONS

A computational framework of literature has been applied for the estimation of the safety margin on the maximum outlet water temperature of the passive RHRs reached during some accident scenarios of the HTR-PM.

The procedure exploits non-parametric Order Statistics performed on a limited number of BE code calculations for providing confidence intervals on the estimated percentiles. An analysis has been performed on the effects of some key parameters, related to the size of the statistical sample and on the number of uncertain input variables considered in the analysis.

The procedure has been demonstrated to give reliable (the estimates are similar despite of the increase of the number of simulation), robust (confidence intervals are very narrow)

and conservative (increasing the number of simulations, the estimates tend to narrow down to the true value) estimates of the 95 th percentiles of the safety parameters distributions.

The method has been demonstrated effective in that it is capable of indicating the passive system safety conditions, accounting for the uncertainties in the model parameters and in the estimate itself.

  element of the representative sample of independent input vectors ( ) i y i th element of the representative sample of independent output vectors k Index of the accident scenario M(y j ,k) Safety margin relative to the j th safety parameter for the k th accident scenario safety parameter relative to the k th accident scenario N Number of simulations r Position of the ordered sample of simulations Y Ordered set of values resulting from running the code N times β Confidence value γ Coverage value m Number of values that lie beyond the γ coverage extent Lower bound of random interval covering the median of the distribution of ŷ Upper bound of random interval covering the median of the distribution of ŷ 1 INTRODUCTION Conservative calculations are traditionally performed for the verification of the safety performance of a Nuclear Power Plant (NPP) under Design Basis Accidents (DBAs) conditions, in terms of the values reached by selected safety parameters in comparison to threshold values enforced by regulation to ensure that sufficient safety margins are maintained for the integrity of the defense-in-depth barriers. The differences between the conservatively computed safety parameter values and the thresholds give the so called safety margins. Conservatism is introduced in the calculations to account for the uncertainties in the model representation of the actual plant behavior.

  y ,..., y = be the ordered set of values resulting from running the code N times for N different input vectors

∑

  size N by fixing a positive integer m. The probability that at least m observations within a random sample of size N are greater than the γ-percentile of the distribution generating the sample Set p=β and compute the sample size N by solving the previous equation in terms of N.ii)Sorting the observations in the sample by increasing value, the element in the r th place being the statistic of order r.iii) Estimating the γ-percentile by setting ŷ γ equal to the statistic of order N-m+1, i.e. the m th largest observation in the sample; then

Figure--

  Figure2):

  as point estimate of the γpercentile of the output distribution for safety margin evaluations.-Step 5: Confidence interval calculation. In alternative to the point estimate ofStep 4, we generate a confidence interval estimate of the median of the distribution of the estimator ŷ

  median of the distribution of ŷ γ . Note that, given the fact that ŷ γ is the (β|γ)-estimator of the γ-percentile of the output distribution, we expect both the point estimate of Step 4 and the interval estimate of Step 5 to cover values larger than the true value of y γ (see Figure2for an illustration of the analysis setting). The method is efficient if these estimates will not be too conservative while guaranteeing the required level of confidence.

Figure 2

 2 Figure 2 Sketch of the coverage value γ, the confidence β and the confidence interval level α; representation of the (unknown) safety parameter probability distribution ( ) f y and its γ th percentile

Figure 3

 3 Figure 3 Schematics of 1 loop of the RHRs in the HTR-PM [Zhao et al., 2008]

  is considered as the safety parameter with respect to which the success or failure of the system is defined. From the engineering experience, when , w out T exceeds the critical temperature c T , local boiling may occur which can significantly worsen thermal transmission. Although engineering experience recommends c T =95°C, in this work, a value of c T =90°C was conservatively chosen.

  , with X=A, B and C, are different. As a consequence, also the distributions of the 95-th percentiles () results obtained with the G estimates of 0.95 ŷ , are different: the more the accident scenario is unsafe, the more the related distributions shift towards the safety threshold (Figure4).

Figure 4

 4 Figure 4 Conceptual sketch of the results obtained for the accident scenarios A, B and C

Figure 5

 5 Figure 5 Conceptual comparison of the results obtained for the accident scenario C with m=1, 50 and G=50, 150, respectively

Table 1 Parameters which are regarded relevant for the behavior of the passive RHRs.

 1

	N	Parameter	Distribution	Note
	1	W	Normal	Residual heat power
	2	Ta,in	Bi-Normal	Temperature of inlet air in the air cooled tower
	3	xi1	Uniform	Resistance coefficient of elbow
	4	xi2	Uniform	Resistance coefficient of header channel
	5	xiw	Uniform	Resistance coefficient of the water tank walls
	6	xia,in	Uniform	Sum of the resistance coefficients of inlet shutter and air cooling tower and silk net
	7	xia,out	Uniform	Sum of the resistance coefficients of outlet shutter and air cooling tower and silk net
	8	xia,narrow	Uniform	Resistance coefficient of the narrowest part of the tower
	9	Pa,in	Uniform	Pressure of the inlet air in the cooler tower
	10	dx	Uniform	Roughness of pipes
	11	Ha	Normal	Height of chimney
	12	La	Normal	Length of pipes in the exchanger
	13	Na	Normal	Total number of pipes in the air cooler
	14	Af	Normal	Air flow crossing are in the narrowest part of the tower
	15	Af,in	Normal	Inlet air flow crossing area in the tower
	16	Af,out	Normal	Outlet air flow crossing area from the tower
	17	Af,narrow	Normal	Crossing area in the narrowest part of the tower
	18	S1	Normal	Distance between centers of adjacent pipes in horizontal direction
	19	S2	Normal	Distance between centers of adjacent pipes in vertical direction
	20	S	Normal	Distance between fins in the ribbed pipe
	21	Da	Normal	Pipes inner diameter in the air cooling exchanger
	22	Do	Normal	Pipes outer diameter

Table 2 Median of the 95 th percentile distribution and the corresponding confidence interval, obtained for the accident scenario A, m=1, N=50 and G=50, 100.

 2

									Scenario B
		m									1
		N									59
	G Median safety margin	γ	( ) ŷ B 0.5		50 37.78	100 37.29
	Confidence interval	 	0.95	( ) ( ) 3 ˆ, y B	0.95	y	( ) ( ) 97 B	 	[40.72, 34.01]	[40.88, 33.19]

Table 3 Median of the 95 th percentile distribution and the corresponding confidence interval, obtained for the accident scenario B, m=1, N=50 and G=50, 100.

 3

								Scenario C
		m								1
		N								59
	G Median safety margin	( ) ŷ C 0.95 0.5		50 56.49	100 56.55
	Confidence interval	 	0.95	( ) ( ) 3 ˆ, y C	0.95	y	( ) ( ) 97 C	 	[58.92, 52.59]	[58.79, 52.70]

Table 4 Median of the 95 th percentile distribution and the corresponding confidence interval, obtained for the accident scenario C, m=1, N=50 and G=50, 100.

 4

ACKOWLEDGEMENTS

This work was supported in part by the Science & Technology Fellowship Programme in China financed by the European Commission (EuropeAid/127024/L/ACT/CN) and by the Important National Science & Technology Specific Projects in China (No. 2008ZX06902-009).

with m=1, N=59, G=100, α=0.95, r=3 and s=97 and reported in the third column of Tables 2, 3 and 4, respectively:

-despite the higher computational time -for all 3 accident scenarios, the application of the procedure with G=100 reduces the conservatism in the results since it provides the lowest point-estimate values of the maximum outlet water temperature (larger positive safety margins for scenarios B and C and smaller negative safety margin for Scenario A).

For comparison, the procedure detailed in Section 2.5 has been repeated increasing the number of simulations N: the value m=50 leads to a sample size N=1228. The results are provided in Tables 5,6 and 7. It can be seen that:

-increasing m and correspondingly increasing the number of values that are requested at least to lie beyond the "extent" γ of the cumulative probability, the estimation 0.95 ŷ of the percentile tends to narrow the true 0.95 y , which have been evaluated running the code 100000 times for each accident scenario. By testing the 100000 outputs for Normality by means of the Lilliefors Test [Lilliefors, 1967] and then applying a parametric approach for the 95-th percentile estimations, the estimate of 0.95 y has turned out to be equal to -3.22, 40.66 and 58.70 for scenarios A, B and C, respectively.

-for all the 3 accident scenarios the application of the procedure with m=50 increases the confidence on the estimated percentile value with respect to that with m=1, as demonstrated by the shrinking of the confidence intervals.

As an example, Figure 5 shows the combined effect of larger values of m and G with reference to the results provided for the accident scenario C. Finally, the physical conclusion that can be drawn from the analysis is that, for the safe operation of the plant, two RHR loops are enough; the 3 rd loop can be considered as a redundancy in the RHRs to guarantee high availability of the safety function. The large safety margins computed in case of the safe scenarios (B and C) suggest a possible improvement of the whole system design to avoid excessive conservatism leading to a more efficient plant design and operation.