Corresponding author: Francesco Di Maio, Politecnico di Milano, Energy Department, Via La Masa 34, 20156, Milano, Italy. Tel: +39 02 23996372, Fax: +39 02 2399 8566. Email: francesco.dimaio@polimi.it A REGIONAL SENSITIVITY ANALYSIS-BASED EXPERT SYSTEM FOR SAFETY MARGINS CONTROL Francesco Di Maio1, Alessandro Bandini1, Michele Damato1, Enrico Zio1,2 1Energy Department, Politecnico di Milano, Milan, Italy 2Chair on System Science and Energetic Challenge, Fondation – EDF, Centrale Supélec, Paris, France Abstract Thermal-Hydraulic (TH) codes are used to simulate the response of nuclear safety systems under transient and accident conditions. The outcomes of the simulations are used to verify the safety margins required for safe operation and make decisions on how to maintain them. In this work, a novel Expert System (ES) based on Regional Sensitivity Analysis (RSA) is developed to guide a system undergoing an accident scenario towards the safest conditions in the optimal number of operation. The ES proceeds by firstly identifying the (uncertain) system controllable variables (i.e., control rods position, feed- water flow rate, void fraction inside the steam generator, etc.) that most affect the system response by RSA; then, the limit-state function is calibrated on a dataset of outcomes of TH code runs and the system failure boundary (i.e., the limit surface) is defined on the set of (uncertain) TH input variables. Application of the ES is firstly shown with respect to an analytical case study that artificially simulates the response of a NPP to an accident scenario and, then, to a practical case study concerning the response of the pressurizer of a Pressurized Water Reactor (PWR). Keywords: Nuclear System, Expert System, Best-Estimate Thermal-Hydraulic Code, Safety Margins, Limit-State Function, Regional Sensitivity Analysis, Pressurizer. 1. INTRODUCTION Safety remains a priority for Nuclear Power Plants (NPPs) design and operation, as the release of radioactive material can result in catastrophic consequences in terms of casualties, environmental pollution and financial losses [Hsieh et al., 2012]. For the safety of NPPs, accident-preventive design and operation, and effective consequence mitigation plans are developed [Ma et al., 2011]. In practice, Emergency Operating Procedures (EOPs) are defined to provide the technical basis for suitable operator response to Design Basis Accidents (DBAs) and, nowadays also, Beyond Design Basis Accidents (BDBAs) [IAEA, 2006]. Fault Detection and Diagnosis (FDD) methods have been developed for detecting different types of faults and supervising the plant behavior for accident prevention [Ma et al., 2011; Park et al., 2015; IAEA, 2000; IAEA, 2004; IAEA, 2005]. Upon localization and isolation, under certain abnormal conditions, manual actions are conducted by the plant operators for restoring mailto:francesco.dimaio@polimi.it 2 normal operating conditions [Hsieh et al., 2012; Park et al., 2015]. For this, Abnormal Operating Procedures (AOPs) are provided, where the sequence of actions to undertake are given for different Initiating Events (IEs). Two practical issues are: i. it is difficult to ensure sufficient coverage of all possible IEs through AOPs, as they are developed based on historical data and there may exist significant IEs that have not yet been experienced [Park et al., 2015]; ii. human errors may occur during abnormal situations and incorrect AOPs may, then, be followed [Hsieh et al., 2012]. In abnormal situations, time is critical, and the amount of information and data to be examined is large. Decision Support Systems (DSSs) can aid the operators and reduce the possibility of errors [Ma et al., 2011; Hsieh et al., 2012; Park et al., 2015]. Examples of DSSs are: the Alarm and Diagnosis-Integrated Operator Support (ADIOS) system that attempts to avoid that too many alarms influence the operators judgement in a wrong way [Kim et al., 2001]; the Hidden Markov Model (HMM) for recognizing accidents in NPPs, proposed in [Kwon et al., 1999]; the Fault Diagnosis Advisory System (FDAS), based on dynamic neural networks [Lee et al., 2007]; the Dynamic System Doctor (DSD), a system-independent interactive software for on- line state/parameter estimation in dynamic system [Aldemir et al., 2001]; the Analysis of Dynamic Accident Progression Trees (ADAPT) methodology [Hakobyan et al, 2008]; the unsupervised clustering technique for NPP components fault diagnosis, based on Haar Wavelet Transform (HWT) and Fuzzy C-Means (FCM) algorithm, in [Baraldi et al., 2013]; the hybrid approach for balancing false and missed alarms, based on Correlation Analysis (CA), Genetic Algorithms (GAs) and Sequential Probability Ratio Tests (SPRTs) in [Di Maio et al., 2013]; the non-parametric decision strategy in [Al-Dahidi et al., 2014] to detect whether NPP components are in abnormal conditions, using Prediction Intervals (PIs) and Auto-Associative Kernel Regression (AAKR) [Hines et al., 1998]. However, there is no guarantee of improvement in the operators’ performance and sometimes the result could be the increase of operator workload, with negative implications on performance [Yoshikawa 2005; Kim et al., 2007; Hsieh et al., 2012]. In this paper, we propose an Expert System (ES) [McBride et al., 1993; Liao, 2005; Ikram et al., 2015] for operator aid, based on the system response outcomes obtained from a Thermal- Hydraulic (TH) code and Regional Sensitivity Analysis [Wei et al., 2014]. The TH code 3 reproduces the system physical behavior according to a mathematical model m that receives an input vector n x  (comprised of n input variables) and generates the system response vector )(xmy  , containing (at least) one safety parameter y . The input variables set n  can be partitioned into two subsets: one with controllable variables ( q  ), i.e. the levers under control of the plant operator, which can be manipulated to increase plant safety (e.g., reactor control rods position, rate of feed-water flow through the plant primary loops, accumulator water temperature and pressure, repair times, etc.), and the other one with non-controllable variables ( qn  ) [Di Maio et al., 2016]. In practice, the inputs x are uncertain [Apostolakis, 1990; Helton, 1996; Oberkampf et al., 2004]. With reference to scenario FE , for which the response of interest Y must be lower than a threshold y γ (imposed by regulation), the limit-state function G is defined as [Bourinet et al., 2011]: y XYXG  )()( (1) The input space can be split into a failure domain }0)(:{  XGXF and a safe domain }0)(:{  XGXS . The system failure boundary (limit surface) }0)(:{  XGXF , which separates S from F , can be projected onto an input controllable subset ( q  ) [Di Maio et al., 2016], which contains the controllable variables q XXX ,,, 21  that can be varied to increase the system safety margin )( XYγy  [Zio et al., 2010]. While the system is operating in its nominal state ex (usually assumed equal to the expected values ][ XE of the input variables X ), the safety parameter )( ee xmy  falls well below the safety threshold y γ , providing a margin of safety with respect to uncertainties in the model [Zio et al., 2008]. In abnormal conditions, the safety margin may decrease and this needs to be kept under control in order to avoid diverging to catastrophic accidents. An additional difficulty in the problem is due to the fact that any generic controllable input variable j X , qj ,,2,1  , can be affected by epistemic and/or aleatory uncertainty. In this paper, this uncertainty is modeled by considering that when j X is set equal to jx̂ , the input variable actually ranges between jj xx ˆ and jj xx ˆ according to the modified Probability Density Function (PDF): 4 )ˆ()ˆ( )( )( * jjXjjX jX jX xxFxxF xt xf jj j j   , ]ˆ,ˆ[ jjjjj xxxxx  (2) where )( jX xt j is the truncated PDF of jX over range ]ˆ,ˆ[ jjjj xxxx  :             ),ˆ(,0 ˆ,ˆ),( )ˆ,(,0 )( jjj jjjjjjX jjj jX xxx xxxxxxf xxx xt jj (3) Coherently, the system (uncertain) initial state is: ][][][ 0,0,20,220,210,110,10 qqqq x,xxxx,xxxx,xxxX   (4) In this paper, we present the development of an Expert System (ES) for maximizing the safety margin by exploiting the results of a Regional Sensitivity Analysis (RSA) to guide the search for the range (or state of the system X ) closest to 0 X , where ])([ XG is the smallest. In detail, the ES proposed in this paper is based on the Revised Ratio Functions (RRFs) [Wei et al., 2014]) that measure the impact on the mean (Revised Mean Ratio Function, HM ) and variance (Revised Variance Ratio Function, HV ) of the model output distribution due to the reduction in the range of variability of an individual input. For our purposes, the effect on the variance is of particular interest and it is here exploited to find which input controllable variable allows obtaining the largest decrease in G (i.e., the largest increase in the expected safety margin) when the system is in its current abnormal state X . Assuming that the system is operating at X and that range ]ˆˆ[ jjjj xx,xx  , ,,,2,1 qj  is reduced to jx̂ , the expected value (  ) and variance (Var ) of G are, respectively:               11 11 ,1 * 1 ))(()ˆ|,,(ˆ| xx xx xx xx xx xx n j Xjqj qq qq dxxfxxxgxG       (5) 5                 11 11 ,1 *2 1 ))((ˆ|)ˆ|,,(ˆ| xx xx xx xx xx xx n j Xjjqj qq qq dxxfxGxxxgxGVar      (6) The original RSA technique defines j HV of input j X , qj ,,2,1  , as the ratio between the residual variance ]ˆ|[ jxGVar with respect to the residual mean ]ˆ|[ jxG , and the full variance  GVar [Wei et al., 2014]:  GVar xGVar HV j j ]ˆ|[  (7) where:             11 11 1 *2 1 ))(()][),,((][ xx xx xx xx xx xx n Xq qq qq dxxfGxxgGVar       (8)               11 11 1 * 1 ))((),,( xx xx xx xx xx xx n Xq qq qq dxxfxxgG       (9) As Eq. (7) shows, j HV indicates the actual reduction of the model output variance due to the restriction of the range of j X . The larger j HV is, the more function G varies (i.e., decreases or increases) for a perturbation of the other inputs X ( j ). Therefore, we define the revised * j HV as the ratio between the variance of G , computed over the range ]ˆˆ[ jjjj xx,xx  , qj ,,2,1  , and the reduced range ( x̂ ) of each of the other inputs ,X q,,2,1  , ,j and the full variance  GVar . It is clear that the larger *jHV is, the more function G varies and, thus, a larger increase in the system safety margins can be achieved through a perturbation of the j-th input. Through the approach described, we aim at providing the plant operator with an effective and unequivocal AOP to keep the safety margins and avoid system failure, giving clear indications of which controllable variables to modify and by how much, when an IE has occurred. The rest of the paper is organized as follows. Section 2 illustrates the proposed approach. Section 3 shows an application of the approach to an analytical case study that artificially 6 simulates the response of a NPP to an accident scenario. Section 4 contains the application of the approach to the pressurizer of a Pressurized Water Reactor (PWR). Conclusions are given in Section 5. 2. AN EXPERT SYSTEM FOR THE MAXIMIZATION OF SAFETY MARGINS DURING SYSTEM OPERATION For the sake of clarity, but without loss of generality, in this work the non-controllable input variables (i.e., nqq XXX ,,, 21   ) are considered to be fixed at some known values (i.e., nqq xxx ~,,~,~ 21   ) and )( XG is, consequently, a function of the controllable inputs only (i.e., q XXX ,,, 21  ). Given a current abnormal state of the system X , the most effective way to “escape” from it and reach safer conditions is to “perturb” the χ-th controllable variable X ,  q,1 , that is associated to the larger * j HV value, nj ,,2,1  , in order to attain the largest possible variation of )( XG (i.e., of the safety margin) by adding or subtracting a constant h to its nominal value j x̂ . When a perturbation of X does not result in a reduction of  )( XG (i.e., the expected safety margin value), it is considered that the system has reached normal conditions. Based on these knowledge rules, the ES algorithm performs the following four steps at each υ- th iteration of safety margins control, N,,1,0  (where 1N is the final number of iterations): 1. it computes: i. the realization )(  xgg  of G when the system is in its current nominal state ),,,( ,,2,1  qxxxx  ; ii. the expected value  )( XG over the range of input uncertainty, representing the system actual state       qqqq x,xxxx,xxxx,xxxX   ,,2,22,21,11,1  ; iii. * j HV of each input j X , qj ,,2,1  , over its own range   jjjj,j, xxxxX   ,, , with the initial ranges of all the other inputs X reduced to  ,,,2,1 ,,,,, qxxxx  ( q,,2,1  , j ). 7 2. it identifies the input variable X ,  q,1 , with the largest * j HV value (i.e.,   }{max: * , ,1 * ,  j qj HVHVX   that, once “perturbed”, generates the largest change in G ; 3. it computes two realizations  g and  g of function G in  x and  x , respectively, where: i. ),,,,,( ,,2,1  qxhxxxx    ; ii. ),,,,,( ,,2,1  qxhxxxx    ; iii. h is the distance between the current system nominal state x and the one that follows 1x (see below and Section 3 for further details). 4. it checks the following statements: i. if  gg ˆˆ   and ])([])([  XGXG   , then     xx 1 ; ii. otherwise, if  gg ˆˆ   and ])([])([  XGXG   , then     xx 1 . In other words, once the (χ-th) direction along which the largest change in G can be achieved is found (as in Step 3), Step 4 looks for the system nominal state that results in a decrease in G between  x and  x and it finally checks that ])([ XG is actually decreasing. The algorithm stops when neither one of the two above statements is verified because ])([])([  XGXG   and ])([])([  XGXG   , which means that moving the system from its nominal current state x to  x or to  x would not result in a decrease of ])([ XG . As it can be seen, the algorithm identifies a sequence of nominal states 1,,2,1 }{  Npp x  in S , obtained by shifting at each step υ the system nominal state x of a length h in the direction of the input variable ( X , ],1[ q ) that allows the largest reduction of  )( XGE (computed as in Step 1 Point iii). The choice of the step length h is a delicate task, because of two reasons: 1. the shorter h is, the more accurate is the final nominal state 1Nx , but the larger is the number of iterations N that the algorithm needs to perform to get to a final solution; 8 2. * ,jHV is computed at the υ-th step of the algorithm for each j-th input variable over the range ],[ , jjjj,j, xxxxRG   , whose length equals x2 ; it is, then, suggested that xh  because at step υ there is no information on the behavior of function g outside range j,RG . As ),,,( 21 q xxxg  is generally a complex function of q xxx ,,, 21  , the integrals in Eq. (7) are very complicated and cannot be solved analytically. Conversely, numerical techniques such as Monte Carlo integration [Zio, 2013], can allow approximating such integrals and computing ])([ XGE and ,jHV in a reasonable time at each υ-th iteration of the algorithm. In all generality, if we are trying to estimate the integral I of a function )( xvv  over an integration domain n  (i.e., xdxvI    )( ), where n x  , and a PDF )(xf X can be defined over  , then I is equivalent to [Zio, 2013]:          )( )( )( )( )( xf xv xdxf xf xv I X X X (10) where        )( )( xf xv is the expected value of the ratio )( )( xf xv computed over  with respect to a random variable distributed according to )(xf X . Eq. (10) is true for any PDF )(xf X on  , as long as 0)( xf X whenever 0)( xv . The value of        )( )( xf xv can, then, be estimated by: i) generating SN random samples )(c x (with SNc ,,2,1  ) of x according to )(xf X , ii) computing the ratio )( )( )( )( c X c xf xv for each sample and iii) finding the average SN I of these values:    S S N c c X c S N xf xv N I 1 )( )( )( )(1 (11) 9 According to the law of large numbers, as more and more samples are taken the Monte Carlo estimator SN I is guaranteed to converge to I , with a standard deviation that decreases as a function of the number of samples SN , by a rate of SN [Rubinstein, 1981; Kalos et al., 2008; Zio, 2013]. In this work, the pseudo-random samples )(c x are generated from the uniform distribution )(xu X , x , according to the algorithm originally proposed by von Neumann [Zio, 2013]. 3. ANALYTICAL CASE STUDY 3.1 Analytical Model Description Let us assume that the safety-significant response variable Y of the system is given by an analytical function m of two controllable inputs 1X and 2X as: 3 2 2 21 2 121 )1()2()3(8),(  XXXXXXmY (12) where  10,101 X and  10,102 X are random variables (such as, reactor control rods position, reactor coolant mass flow rate, steam generator feed-water mass flow rate, surge line mass flow rate, etc. [Cammi et al., 2006; Baraldi et al., 2013; Di Maio et al., 2015]), obeying two independent normal distributions )4,2(1N and )25.6,0(2N , respectively. The safety threshold limit to be applied to Y is 125 y γ and the resulting limit-state function of the system is [Bourinet et al., 2011; Zio et al., 2008]: 125)1()2()3(8),(),( 3 2 2 21 2 12121  XXXXγXXYXXG y (13) In Figures 1 and 2, the limit-state function G (continuous grid) and the system failure boundary F (solid line) are shown as functions of variables 1X and 2X . According to the definition of failure boundary given in Section 1, in this simplified case, F (Figures 1 and 2) can be easily found by the intersection between surface ),( 21 xxg (Figure 1) and plane 0g . In Figure 2, the safe ( S ) and failure ( F ) domains of the system are also shown. 10 Figure 1: limit-state surface G and failure boundary F of the analytical model m considered. Figure 2: failure boundary F of the analytical model m represented in the input space ]10,10[]10,10[  . 3.2 Assumptions As already discussed in Section 1, the uncertainty affecting input j X , 2,1j , with j X being fixed at some value  10,10ˆ  j x , is here modeled by considering j X to be ranging in the interval ]ˆ,ˆ[ jjjj xxxx  . For the sake of simplicity and without loss of generality, j x is independent of j x~ and it is computed as follows, for both 1X and 2X : 11 5.0   I jj j N ab x , 2,1j (14) where 10 j a and 10 j b are the lower and upper bounds of the domain ],[ jj ba of input j X , respectively, while 20IN is an arbitrary number of subintervals in which ],[ jj ba is partitioned. In a real case, the magnitude of j x would depend on the actual variability of j X and represents the confidence on how close j X actually is to its nominal value j x̂ (for instance, j X may describe the mass flow rate removed from the Steam Generator (SG) of a NPP by a set of safety valves during a “swell and shrink” accident [Di Maio et al., 2015], with j x̂ and j x equal to 25 [kg/s] and 5 [kg/s], respectively). Furthermore, the step length h is set equal to its maximum allowed value (Section 2) that is 0.5 in this case (for example, if j X represents the mass flow rate through the SG safety valves, h is the increase or decrease in j X resulting from regulating the valves operating position). 3.3 Results The methodological steps illustrated in Section 3.2 have been applied to five different initial nominal positions of the system )( 0  x , 5,,2,1  , in S . After applying on each )( 0  x the set of rules that constitutes the ES detailed in Section 2, five final nominal positions )( 1   N x have been obtained as reported in Table 1, where 1N is the final number of iterations needed to identify )( 1   N x . x 0,1 (ς) x 0,2 (ς) x Nς-1,1 (ς) x Nς-1,2 (ς) 1 -3.0 -2.0 0.0 -0.5 11 2 3.0 -2.0 0.0 -0.5 11 3 0.0 3.0 0.0 -0.5 15 4 8.0 6.0 9.0 10.0 12 5 0.0 6.0 9.0 10.0 28 N ς +1 Initial nominal states x 0 (ς) ς Final nominal states x Nς-1 (ς) 12 Table 1: initial ( )( 0  x ) and final ( )( 1   N x ) nominal states of the system in S , and number of performed iterations ( 1N ). From Table 1, one can read that )5.0,0.0( )3( 1 )2( 1 )1( 1 321   NNN xxx and ).0.10,0.9()5( 1 )4( 1 54   NN xx Indeed, depending on the initial nominal state of the system )( 0  x , )( 1   N x is either )5.0,0.0( P or )0.10,0.9(P , as )625.0,118.0( Q and P are two local minima of function g (see Figure 1). As an example, Figure 3 shows the optimal sequence 1,,2,1 )( }{    Npp x  of the system nominal states computed for 1 , where the system initial nominal position is )2,3()1( 0 x and the total number of iterations is 1111 N . To explain in more detail, Figure 4 illustrates each computed system state )( p x , 9,,2,1 p , coupled with its corresponding realization )( )1( pp xgg  of function G . As can be noticed, at each step the algorithm successfully identifies the direction along which function G varies the most between the 1X and 2X axes; then, it brings the system into a new optimal state, where p g is lower; finally, the algorithm stops when ])([ pXG cannot be reduced any further and the system is ultimately located in position )5.0,0.0( P . 13 Figure 3: sequence of system optimal nominal states 1,,2,1 )( }{    Npp x  for 1 . Figure 4: sequence of system optimal nominal states 1,,2,1 )( }{    Npp x  coupled with the corresponding realizations 1,,2,1 )( )}({     Nppp xgg  . For the sake of completeness, in Figure 5 the entire sequence of optimal nominal states of the system 1,21 )( }{   …,N,p=p x is illustrated for 5,,2,1  . Figure 5: sequences of system optimal nominal states 1,21 )( }{   …,N,p=p x computed for 5,,2,1  . 14 These results can be usefully compared with those obtained in [Di Maio et al., 2016] for the same analytical case study of Eq. (11). However, it should be honestly pointed out that these two works rest on two different rationales: i. in [Di Maio et al., 2016], the objective is to identify the farthest system state * x from the failure boundary F , regardless of the system working conditions when the accident scenario FE is triggered. Since * x is the farthest from F , it is also considered the optimal position for the system, i.e., the one where the final failure event is less likely to occur; ii. in this work, the initial abnormal state of the system 0 X is of paramount importance, as it is required to assess the feasibility of “escaping away” from failure starting from it. Indeed, the proposed ES evaluates where the system should go step-by-step in the input space depending on the behavior of function G , which is representative of the system safety margin. In [Di Maio et al., 2016], the only non-controllable variable y  (i.e., the system safety threshold) is assumed to be normally distributed on the range [−500,2500], with mean 500 y  and variance 2500 2  y  . Based on the randomness of variable y  , various safest operating conditions (i.e., states) ),( * 2 * 1 * xxx  are recommended for the system, depending on the operator risk attitude: i. risk-averse: )00.10,63.4( * x ; ii. risk-prone: )00.10,00.8( * x ; iii. mean: )95.9,91.5( * x ; iv. median: )00.10,25.4( * x . where the risk-prone and risk-averse conditions are defined as the 80-th and 20-th percentiles of a population of safest operating conditions computed for different values of y  , respectively. 15 In this case study, y  is supposed to be known and set equal to 125. As Table 1 and Figure 5 show, for 3,2,1 the final optimal nominal state )( 1   N x of the system (i.e., )5.0,0.0( P ) is far from all the solutions * x of [Di Maio et al., 2016]. This is due to the different criterion adopted in determining 1,21 )( }{   …,N,p=p x and ),( * 2 * 1 * xxx  : in [Di Maio et al., 2016] the (possibly unreachable) safest state * x is identified as the input variables setting for which the system is best sheltered against failures, whereas, in this work we provide the plant operator with instructions on how to change step by step the abnormal variables from the setting 0x to a normal state 1N x , which may (or may not) correspond to * x (that might not be realistic to reach for some accident scenarios FE ), depending on the actual behavior of function G . 4. PRESSURIZER CASE STUDY 4.1 Pressurizer Model Description The pressurizer of a PWR NPP, whose scheme is shown in Fig. 6, has been simulated through a Matlab SIMULINK model that is here taken as the TH model of the component [Baraldi et al., 2013]. The developed TH model is based on the application of the mass and energy conservation equations to the two regions of vapor and liquid. The exchanges between the two regions, due to evaporation of the liquid and condensation of the gas, have been taken into account [Kuridan et al., 1998, Todreas et al., 1990]. The system of non-linear differential equations describing the TH model is detailed in [Baraldi et al., 2010]. 16 Fig. 6: Scheme of a PWR NPP pressurizer. Let us assume the safety-relevant response Y of the system to be the pressure reached within the pressurizer when two controllable variables X1 and X2 (i.e., sprayers mass flow rate and heaters power, respectively) are the levels under control of the operators to counteract any developing surge line flow rate excursion scenario. To estimate g, a data set has been built in order to represent a realistic situation simulating 100 surge line transients, whose initial conditions have been taken to represent the realistic situation of a standard PWR pressurizer, (Table 2), whereas the surge mass flow rate excursion has been randomly sampled from the uniform distribution in the range [-10,10] kg/s. Initial condition Level 7.221 m Liquid temperature 342.1 °C Vapor temperature 342.3 °C Pressure 150.0 bar Table 2: Initial conditions of the pressurizer. 17 If 𝑋1 ∈ [0,7] kg/s is distributed like a N~(3.5, 2.4) and 𝑋2 ∈ [0,320] kW is distributed like a N~ (160,107.650), the resulting Y fitted to the available dataset is equal to Eq. (15). 𝑌 = 𝑚(𝑋1,𝑋2) = 160.6 − 11.2𝑋1 + 10 −2𝑋2 + 2.5 · 𝑋1 2 + 2.9 · 10−3𝑋1𝑋2 − 0.1 · 𝑋1 3 (15) The safety threshold γy is set equal to 155 bar and the resulting limit-state function G is given in Eq. (16). 𝐺(𝑋1,𝑋2) = 𝑌(𝑋1,𝑋2) − 𝛾𝑦 = 5.6 − 11.2𝑋1 + 10 −2𝑋2 + 2.5 · 𝑋1 2 + 2.9 · 10−3𝑋1𝑋2 − 0.1 · 𝑋1 3 (16) In Figures 7 and 8, the limit-state function G (continuous grid) and the system failure boundary 𝜕𝐹 (solid line) are shown as functions of variables X1 and X2. In Figure 8, the safe (S) and failure (F) domains are also shown. Figure 7: limit-state function G and failure boundary F of the pressurizer model m. 18 Figure 8: failure boundary F of the pressurizer model m represented in the input space ]320,0[]7,0[  . 4.2 Assumptions The uncertainty affecting the input Xj, j = 1,2, is here modelled considering Xj to be ranging in the interval [𝑥�̂� − ∆𝑥𝑗,𝑥�̂� + ∆𝑥𝑗]. The ∆𝑥𝑗 values to be used in what follows are listed in Table 3. j aj bj NI Δxj 1 0 7 18 0.39 2 0 320 18 17.8 Table 3: Values of Δxj. The step length h is set equal to its maximum allowed value (Section 2) for each one of the variables, i.e. h1 = 0.78 kg/s and h2 = 35.6 kW (with jX representing the mass flow rate through 19 the sprayers of the pressurizer and the power of its heaters respectively, j h is the increase or decrease in j X resulting from regulating the two different systems under operator’s control). 4.4 Results The methodological steps have been applied to three different initial nominal positions of the system �̅�0 (𝜍) , ς = 1, 2, 3, that belong to S and are close to 𝜕𝐹. The ES has been applied to these initial positions and the final system nominal states �̅�𝑁𝜍−1 (𝜍) have been obtained as shown in Table 4, where Nς+1 is the number of iterations needed to reach �̅�𝑁𝜍−1 (𝜍) . ς Initial nominal states �̅�0 (𝜍) Final nominal states �̅�𝑁𝜍−1 (𝜍) Nς+1 �̅�0,1 (𝜍) �̅�0,2 (𝜍) �̅�𝑁𝜍−1,1 (𝜍) �̅�𝑁𝜍−1,2 (𝜍) 1 0.79 0 3.89 0 5 2 0.79 35.6 3.89 0 6 3 1.56 320 3.89 0 13 Table 4: initial (�̅�0 (𝜍) ) and final (�̅�𝑁𝜍−1 (𝜍) ) states of the system in S, and number of performed iterations (Nς+1) At each step, the algorithm identifies the direction X1 or X2 along which the function G mostly decreases, leading the system into a new state where gp is lower. The algorithm automatically stops in the state in which the minimum value of gp is found. From Table 4, it can be seen that, irrespective of the initial nominal state, the ES leads the operator to the same final state �̅�𝑁𝜍−1 (𝜍) ≡ (3.89,0) ≡ 𝑃𝜗, which corresponds to a local minima of the limit-state function G. Without loss of generality, Table 5 shows the sequence of optimal nominal states for ς = 3, in which the final number of iterations is N3+1=13. p �̅�1 (3) �̅�2 (3) 1 1.56 320 2 2.33 320 3 3.11 320 4 3.89 320 20 5 3.89 284 6 3.89 249 7 3.89 213 8 3.89 178 9 3.89 142 10 3.89 107 11 3.89 71.1 12 3.89 35.6 13 3.89 0 Table 5: sequence of system optimal nominal states 1,,2,1 )( }{    Npp x  for 3 . In Figure 9, it is noteworthy comparing the pressure transient arising from the ES with the one deriving from the application of a PID control system, when the same in-surge accident scenario occurs (e.g., a total water mass entering in the pressurizer from the primary system of 1466 kg, with flow rate of 9 kg/s). Figure 9: Pressure transient in an in-surge accident scenario using the ES (dashed line) and the PID controller (thick solid line) to control the system. The failure boundary is represented by the horizontal thin line. In the first phase of the transient (up to 30 s), the PID system shows a constant pressure of 150 bar, while the ES produces a drop of 0.5 bar. The reason is that the safest position is 147.9 bar, 21 as it can be seen from the data at the end of the transient; thus, the ES works to reach this condition. At t = 30 s, the in-surge accident takes place and it is possible to notice a rise in pressure in both transients, steeper for the PID controller. But the main observation is that in the PID transient the pressure overcomes the failure boundary of 155 bar, whereas the ES manages to keep the pressure well below the failure boundary for the entire duration of the in-surge accident scenario. Lastly, as the accident scenario is recovered, the actions driven by the ES lead to a reduction in pressure, till it reaches its safest optimal working level of 147.9 bar. Finally, it is worth mentioning that the capability of the proposed ES in providing the plant operator with a live, effective and unequivocal AOP to keep the safety margins and avoid system failure only depends on the computational demand of the TH code utilized, and not on the ES algorithm itself: more complex simulation codes, that may become computationally intensive, may challenge the practical scalability of the proposed procedure and endanger the live guidance feature offered by the ES, despite that these could still be utilized for an off-line design of accurate, effective and unequivocal AOPs. 5. CONCLUSIONS In this work, an ES based on regional sensitivity indices (i.e., the Revised Variance Ratio Functions) has been proposed for improving the safety margin of a system when an IE occurs during operation. A RSA is performed to define the strategy for keeping under control the uncertainty affecting the controllable input variables. Such uncertainty has been modeled by centering the nominal value of the input variable on a reduced range of its domain. The ES successfully achieves the stepwise maximization of the system safety margin by acting on the model controllable inputs. The proposed approach has been proved effective on an analytical case study, which mimics the response of a NPP to an accident scenario, and on a case study of a pressurizer exposed to an accidental surge line transient. 22 REFERENCES [Aldemir et al., 2001] T. Aldemir, P. Wang, D.W. Miller. Parameter and State Estimation Using DSD, Transactions of the American Nuclear Society, volume 84; Jun 2001. [Apostolakis, 1990] G. E. Apostolakis. The Concept of Proability in Safety Assessments of Technological Systems. In Science, Vol. 250, pp. 1359-1364. December 1990. [Baraldi et al., 2010] P. Baraldi, A. Cammi, F. Mangili and E. Zio. An ensemble approach to sensor fault detection and signal reconstruction for nuclear system control. In Annals of Nuclear Energy, Vol. 37, pp. 778-790. Elsevier, April 2010. [Baraldi et al., 2013] P. Baraldi, F. Di Maio and E. Zio. Unsupervised Clustering for Fault Diagnosis in Nuclear Power Plant Components. In International Journal of Computational Intelligence Systems, pp. 764-777. Taylor and Francis, Ltd., 2013. [Bolado-Lavin et al., 2009] R. Bolado-Lavin, W. Castaings and S. Tarantola. Contribution to the Sample Mean Plot for Graphical and Numerical Sensitivity Analysis. In Reliability Engineering & System Safety, Vol. 93, pp. 1563-1573. Elsevier, 2009. [Bourinet et al., 2011] J. M. Bourinet, F. Deheeger and M. Lemaire. Assessing Small Failure Probabilities by Combined Subset Simulation and Support Vector Machines. In Structural Safety, Vol. 33, pp. 343-353. Elsevier, June 2011. [Buchanan et al., 1988] B. G. Buchanan and R.G. Smith. Fundamentals of Expert Systems. In Annual Review of Computer Science, Vol. 3, pp. 23-58. Annual Reviews Inc., June 1988. [Cammi et al., 2006] A. Cammi, L. Luzzi, A. A. Porta and M. E. Ricotti. Modelling and Control Strategy of the Italian LBE-XADS. In Progress in Nuclear Energy, Vol. 48, pp. 578- 589. Elsevier, 2006. [Di Maio et al., 2015] F. Di Maio, M. Vagnoli and E. Zio. Risk-Based Clustering for Near Misses Identification in Integrated Deterministic and Probabilistic Safety Analysis. In 23 Science and Technology of Nuclear Installations. Hindawi Publishing Corporation, February 2015. [Di Maio et al., 2016] F. Di Maio, A. Bandini, E. Zio, A. Alfonsi and C. Rabiti. An Approach Based on Support Vector Machines and a K-D Tree Search Algorithm for Identification of the Failure Domain and Safest Operating Conditions in Nuclear Systems. In Progress in Nuclear Energy, Volume 88, pp. 297-309, 2016. [Hakobyan et al., 2008] A. Hakobyan, T. Aldemir, R. Denning, S. Dunagan, D. Kunsman, B. Rutt, U. Catalyurek, "Dynamic generation of accident progression event trees," Nuclear Engineering and Design, 238, 3457-3467 (2008) [Helton et al., 1996] J.C. Helton and D. E. Burmaster. Guest Editorial: Treatment of Aleatory and Epistemic Uncertainty in Performance Assessments for Complex Systems. In Reliability Engineering & System Safety, Vol. 54, pp. 91-94. Elsevier, 1996. [Hines et al., 1998] J.W Hines, R.E. Uhrig, “Use of autoassociative neural networks for signal validation”, Journal of Intelligent and Robotic Systems, vol. 21, pp. 143–154, 1998. [Hsieh et al., 2012] M. H. Hsieh, S. L. Hwang, K. H. Liu, S. F. M. Liang and C. F. Chuang. A Decision Support System for Identifying Abnormal Operating Procedures in a Nuclear Power Plant. In Nuclear Engineering and Design, Vol. 249, pp. 413-418. Elsevier, April 2012. [Hurtado et al., 2004] J. E. Hurtado. Structural Reliability – Statistical Learning Perspectives. Lecture Notes in Applied and Computational Mechanics. Springer, 2004. [IAEA, 2000] VV. AA. IAEA Management of Aging of I&C Equipment in Nuclear Power Plant. IAEA-TECDOC-1147. International Atomic Energy Agency, 2000. [IAEA, 2004] VV. AA. Management Life Cycle and Aging at Nuclear Power Plants: Improved I&C Maintenane. IAEA-TECDOC-1402. International Atomic Energy Agency, 2004. 24 [IAEA, 2005] VV. AA. Effective Corrective Actions to Enhance Operational Safety of Nuclear Installations. IAEA-TECDOC-1458. Internation Atomic Energy Agency, July 2005. [IAEA, 2006] VV. AA. Development and Review of Plant Specific Emergency Operating Procedures. IAEA Safety Reports Series, No. 48. International Atomic Energy Agency, February 2006. [Ikram et al., 2015] A. Ikram and U. Qamar. Developing an Expert System Based on Association Rules and Predicate Logic for Earthquake Prediction. In Knowledge-Based Systems, Vol. 75, pp. 87-103. Elsevier, 2015. [Kalos et al., 2008] M. H. Kalos and P. A. Whitlock. Monte Carlo Methods. Wiley, 2008. [Kim et al., 2001] J. T. Kim, K. C. Kwong, I. K. Hwang, D. Y. Lee, W. M. Park, J. S. Kim and S. J. Lee. Development of Advanced I&C in Nuclear Power Plants: ADIOS and ASICS. In Nuclear Engineering and Design, Vol. 207, pp. 105-119. Elsevier, 2001. [Kim et al., 2007] J. H. Kim and P. H. Seong. The Effect of Information Types on Diagnostic Strategies in the Information Aid. In Reliability Engineering & System Safety, Vol. 92, pp. 171-186. Elsevier, February 2007. [Kwon et al., 1999] K. C. Kwon and J. H. Kim. Accident Identification in Nuclear Power Plants Using Hidden Markov Models. In Engineering Applications of Artificial Intelligence, Vol. 12, pp. 491-501. Elsevier, 1999. [Lee et al., 2007] S. J. Lee, K. Mo, P. H. Seong. Development of an Integrated Decision Support System to Aid the Cognitive Activities of Operators in Main Control Rooms of Nuclear Power Plants. In Proceedings of IEEE Symposyum on Computational Intelligence in Multicriteria Decision Making (MCDM), pp. 146-152. IEEE, 2007. [Liao, 2005] S. H. Lao. Expert System Methodologies and Applications – A Decade Review from 1995 to 2004. In Expert Systems with Applications, Vol. 28, pp. 93-103. Elsevier, 2005. 25 [Lucia et al., 2004] D. J. Lucia, P. S. Beran and W. A. Silva. Reduced-Order Modeling: New Approaches for Computational Physics. In Progress in Aerospace Sciences, Vol. 40, pp. 51-117. Elsevier, February 2004. [Ma et al., 2011] J. Ma and J. Jiang. Applications of Fault Detetion and Diagnosis Methods in Nucelar Power Plants: A Review. In Progress in Nuclear Energy, Vol. 53, pp. 255- 266. Elsevier, January 2011. [McBride et al., 1993] R. D. McBride and D. E. O’Leary. The Use of Mathematical Programming with Artificial Intelligence and Expert Systems. In European Journal of Operational Research, Vol. 70, pp. 1-15. Elsevier, October 1993. [Oberkampf et al., 2004] W. L. Oberkampf, J. C. Helton, C. A. Joslyn, S. F. Wojtkiewicz and S. Ferson. Challenge Problems: Uncertainty in System Response Given Uncertain Parameters. In Reliability Engineering & System Safety, Vol. 85, pp. 11-19. Elsevier, September 2004. [Park et al., 2015] J. Park and W. Jung. A Systematic Framework to Investigate the Coverge of Abnormal Operating Procedures in Nuclear Power Plants. In Reliability Engineering & System Safety, Vol. 138, pp. 21-30. Elsevier, February 2015. [Rubinstein, 1981] R. Y. Rubinstein. Simulation and the Monte Carlo Method. Wiley, 1981. [Tarantola et al., 2012] S. Tarantola, V. Kopustinskas, R. Bolado-Lavin, A. Kaliatka, E. Uspuras and M. Vaisnoras. Sensitivity Analysis Using Contribution to Sample Variance Plot: Application to a Water Hammer Model. In Reliability Engineering & System Safety, Vol. 99, pp. 62-73. Elsevier, March 2012. [Todreas and Kazimi, 1990] N. Todreas, M. S. Kazimi. Nuclear Systems 1, pp. 239- 287. Taylor & Francis, 1990. 26 [Wei et al., 2014] P. Wei, Z. Lu, W. Ruan and J. Song. Regional Sensitivity Analysis Using Revised Mean and Variance Ratio Functions. In Reliability Engineering & System Safety, Vol. 121, pp. 121-135. Elsevier, 2014. [Yoshikawa, 2005] H. Yoshikawa. Human-Machine Interaction in Nuclear Power Plants. In Nuclear Engineering and Technology, Vol. 79, pp. 151-158. Elsevier, 2005. [Zio et al., 2008] E. Zio and F. Di Maio. Bootstrap and Order Statistics for Quantifying Thermal-Hydraulic Code Uncertainties in the Estimation of Safety Margins. In Science and Technology of Nuclear Installations, Vol. 2008, Article ID 340164. Hindawi Publishing Corporation, 2008. [Zio et al., 2010] E. Zio, F. Di Maio and J. Tong. Safety Margins Confidence Estimation for a Passive Residual Heat Removal System. In Reliability Engineering & System Safety, Vol. 95, pp. 828-836. Elsevier, August 2010. [Zio, 2013] E. Zio. The Monte Carlo Simulation Method for System Reliability and Risk Analysis. Springer Series in Reliability Engineering. Springer, 2013. [F. Di Maio et al., 2013] F. Di Maio, P. Baraldi, E. Zio and R. Seraoui. Fault Detection in Nuclear Power Plants Components by a Combination of Statistical Methods. In IEEE Transactions on Reliability, Vol. 62, pp. 833-845. IEEE, October 2013. [S. Al-Dahidi] S. Al-Dahidi, P. Baraldi, F. Di Maio and E. Zio. A Novel Fault Detection System Taking Into Account Uncertainties in the Resconstructed Signals. In Annals of Nuclear Energy, Vol. 73, pp. 131-144. Elsevier, November 2014.