Biasing transition rate method for probabilistic safety assessment by direct M-C simulation

Biasing transition rate method for probabilistic safety assessment by direct M-C simulation 1 first-author 2 2 2 2 2 2 corresp Corresponding author. E-mail address: jin.wang@fds.org.cn Corresponding author. E-mail address: jin.wang@fds.org.cn jin.wang@fds.org.cn University of Science and Technology of China, Hefei 230027, China University of Science and Technology of China, Hefei 230027, China Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academic of Sciences, Hefei 230031, China Key Laboratory of Neutronics and Radiation Safety, Institute of Nuclear Energy Safety Technology, Chinese Academic of Sciences, Hefei 230031, China Direct Monte Carlo (MC) simulation is a powerful probabilistic safety assessment (PSA) method for accounting dynamics of the system. But it is not efficient at simulating rare events. A biasing transition rate method based on direct MC simulation is proposed to solve the problem in this paper. This method biases transition rates of the components by adding virtual components to them in series to increase the occurrence probability of the rare event, hence the decrease of the variance of MC estimator. Several cases are used to benchmark this method. The results show that the method is effective at modeling system failure and is more efficient at collecting evidence of rare events than the direct MC simulation. The performance is greatly improved by the biasing transition rate method. Direct Monte Carlo simulation Probabilistic safety assessment Biasing transition rate method 1 Introduction 1 Probabilistic safety assessment (PSA) is performed by regulatory bodies to check whether the designs of nuclear power plants comply with regulatory requirements, and by industry for identifying key vulnerabilities [ 1 1 - 4 4 1-4 ]. Traditional PSA methods, i.e. fault trees (FTs) and event trees (ETs) have critical limitations in practice, and the results may differ widely from real values due to imprecise descriptions of component aging and maintenance and binary modeling of component behavior (only faulty/safe states are considered) and the neglect of dynamics of the system (i.e., effect of the order and timing of failure events on the accident progression) [ 2 2 , 5 5 , 6 6 ]. In the PSA analysis of some advanced nuclear systems, such as the China Lead-based Research Reactor (CLEAR) [ 7 7 ], the fusion-driven subcritical system [ 8 8 ] and the International Thermonuclear Experimental Reactor (ITER) [ 9 9 ], high accuracy is required and such factors cannot be neglected. Monte Carlo (MC) simulation is a popular dynamic method developed to overcome these limitations, [ 5 5 , 10 10 - 12 12 10-12 ], but it is not efficient at simulating rare events in a complex system [ 13 13 ]. Although sampling algorithms of antithetic variable sampling [ 14 14 ], dagger sampling [ 15 15 ], and stratified sampling [ 16 16 ], can be applied to MC simulation to reduce the estimate error and the computational cost, the results are still unsatisfactory. Biasing techniques can solve these problems effectively by forcing the events of interest to occur more frequently in the MC simulation. Some biasing techniques, such as importance sampling [ 17 17 - 20 20 17-20 ] and multicanonical Monte Carlo (MMC) [ 21 21 , 22 22 ], can reduce the estimate error and the computational cost. To improve the simulation efficiency, in this paper we propose a new biasing method, the biasing transition rate method. Its idea is to bias transition rates of the components by adding virtual components to them in series to increase the probability of occurrence of the rare event. The performance is compared to that of the sole use of direct MC method. 2 The method 1 2.1. Assumptions 2 The biasing transition rate method is based on the following assumptions: 1) The system consists of l components, 1,..., l ; 2) The components’ failure time and repair time follow an exponential distribution or a Weibull distribution, or linear aging based on these two distributions can be considered; 3) Every system state is possible at any time; 4) The system is s-coherent, each of its working components is monotonously beneficial to the system; and 5) The unexpected event is a rare event, 0 < P u < 1/2 is the probability of the unexpected event. 2.2. Biasing transition rate method 2 Let us consider a simple case in which a component transition occurs with a low transition rate λ . The idea of biasing transition rate is to add a virtual component with a transition rate nλ to the component in series, so that transition rate of the integration of the real component and virtual component increases, and with Monte Carlo sampling it is easier to achieve the transition within the mission time of T m in each trial. The sequence branches into two at the transition time <math id="M1"><mrow><msub><mi>T</mi><mi>t</mi></msub></mrow></math> of the integration: in the first branch, the real component changes its state (i.e. from safe to faulty) and the probability (weight) of this sequence equals to the product of the contribution rate of the real component (1/( n +1)) and the probability of the mother sequence P m ; the second branch keeps the original state, and the probability is the product of the contribution rate of the virtual component ( n /( n +1)) and the probability of the mother sequence. The sequences are simulated continuously until the absorbing state, or T m , is reached. The transition of the integration occurs with much greater probability in each trial. Fig. 1 Fig. 1 is schematics of this method of component simulation. The red dash lines denote the sequences contributed by the real component, and the probability of each sequence is shown on the right. The flowchart of this method is shown in Fig. 2 Fig. 2 . In a multi-component system simulation, the probability that the rare event occurs can be significantly increased by biasing the transition rates of some selected components. The probability (weight) of the j th sequence generated in the i th trial at mission time T m , can be expressed as <math display="block" id="M2"><mrow><msub><mi>p</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mtext> </mtext><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><msubsup><mo>∑</mo><mrow><mi>s</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>q</mi><mi>i</mi><mi>j</mi></mrow></msubsup><mi>θ</mi><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>,</mo><mtext> </mtext><msub><mi>q</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>≥</mo><mn>1</mn></mrow></mtd></mtr><mtr><mtd><mrow><mn>0</mn><mo>,</mo><mtext>              </mtext><msub><mi>q</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mn>0</mn></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> where q ij is the number of branch points in the j th sequence generated in i th trial, and <math display="block" id="M3"><mrow><mo> </mo><mo> </mo><mi>θ</mi><mrow><mo>(</mo><mi>s</mi><mo>)</mo></mrow><mo>=</mo><mo> </mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mo> </mo><mfrac><mn>1</mn><mrow><mrow><mo>(</mo><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo><mtext> the </mtext><mi>s</mi><mo>−</mo><mtext>th transition caused by the real component</mtext></mrow></mtd></mtr><mtr><mtd><mrow><mfrac><mi>n</mi><mrow><mrow><mo>(</mo><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo><mtext> otherwise</mtext></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> where n is the biasing factor, a non-negative number. A greater n means a greater number of branches created in each trial and longer time of computation. Given the increase in computational cost, it is suggested that the sum of the failure probabilities of all the biased components with the biased rates in the mission time lies below 0.1. Four transition rates are considered in this paper: 1) An exponential distribution for the component transition ( λ =λ 0 , λ 0 is the design transition rate, a constant) 2) A Weibull distribution for the component transition ( λ =αλ 0 t α −1 ) 3) A linear aging model based on an exponential distribution ( λ = λ 0 + k t , k is the aging factor) 4) A linear aging model based on a Weibull distribution ( λ = αλ 0 t α −1 + k t , <math id="M4"><mrow><mtext> </mtext><mi>α</mi></mrow></math> is a factor in the Weibull distribution) As an example, assuming that the transition time of a component follows an exponential distribution with failure rate λ , its probability of failure before time t can be expressed as <math display="block" id="M5"><mrow><mi>F</mi><mo>=</mo><mn>1</mn><mo>−</mo><msup><mi>e</mi><mrow><mo>−</mo><mi>λ</mi><mi>t</mi></mrow></msup></mrow></math> Assuming that the failure rate of the added virtual component is nλ , taking the two components as an integration, the probability of the integration failing before time t is equal to <math display="block" id="M6"><mrow><mi>F</mi><mo>=</mo><mn>1</mn><mo>−</mo><msup><mi>e</mi><mrow><mo>−</mo><mrow><mo>(</mo><mrow><mi>n</mi><mo>+</mo><mn>1</mn></mrow><mo>)</mo></mrow><mi>λ</mi><mi>t</mi></mrow></msup></mrow></math> The mean time to transition is reduced from 1/ λ to1/[( n +1) λ ]. The biasing parameters for the four transition rates are given in Table 1 Table 1 . * For Cases 1–4, the contribution rate θ is 1/( n +1) and n /( n +1) for the real and virtual components, respectively. 1 1 Cases 1 1 Real component 1 1 Virtual component 1 1 Integration 1 1 1 1 1 λ 0 1 1 nλ 0 1 1 ( n +1) λ 0 1 1 2 1 1 αλ 0 t α −1 1 1 nαλ 0 t α −1 1 1 ( n +1) αλ 0 t α −1 1 1 3 1 1 λ 0 + k t 1 1 n ( λ 0 + k t ) 1 1 ( n +1) ( λ 0 + k t ) 1 1 4 1 1 αλ 0 t α −1 + k t 1 1 n ( αλ 0 t α −1 + k t ) 1 1 ( n +1) ( αλ 0 t α −1 + k t ) To state advantages of this method, the estimator variation is theoretically analyzed in Section 2.3. Some parameters to evaluate performance of this method, such as the root mean square deviation (RMSD), the efficiency in collecting evidence of system failure, and the figure of merit, are analyzed in benchmark cases in Section 3. 2.3. Estimator based on the biasing transition rate method 2 Let <math id="M7"><mrow><msub><mi>z</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> be the j th sequence generated in the i th trial of the biasing transition rate simulation. The failure probability of the system Q can be estimated by the estimator <math id="M8"><mrow><msub><mi>Q</mi><mtext>b</mtext></msub></mrow></math> based on the biasing transition rate method, <math display="block" id="M9"><mrow><msub><mi>Q</mi><mtext>b</mtext></msub><mo>=</mo><mo> </mo><mfrac><mn>1</mn><mi>N</mi></mfrac><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><msubsup><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>M</mi><mi>i</mi></msub></mrow></msubsup><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>z</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow></math> where, N is the number of MC cycles; <math id="M10"><mrow><msub><mi>M</mi><mi>i</mi></msub></mrow></math> is the sequence number of the i th trial; and <math id="M11"><mi>ψ</mi></math> is a discrete function expressing the unexpected event, which can be expressed as Eq.(6) for the biasing transition rate method, <math display="block" id="M12"><mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>z</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow><mo>=</mo><mo> </mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>p</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>,</mo><mtext> if the unexpected event occurs</mtext></mrow></mtd></mtr><mtr><mtd><mrow><mo> </mo><mo> </mo><mtext>0, otherwise</mtext></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> The variation of Q b can be expressed by <math display="block" id="M13"><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mrow><mo>{</mo><mrow><msub><mi>Q</mi><mi>b</mi></msub></mrow><mo>}</mo></mrow><mo>=</mo><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><msubsup><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>M</mi><mi>i</mi></msub></mrow></msubsup><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>z</mi><mi>i</mi></msub><msub><mrow/><mi>j</mi></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mo stretchy="false">(</mo><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><msubsup><mo>∑</mo><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>M</mi><mi>i</mi></msub></mrow></msubsup><mi>ψ</mi><mo stretchy="false">(</mo><msub><mi>z</mi><mi>i</mi></msub><msub><mrow/><mi>j</mi></msub><mo stretchy="false">)</mo><mo stretchy="false">)</mo><mo stretchy="false">)</mo></mrow></math> The failure probability of the system Q can be estimated by the direct MC estimator <math id="M14"><mrow><msub><mi>Q</mi><mtext>d</mtext></msub></mrow></math> , <math display="block" id="M15"><mrow><msub><mi>Q</mi><mtext>d</mtext></msub><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow></math> where <math id="M16"><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow></math> donates the sequence in the i th trial. The variance of the estimator <math id="M17"><mrow><mtext> </mtext><mi>V</mi><mi>a</mi><mi>r</mi><mrow><mo>{</mo><mrow><msub><mi>Q</mi><mtext>d</mtext></msub></mrow><mo>}</mo></mrow></mrow></math> is <math display="block" id="M18"><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mrow><mo>{</mo><mrow><msub><mi>Q</mi><mtext>d</mtext></msub></mrow><mo>}</mo></mrow><mo>=</mo><mi>V</mi><mi>a</mi><mi>r</mi><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mo> </mo><mo> </mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mrow><mo>(</mo><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow></mrow></math> For the direct MC method, the discrete function expressing the unexpected event can be expressed as <math display="block" id="M19"><mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mo> </mo><mrow><mo>{</mo><mrow><mtable><mtr><mtd><mrow><mn>1</mn><mo>,</mo><mtext> if the unexpected event occurs</mtext></mrow></mtd></mtr><mtr><mtd><mrow><mn>0</mn><mo>,</mo><mtext> otherwise</mtext></mrow></mtd></mtr></mtable></mrow></mrow></mrow></math> Since E <math id="M20"><mrow><mo stretchy="false">(</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>M</mi><mi>i</mi></msub></mrow></munderover><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>z</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow><mo stretchy="false">)</mo><mtext> </mtext><mo>=</mo><mtext> </mtext><mi>E</mi><mo stretchy="false">(</mo><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo stretchy="false">)</mo></mrow></math> , and 0 < P u < 1/2, we have <math display="block" id="M21"><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mrow><mo>(</mo><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>M</mi><mi>i</mi></msub></mrow></munderover><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>z</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>≤</mo><mi>V</mi><mi>a</mi><mi>r</mi><mrow><mo>(</mo><mrow><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow></mrow></math> and <math display="block" id="M22"><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mrow><mo>{</mo><mrow><msub><mi>Q</mi><mi>b</mi></msub></mrow><mo>}</mo></mrow><mo>=</mo><mo> </mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mrow><mo>(</mo><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>j</mi><mo>=</mo><mn>1</mn></mrow><mrow><msub><mi>M</mi><mi>i</mi></msub></mrow></munderover><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>z</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mo>)</mo></mrow><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow><mo>≤</mo><mfrac><mn>1</mn><mrow><msup><mi>N</mi><mn>2</mn></msup></mrow></mfrac><mrow><mo>(</mo><mrow><munderover><mstyle mathsize="140%" displaystyle="true"><mo>∑</mo></mstyle><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></munderover><mi>V</mi><mi>a</mi><mi>r</mi><mo stretchy="false">(</mo><mi>ψ</mi><mrow><mo>(</mo><mrow><msub><mi>c</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo stretchy="false">)</mo></mrow><mo>)</mo></mrow><mo>=</mo><mi>V</mi><mi>a</mi><mi>r</mi><mrow><mo>{</mo><mrow><msub><mi>Q</mi><mi>d</mi></msub></mrow><mo>}</mo></mrow></mrow></math> This proves that the biasing transition rate method can decrease the variance of the MC estimator. 3 Benchmark cases 1 3.1. Description of the system 2 To illustrate performance of the method, a system consisting of three components, as shown in Fig. 3 Fig. 3 , is used to benchmark the biasing transition rate method. The following four cases are considered: 1) an exponential distribution for the component life; 2) a Weibull distribution for the component life; 3) linear aging based on an exponential distribution; and 4) an exponential distribution for the component life and for the component repair time. The mission timeis T m = 1000 h. For all the cases, the design failure rates are λ 1 = λ 2 = 10 −4 and λ 3 = 10 −5 ; the biasing factor is n = 9; and number of MC cycles is N = 1000. For Case 2, the parameter factors are α 1 = α 2 = α 3 =1.1. For Case 3, the aging factors are k 1 = k 2 =10 −8 and k 3 =10 −9 . For Case 4, design repair rates are u 1 = u 2 =10 −3 and u 3 =0. 3.2 Results and discussion 2 Figs. 4 Figs. 4 shows the results of the biasing transition rate simulation and the direct MC simulation, compared with the reference results (i.e. the results using the minimal cut set method and the result of a direct MC simulation with a huge sampling size). The results of the two kinds of models are almost identical, which shows that the biasing transition rate method is effective at estimating reliability in the four cases. Fig.4 Fig.4 also shows that the estimates provided by the biasing transition rate method are smoother than those provided by the direct MC simulation, because the former can collect more evidences of system failure than the latter in the same number of MC cycles. Table 2 Table 2 lists the RMSDs and the efficiency for collecting evidence of the failure events for the two methods applied to the four cases. RMSD is a measure of the error between the estimated and reference values. To compute the RMSD for the two methods, the results of minimal cut set models (Cases 1–3) and the result of direct MC models with a huge number of MC cycles (Case 4) are set as the reference values. The results show that the RMSD value of the biasing transition rate model is much smaller than that of the direct MC model, hence the increased closeness to the real value. This matches the proof in Section 2.3: the biasing transition rate method is much more efficient at collecting evidence of system failure, a rare event, than the direct MC method, because the probability increases greatly for the transition of the integrated real and virtual components, and so does the probability that the system failure occurs in each trial in the biasing transition rate simulation. For example, in Case 1), the RMSD values of the direct MC model and the biasing transition rate model are 0.00647 and 0.00034, respectively, and the expected times to ‘ system failure ’ are 0.2031 and 0.0064, respectively. * DMC, direct MC; BTR, biasing transition rate model. 1 1 Case No. 2 1 Computational cost /s 2 1 Number of evidence 2 1 RMSD 2 1 Expected time to ‘ system failure ’ /s 1 1   1 1 DMC 1 1 BTR 1 1 DMC 1 1 BTR 1 1 DMC 1 1 BTR 1 1 DMC 1 1 BTR 1 1 1 1 1 5.6865 1 1 7.7354 1 1 28 1 1 1213 1 1 0.00647 1 1 0.00034 1 1 0.2031 1 1 0.0064 1 1 2 1 1 6.1644 1 1 12.6543 1 1 44 1 1 4640 1 1 0.00482 1 1 0.00204 1 1 0.1401 1 1 0.0027 1 1 3 1 1 4.6761 1 1 7.0110 1 1 11 1 1 1315 1 1 0.00421 1 1 0.00057 1 1 0.4251 1 1 0.0053 1 1 4 1 1 4.2993 1 1 8.5845 1 1 12 1 1 1042 1 1 0.00264 1 1 0.00053 1 1 0.3583 1 1 0.0082 To quantify the efficiency of both methods in these cases, the figure of merit (FOM), a quantity to characterize the performance of a method, is introduced [ 23 23 ]: <math display="block" id="M23"><mrow><mi>F</mi><mi>O</mi><mi>M</mi><mo>=</mo><mtext> </mtext><mn>1</mn><mo>/</mo><mrow><mo>(</mo><mrow><msup><mi>σ</mi><mn>2</mn></msup><mi>T</mi></mrow><mo>)</mo></mrow></mrow></math> where T is the computational cost and <math id="M24"><mrow><msup><mi>σ</mi><mn>2</mn></msup></mrow></math> is average of the squared deviation between the MC and the analytical values. The greater the q value is, the better the method. The figure of merit of both methods for the four cases are listed in Table 3 Table 3 , It can be seen that the performance is improved by the biasing transition rate method. 1 1 Case No. 1 1 MC 1 1 BTR 1 1 1 1 1 9.50 1 1 464 1 1 2 1 1 0.70 1 1 1.89 1 1 3 1 1 1.20 1 1 41.4 1 1 4 1 1 3.35 1 1 40.7 4 Conclusion 1 A biasing transition rate method for safety assessment of a complex system by MC simulation is proposed. The estimator of this method is stated and variance of the MC estimator is decreased. Four cases are used to benchmark this method. It is an effective method for modeling system failure, being more efficient at collecting evidence of rare events than the direct MC method in the same number of MC cycles. This method may be applied to the rare event simulation of a complex system to save computational cost. When applying an MC simulation with a deterministic code in a safety assessment of a nuclear system, its performance advantage may be more prominent. The performance can be further improved by coupling this method with another efficient method, the Monte Carlo dynamic event tree (MCDET). References IAEA , Applications of probabilistic safety assessment (PSA) for nuclear power plants. IAEA-TECDOC-1200 . 1 - 96 ( 2001 ) F. Di Maio , M. Vagnoli , and E. Zio , Transient identification by clustering based on integrated deterministic and probabilistic safety analysis outcomes . Ann. Nucl. Energy . 87 , 217 - 227 ( 2016 ). doi: 10.1016/j.anucene.2015.09.007 http://doi.org/10.1016/j.anucene.2015.09.007 Transient identification by clustering based on integrated deterministic and probabilistic safety analysis outcomes Y.C. Wu , Development of reliability and probabilistic safety assessment program RiskA . Ann. Nucl. Energy . 83 , 316 - 321 ( 2015 ). doi: 10.1016/j.anucene.2015.03.020 http://doi.org/10.1016/j.anucene.2015.03.020 Development of reliability and probabilistic safety assessment program RiskA Y.C. Wu , P. Liu , L.Q. Hu et al ., Development of an integrated probabilistic safety assessment program . Chin. J. Nucl. Sci. Eng . 27 , 270 - 276 ( 2007 ) Development of an integrated probabilistic safety assessment program E. Zio , The Monte Carlo simulation method for system reliability and risk analysis . Springer ( 2013 ) The Monte Carlo simulation method for system reliability and risk analysis L.Q. Hu , and Y.C. Wu , Probabilistic safety assessment of the dual-cooled waste transmutation blanket for the FDS-I . Fusion Eng. Des . 81 , 1403 - 1407 ( 2006 ). doi: 10.1016/j.fusengdes.2005.08.091 http://doi.org/10.1016/j.fusengdes.2005.08.091 Probabilistic safety assessment of the dual-cooled waste transmutation blanket for the FDS-I Y.C. Wu , Y.Q. Bai , Y. Song et al ., Development strategy and conceptual design of China Lead-based Research Reactor . Ann. Nucl. Energy . 87 , 511 - 516 ( 2016 ). doi: 10.1016/j.anucene.2015.08.015 http://doi.org/10.1016/j.anucene.2015.08.015 Development strategy and conceptual design of China Lead-based Research Reactor Y.C. Wu , J.Q. Jiang , M.H. Wang et al ., A fusion-driven subcritical system concept based on viable technologies . Nucl. Fusion . 51 , 1 - 7 ( 2011 ). doi: 10.1088/0029-5515/51/10/103036 http://doi.org/10.1088/0029-5515/51/10/103036 A fusion-driven subcritical system concept based on viable technologies Y.C. Wu , Conceptual design and testing strategy of a dual functional lithium-lead test blanket module in ITER and EAST . Nucl. Fusion . 47 , 1533 - 1539 ( 2007 ). doi: 10.1088/0029-5515/47/11/015 http://doi.org/10.1088/0029-5515/47/11/015 Conceptual design and testing strategy of a dual functional lithium-lead test blanket module in ITER and EAST A.J. Buslik , Monte Carlo methods for the reliability analysis of Markov systems nuclear safety , in proceedings: of EPRI NP-3912-SR , San Francisco, USA , ( 1985 ) Monte Carlo methods for the reliability analysis of Markov systems nuclear safety G. Bucciarelli , D. Franciotti , M. Marseguerra et al ., Monte Carlo reliability analysis of a safety plant of the Italian Gran Sasso high energy physics national laboratory by means of the MARA code . in Proceedings of ESREL2003 , Maastricht, Netherlands , ( 2003 ) Monte Carlo reliability analysis of a safety plant of the Italian Gran Sasso high energy physics national laboratory by means of the MARA code Y.C. Wu , J. Song , H. Q. Zheng et al ., CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC . Ann. Nucl. Energy . 2 , 161 - 168 ( 2014 ). doi: 10.1016/j.anucene.2014.08.058 http://doi.org/10.1016/j.anucene.2014.08.058 CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC M. Marseguerra , E. Zio , F. Cadini , Biased Monte Carlo unavailability analysis for systems with time-dependent failure rates . Reliab. Eng. Syst. Saf . 76 , 11 - 17 ( 2002 ). DOI: 10.1016/S0951-8320(01)00139-9 http://doi.org/10.1016/S0951-8320(01)00139-9 Biased Monte Carlo unavailability analysis for systems with time-dependent failure rates R.Y. Rubinstein and G. Samorodnitsky , Variance reduction by the use of common and antithetic random variables . J. Stat. Comput. Simul . 22 , 161 - 180 ( 1985 ). doi: 10.1080/00949658508810841 http://doi.org/10.1080/00949658508810841 Variance reduction by the use of common and antithetic random variables H. Kumamoto , K. Tanaka , K. Inoue et al ., Dagger-sampling Monte Carlo for system unavailability evaluation . IEEE Trans. Reliab . 29 , 122 - 125 ( 1980 ). doi: 10.1109/TR.1980.5220749 http://doi.org/10.1109/TR.1980.5220749 Dagger-sampling Monte Carlo for system unavailability evaluation G.W. Imbens and T. Lancaster , Efficient estimation and stratified sampling . J. Econom . 74 , 289 - 318 ( 1996 ). doi: 10.1016/0304-4076(95)01756-9 http://doi.org/10.1016/0304-4076(95)01756-9 Efficient estimation and stratified sampling Y. Choe , E. Byon , N. Chen , Importance sampling for reliability evaluation with stochastic simulation models . Technometrics . 57 , 351 - 361 ( 2015 ). doi: 10.1080/00401706.2014.1001523 http://doi.org/10.1080/00401706.2014.1001523 Importance sampling for reliability evaluation with stochastic simulation models L. Campioni , P. Vestrucci , Monte Carlo importance sampling optimization for system reliability applications . Ann. Nucl. Energy . 31 , 1005 - 1025 ( 2004 ). doi: 10.1016/j.anucene.2004.01.004 http://doi.org/10.1016/j.anucene.2004.01.004 Monte Carlo importance sampling optimization for system reliability applications L. Campioni , R. Scardovelli , P. Vestrucci , Biased Monte Carlo optimization: the basic approach . Reliab. Eng. Syst. Saf . 87 , 387 - 394 ( 2005 ). doi: 10.1016/j.ress.2004.06.008 http://doi.org/10.1016/j.ress.2004.06.008 Biased Monte Carlo optimization: the basic approach L. Campioni , P. Vestrucci , Monte Carlo importance sampling optimization for system reliability applications . Ann. Nucl. Energy . 31 , 1005 - 1025 ( 2004 ). doi: 10.1016/j.anucene.2004.01.004 http://doi.org/10.1016/j.anucene.2004.01.004 Monte Carlo importance sampling optimization for system reliability applications A.O. Lima , C.R. Menyuk , I. T. Lima , Comparison of two biasing Monte Carlo methods for calculating outage probabilities in systems with multisection PMD compensators . IEEE Photonics Technol. Lett . 17 , 2580 - 2582 ( 2005 ). doi: 10.1109/LPT.2005.859183 http://doi.org/10.1109/LPT.2005.859183 Comparison of two biasing Monte Carlo methods for calculating outage probabilities in systems with multisection PMD compensators C.R. Menyuk , Statistical errors in biasing Monte Carlo simulations with applications to Polarization-Mode Dispersion compensators . J. Lightwave Technol . 24 , 4184 - 4196 ( 2006 ). doi: 10.1109/JLT.2006.883131 http://doi.org/10.1109/JLT.2006.883131 Statistical errors in biasing Monte Carlo simulations with applications to Polarization-Mode Dispersion compensators E. Zio , Biasing the transition probabilities in direct Monte Carlo . Reliab. Eng. Syst. Saf . 47 , 59 - 63 ( 1995 ). doi: 10.1016/0951-8320(94)00035-M http://doi.org/10.1016/0951-8320(94)00035-M Biasing the transition probabilities in direct Monte Carlo We thank other members of the FDS Team for their help in this research. 10.1007/s41365-017-0255-2 This work was supported by the Special Projects of International Thermonuclear Experimental Reactor (2015GB116000), the Strategic Priority Research Program of Chinese Academy of Sciences (No. XDA03040000), the Informatizational Special Projects of Chinese Academy of Sciences (No. XXH12504-1-09), the Major / Innovative Program of Development Foundation of Hefei Center for Physical Science and Technology (No. 2014FXCX004). 2016-07-15 2016-09-12 2016-09-15 2017-05-26 2017-07-01 2021-04-25 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2017 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-017-0255-2 http://dx.doi.org/10.1007/s41365-017-0255-2 . 2017 Nuclear Science and Techniques 07 28