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

NUCLEAR ENERGY SCIENCE AND ENGINEERING

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

Xiao-Lei Pan，

Jia-Qun Wang，

Run Yuan，

Fang Wang，

Han-Qing Lin，

Li-Qin Hu，

Jin Wang

Nuclear Science and Techniques

Vol.28, No.7

Article number 91

Published in print 01 Jul 2017

Available online 26 May 2017

DOI：10.1007/s41365-017-0255-2

33400

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 simulationProbabilistic safety assessmentBiasing transition rate method

1 Introduction

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-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, 5, 6]. In the PSA analysis of some advanced nuclear systems, such as the China Lead-based Research Reactor (CLEAR) [7], the fusion-driven subcritical system [8] and the International Thermonuclear Experimental Reactor (ITER) [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, 10-12], but it is not efficient at simulating rare events in a complex system [13]. Although sampling algorithms of antithetic variable sampling [14], dagger sampling [15], and stratified sampling [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-20] and multicanonical Monte Carlo (MMC) [21, 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

2.1. Assumptions

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

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 $T_{t}$ 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 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.

Fig. 1

Schematics of biasing transition rate simulation for a component

The flowchart of this method is shown in 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

Fig. 2

The flowchart of the biasing transition rate method.

p_{i j} = {\begin{matrix} \sum_{s = 1}^{q i j} θ (s), q_{i j} \geq 1 \\ 0, q_{i j} = 0 \end{matrix}

(1)

where q_ij is the number of branch points in the j^th sequence generated in i^th trial, and

θ (s) = {\begin{matrix} \frac{1}{(n + 1)}, the s - th transition caused by the real component \\ \frac{n}{(n + 1)}, otherwise \end{matrix}

(2)

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 (λ =λ₀, λ₀ is the design transition rate, a constant)

2) A Weibull distribution for the component transition (λ =αλ₀t^α−1)

3) A linear aging model based on an exponential distribution (λ = λ₀+ k t, k is the aging factor)

4) A linear aging model based on a Weibull distribution (λ = αλ₀ t^α−1+ k t, $α$ 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

F = 1 - e^{- λ t}

(3)

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

F = 1 - e^{- (n + 1) λ t}

(4)

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.

Parameter settings for the biasing transition rate method

Cases	Real component	Virtual component	Integration
1	λ₀	nλ₀	(n+1) λ₀
2	αλ₀ t^α−1	nαλ₀ t^α−1	(n+1)αλ₀ t^α−1
3	λ₀+ k t	n (λ₀+ k t)	(n+1) (λ₀+ k t)
4	αλ₀ t^α−1+ k t	n (αλ₀ t^α−1+ k t)	(n+1) (αλ₀ t^α−1+ k t)

*For Cases 1–4, the contribution rate θ is 1/(n+1) and n /(n+1) for the real and virtual components, respectively.

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

Let $z_{i j}$ 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 $Q_{b}$ based on the biasing transition rate method,

Q_{b} = \frac{1}{N} \sum_{i = 1}^{N} \sum_{j = 1}^{M_{i}} ψ (z_{i j})

(5)

where, N is the number of MC cycles; $M_{i}$ is the sequence number of the i^th trial; and $ψ$ is a discrete function expressing the unexpected event, which can be expressed as Eq.(6) for the biasing transition rate method,

ψ (z_{i j}) = {\begin{matrix} p_{i j}, if the unexpected event occurs \\ 0, otherwise \end{matrix}

(6)

The variation of Q_b can be expressed by

V a r {Q_{b}} = V a r (\frac{1}{N} \sum_{i = 1}^{N} \sum_{j = 1}^{M_{i}} ψ (z_{i}_{j})) = \frac{1}{N^{2}} (\sum_{i = 1}^{N} V a r (\sum_{j = 1}^{M_{i}} ψ (z_{i}_{j})))

(7)

The failure probability of the system Q can be estimated by the direct MC estimator $Q_{d}$ ,

Q_{d} = \frac{1}{N} \sum_{i = 1}^{N} ψ (c_{i})

(8)

where $c_{i}$ donates the sequence in the i^th trial. The variance of the estimator $V a r {Q_{d}}$ is

V a r {Q_{d}} = V a r \sum_{i = 1}^{N} ψ (c_{i}) = \frac{1}{N^{2}} (\sum_{i = 1}^{N} V a r (ψ (c_{i})))

(9)

For the direct MC method, the discrete function expressing the unexpected event can be expressed as

ψ (c_{i}) = {\begin{matrix} 1, if the unexpected event occurs \\ 0, otherwise \end{matrix}

(10)

Since E $(\sum_{j = 1}^{M_{i}} ψ (z_{i j})) = E (ψ (c_{i}))$ , and 0< P_u <1/2, we have

V a r (\sum_{j = 1}^{M_{i}} ψ (z_{i j})) \leq V a r (ψ (c_{i}))

(11)

and

V a r {Q_{b}} = \frac{1}{N^{2}} (\sum_{i = 1}^{N} V a r (\sum_{j = 1}^{M_{i}} ψ (z_{i j}))) \leq \frac{1}{N^{2}} (\sum_{i = 1}^{N} V a r (ψ (c_{i}))) = V a r {Q_{d}}

(12)

This proves that the biasing transition rate method can decrease the variance of the MC estimator.

3 Benchmark cases

3.1. Description of the system

To illustrate performance of the method, a system consisting of three components, as shown in Fig. 3, is used to benchmark the biasing transition rate method. The following four cases are considered:

Fig. 3

Diagram of the benchmark system.

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 λ₁ = λ₂ = 10⁻⁴ and λ₃ = 10⁻⁵; the biasing factor is n = 9; and number of MC cycles is N = 1000. For Case 2, the parameter factors are α₁ = α₂ = α₃ =1.1. For Case 3, the aging factors are k₁ = k₂ =10⁻⁸ and k₃ =10⁻⁹. For Case 4, design repair rates are u₁ = u₂ =10⁻³ and u₃ =0.

3.2 Results and discussion

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 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.

Fig. 4

Cumulative probabilities for system failure at T_m=1000 h, N =1000, and different distributions of the component failure time, repair time, and linear aging.

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.

RMSD and efficiency of collecting evidence of failure events for the direct MC and biasing transition models*

Case No.	Computational cost /s		Number of evidence		RMSD		Expected time to ‘system failure’ /s
	DMC	BTR	DMC	BTR	DMC	BTR	DMC	BTR
1	5.6865	7.7354	28	1213	0.00647	0.00034	0.2031	0.0064
2	6.1644	12.6543	44	4640	0.00482	0.00204	0.1401	0.0027
3	4.6761	7.0110	11	1315	0.00421	0.00057	0.4251	0.0053
4	4.2993	8.5845	12	1042	0.00264	0.00053	0.3583	0.0082

* DMC, direct MC; BTR, biasing transition rate model.

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]:

F O M = 1 / (σ^{2} T)

(13)

where T is the computational cost and $σ^{2}$ 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, It can be seen that the performance is improved by the biasing transition rate method.

FOMs (×10⁴) for MC and BTR methods

Case No.	MC	BTR
1	9.50	464
2	0.70	1.89
3	1.20	41.4
4	3.35	40.7

4 Conclusion

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

[1]

IAEA, Applications of probabilistic safety assessment (PSA) for nuclear power plants. IAEA-TECDOC-1200. 1-96 (2001)

[2]

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