Reliability Eng~ineering 10 ~1985) 233 243
The Monte-Carlo Method Without Sorting for Uncertainty Propagation Analysis in PRA
Soon Heung Chang, Joo Young Park and Myung Ki Kim Nuclear Engineering Department, Korea Advanced Institute of Science and Technology, PO Box 150, Chongryang, Seoul, South Korea (Received: 21 April 1984)
ABSTRACT The Monte-Carlo method has been widely used Jor the analysis of an uncertainty propagation in Jault tree analysis. The conventional Monte-Carlo method is time-consuming because of the sorting process. In this stud),, the Monte-Carlo method without sorting is developed by using the segmentation between the sujficient upper bound and the sufficient lower bound of the evaluated top event jrequencies into meshes and interpolation within the mesh. It is" shown that the results obtainedJrom the new method are almost equivalent to those from the conventional method, and that the new method reduces greatly the computational running time.
1.
INTRODUCTION
Fault-tree analysis is one of the most important procedures in the safety study of the nuclear power plant. The top event in a fault tree has uncertainty because of the imperfection of a fault tree construction and the uncertainties of the basic event probabilities.1 Many methods have been developed to obtain the uncertainty of the top event probability 233 Reliability Engineering 0143-8174/85/$03-30 © Elsevier Applied Science Publishers Ltd, England, 1985. Printed in Great Britain
Soon HeungChang, Joo Young Park, Myung Ki Kim
234
Specify parameters of component distribution
rI
Generate random number for each component characteristic distribution
Combine random component values according to system function to yield value of system characteristic distribution
NO
System characteristic function of component characteristics
Sampling complete'? YES Sort distribution values
Calculate cumulative ft. and density function
Fig. 1. Flow chart of conventional Monte-Carlo method.
Monte-Carlo method without sorting in PRA
235
propagated from the uncertainties of the basic event probabilities. Wellknown methods in these analyses can be divided into three groups, as follows. (1) Methods to obtain the cumulative distribution function of the top event using the discretizations of basic event distributions. (2) Methods to obtain the top event distribution from the moments of the top event probabilities which are obtained using the moments of the basic event probabilities. (3) Methods using the Monte-Carlo technique and sorting. Each method has characteristics of its own. In the first group are the histogram method and D P D method. 2 --4 The histogram method provides the desired accuracy in an uncertainty analysis at the expense of high CPU-time requirements. The D P D method may also require comparatively high CPU-time to produce the high accuracy when there are many basic events. The moment method and the method of using fitted response surfaces belong to the second group. 5'6 In the moment method, an empirical distribution is obtained from the first and second moments of the top event. This method has been known to give comparatively reasonable results. The method of using fitted response surfaces has the problem of accuracy and requires high C PU-time in the case of many basic events. However, this method has the advantage of getting the information on which basic event is more sensitive to the uncertainty of the top event probabilities. Methods in the third group are the most powerful that can be used in any case if the basic event distribution is once assumed. The general procedures in these methods are shown in Fig. 1. However, the sorting which this group used is a very time-consuming process. Based on the review of several uncertainty propagation methods, a new method combining the three groups is developed. By using the new method, the computer code M O M O D is developed. The output and the run-time of M O M O D is compared with those of SAMPLE, which is a computer code using the conventional Monte-Carlo method and was used in WASH-1400.1 2.
METHOD DESCRIPTION
The new method developed can be described as the following six steps. First, generate the randomly-sampled probability for each basic event
236
Soon Heung Chang, Joo Young Park, Myung Ki Kim
based on the characterized distribution for the Monte-Carlo simulations. The new method described in this study can treat the basic events which are characterized by any combination of the following probability distributions: (1) (2) (3) (4) (5) (6) (7) (8)
normal, log-normal, log-uniform, uniform, exponential, gamma, Weibull and discrete distribution defined by user.
Random number generation from each basic event distribution was done by the subroutine packages in IMSL. v So this part of M O M O D can be modified according to the computer used and the number of subroutines supplied. Second, from the given relationship between the top event and basic events as the sum of minimal cut sets and sampled probabilities of the basic events, obtain the probability of the top event Qr. The first and second steps are repeated for the sufficient Monte-Carlo runs. The above steps are typical procedures which the conventional MonteCarlo method takes before the sorting process. In other words, after the second step, the new method developed takes a different procedure from that of the conventional Monte-Carlo method, sorting. Third, the mean rn and the variance V of the probability of the top event Qr are used in order to obtain the sufficient confidence bounds in Qr. Assuming that the top event has an ideal log-normal distribution for this purpose, the distribution parameters p and a and the sufficient confidence bounds, AAA and ZZZ, can be obtained from the following equations: a2 = log ( ~ 2 + 1)
(1)
0-2
/.t = log (m) + -2-
(2)
AAA = exp (/~ - 3.0a)
(3)
ZZZ = exp (/~ + 3.0a)
(4)
Monte-Carlo method without sorting in PRA
237
The values AAA and ZZZ, which are 99.73 7o confidence bounds in an ideal log-normal distribution, are used as the sufficient confidence bounds in the new method. Fourth, the distance between AAA and ZZZ is divided into the meshes. There are two choices in the segmentation. One method is a linear scale segmentation with equal distances, and the other is a log-scale segmentation. Monte-Carlo simulations for obtaining the probabilities of the top event
Obtain the mean m and the variance V of the probabilities of the top event
Obtain the sufficient confidence bounds AAA, ZZZ using the values of m and V
Segment the distance between AAA and ZZZ into meshes
Determine the number of probabilities of the top event in each mesh
I Obtain the cumulative distribution function of the top event and plot Fig. 2.
Flow chart of the method developed.
Soon Heung Chang, Joo Young Park, Myung Ki Kim
238
Fifth, after segmentation the mesh that each obtained top event probability belongs to is determined, and the number of top event probabilities obtained in each mesh is counted. Finally, from the information obtained from the above step, an approximate cumulative distribution function of the top event can be obtained. The method of obtaining the cumulative distribution function of the top event is interpolation. After all these procedures, the relative frequency and the cumulative distribution function of the top event is plotted. The procedures in M O M O D which are described above are shown schematically in Fig. 2.
3.
S A M P L E C A L C U L A T I O N S A N D RESULTS
The code developed, M O M O D , was tested in four cases. In the first case, the basic events were characterized by the log-normal distribution, which is most widely used in nuclear power plant failure data. 1 This case dealt with the fault tree in Fig. 3 and the failure data in Table 1. The cumulative distribution functions of the top event obtained from S A M P L E and M O M O D , and the computer running time, are shown in Table 2. In the second case, with the same fault-tree used in the above case, the case when the failure data were characterized by the normal distributions in Table 3 was considered. The results are shown in Table 4. T System Failure
+ Fig. 3.
+ Fault tree.
TABLE 1 Failure Data
Fault
Distribution
X(I) X(2) X(3) X(4) X(5) X(6)
log-normal log-normal log-normal log-normal log-normal log-normal
Median 1.0 3.0 1-0 3.0 1.0 3.0
x x x x x x
Error factor
10- 3 10 2 10 2 10- 2 10 : 10 3
3 3 3 3 3 3
C o m m o n Mode Xc~~= )((3). X(5) QT= X(l) +
(X(2) + X(3)). (X(4) + X(5)) + X(6) +,/X(3)5. Xi5).
TABLE 2 Results of Test Case 1
Code
SAMPLE
SAMPLE
MOMOD (linear)+
MOMOD (linear)
MOMOD (log)t+
MOMOD (log)
imax*
1200
3 600
1 200
3 600
1 200
3 600
100
100
100
isegt
100 _
_
o
Percentile l 2.5 5 10 20 30 40 50 60 70 80 90 95 97.5 99 C PU-time (s) * t + tt
3.190E-3 3.868E-3 4-255E-3 4.954E-3 5.794E-3 6.602E-3 7.338E-3 8.034E-3 8.912E-3 1-003E-2 1.169E-2 1.385E-2 1.655E-2 1.847E-2 2-150E-2
3.154E-3 3.683E-3 4.154E-3 4-778E-3 5.704E-3 6.475E-3 7.192E-3 8.023E-3 8.936E-3 1.002E-2 1.149E-2 1.395E-2 1.638E-2 1.847E-2 2.171E-2
3-212E-3 3.632E-3 4.043E-3 4.707E-3 5.610E-3 6.391E-3 7.150E-3 7.864E-3 8-771E-3 9.938E-3 1.127E-2 1.369E-2 1.720E-2 2.032E-2 2.455E-2
3.189E-3 3.619E-3 4.085E-3 4.761E-3 5.670E-3 6.450E-3 7.218E-3 8.000E-3 8.966E-3 1.012E-2 1-157E-2 1-408E-2 1.698E-2 1.993E-2 2.402E-2
3.189E-3 3.642E-3 4.004E-3 4.727E-3 5.614E-3 6.385E-3 7.518E-3 7.858E-3 8-765E-3 9-934E-3 1.124E-2 1.373E-2 1-721E-2 2.043E-2 2.440E-2
3.177E-3 3.623E-3 4.086E-3 4.759E-3 5.665E-3 6.455E-3 7.218E-3 8.014E-3 8.960E-3 1.013E-2 1.157E-2 1.398E-2 1-698E-2 1.994E-2 2.390E-2
18.31
98.84
5.49
13.67
5.63
14.14
imax = number of Monte-Carlo runs. iseg = number of meshes. M O M O D (linear) = M O M O D with linear segmentation. M O M O D (log) = M O M O D with log-scale segmentation.
Soon HeungChang, Joo Young Park, M y u n g K i K i m
240
TABLE 3 Failure Data
Basic event
Mean
90 '"o Error spread
X(I)
0-001
0.000 3
X(2)
0.03
0.01
Y(3)
0.01
0.003
X(4)
0.03
0.01
X(5) X(6)
0.01 0.003
0.003 0.001
Qr=Y(l)
+ (X(2) + X ( 3 ) ) . (X(4) + X(5)) + X(6) + .v/X(3)2. X(5).
TABLE 4 Results of Test Case 2
Percentile
SAMPLE
MOMOD (log)
1
4.901E-
5
5.382E - 3
3
4.90lE.-
5.400E - 3
MOMOD (linear)
3
4-899E-
3
5.404E - 3
10
5.646E - 3
5.633E - 3
5-637E -- 3
50
6.605E - 3
6.583E - 3
6.580E - 3
90
7-606E - 3
7-60 4E - 3
7.609E - 3
95
7.942E - 3
7.873E - 3
7.865E - 3
99
8.412E - 3
8.393E - 3
8.402E - 3
17- 83
4.15
4.00
C P U - t i m e (s)
TABLE 5 Results of Test Case 3
Uncertainty method MOMOD
(log)
MOMOD SAMPLE
(linear)
Moment method DPD
5 '!;,
Median
95 i'~,
2'666E - 3 2 . 6 4 2 E -- 3 2.581E- 3
5.222E - 3 5.219E - 3 5.350E- 3
1.167E - 3 1.159E - 3 1.171E- 3
2.67E - 3 2.47E - 3
5.27E - 3 5.30E - 3
t.32E - 3 1' 46E - 3
Q r = X1 + X2. X 3 q- X 4 . X 5 -~- X 6. E r r o r f a c t o r = 3.
241
Monte-Carlo method without sorting in P R A TABLE 6 Results of Test Case 3 Uncertainty method
M O M O D (log) M O M O D (linear) SAMPLE Moment method DPD
5 04,
Median
95 ~~
1.832E - 3 1.042E - 3 1.776E - 3 8.36E - 4 1.87E - 3
8.143E - 3 8.187E - 3 8.530E - 3 6.07E - 3 8.75E - 3
4-303E - 2 4.350E - 2 4.533E - 2 5-36E - 2 4-542E - 2
QT = XI + X2 • )(3 nt-X4- X5 -b X6. Error factor = 10.
I n t h e t h i r d case, u s i n g the p r e v i o u s s e n s i t i v i t y s t u d y d o n e ip. ref. 4, MOMOD
w a s c o m p a r e d w i t h the o t h e r m e t h o d s . T h e r e s u l t s a r e s h o w n
in T a b l e s 5 a n d 6. F i n a l l y , t h e case w h e n t h e b a s i c e v e n t s a r e n o t c h a r a c t e r i z e d b y t h e s a m e k i n d o f d i s t r i b u t i o n w e r e c o n s i d e r e d w i t h the s a m e f a u l t tree as t h e first case a n d t h e f a i l u r e d a t a i n T a b l e 7. T h e r e s u l t s a r e s h o w n in T a b l e 8 a n d c o m p a r e d w i t h t h o s e o f the m o d i f i e d S A M P L E . S A M P L E c a n n o t t r e a t t h e case w h e n t h e b a s i c e v e n t s a r e n o t c h a r a c t e r i z e d b y t h e s a m e k i n d o f d i s t r i b u t i o n , so S A M P L E
w a s m o d i f i e d for this p r o c e d u r e .
TABLE 7 Failure Data Fault
Distribution
X(1) X(2) X(3) X(4) X(5) X(6)
normal normal gamma normal gamma normal
Distribution parameter
mean = mean = A* = mean = A* = mean =
0.001, 90 '!o error 0.03, 90130error 3.0, B* = 0'000 4 0.03, 90"; error 3"0, B* = 0"0004 0.003, 90 ?o error
Common mode XcM = X(3). X(5). * Gamma distribution is given as follows f X a ~ I e -X.B
J(x)=IoBgF~A)
elsewhere.X>0(A>0'B>0)'
spread = 0.000 3 spread = 0 0 t spread = 0-01 spread = 0-001
242
Soon Heung Chang, Joo Young Park, Myung Ki Kim TABLE 8 R e s u l t s o f Test C a s e 4
Percentile
Modified SAMPLE
MO MOD (log)
MO MOD ( linea r)
0.349 3E - 02 0-395 7E - - 0 2 0.421 0E - 02 0.501 7F 02 05924E02 06176E02
0-349 3E - 02 0'395 5E - 02 0"420 6E - 02 0"501 4E - 02 0.593 I E - 0 2 0.6173E -02
0.6654E-
0.6622E-02
+
1 5 10 50 90 95
0"348 0-395 0.420 0.501 0-593 0617
99
0.6623E
C P U - t i m e (s)
3E - 02 8E - 02 4E - 02 8E -- 02 5E-02 IE--02
1533
+02
02
4+25
4.09
In the results of the first case, it was shown that the output of M O M O D is almost equal to that of S A M P L E when the basic events were characterized by the log-normal distributions, and that the computer running time required in M O M O D was much less than that of SAM PLE. This was true especially when the larger Monte-Carlo runs were used. F r o m the results of the second case, it was shown that the output of M O M O D was reasonable when another distribution (such as normal distribution versus log-normal distribution) was employed for the failure data. In the results of the third case, it was shown that M O M O D gave results compatible with those of other methods for both the narrow spread and the large spread. It is also seen from Table 6 that M O M O D with log-scale segmentation gave the more reasonable results for the large spread compared with S A M P L E than did M O M O D with linear scale segmentation. Finally, in the results of the fourth case, it was shown that the output of M O M O D was also sufficiently accurate and effective in the computer running time saving when the basic event was characterized by the combination of the various distributions. 4.
CONCLUDING REMARKS
For the uncertainty propagations of probabilistic risk assessment+ the Monte-Carlo method without sorting is developed by using the
Monte-Carlo method without sorting in PRA
243
segmentation between the sufficient upper bound and the sufficient lower bound of the evaluated top event frequencies into meshes and interpolation within the mesh. The method developed was tested when the basic events were characterized by log-normal distribution, normal distribution and the combination of normal distribution and gamma distribution. The method developed gives results equivalent to those of the conventional Monte-Carlo method and other methods of treating the uncertainty propagation in all of test cases. The new method reduces greatly the computational running time, especially when larger MonteCarlo runs were used.
ACKNOWLEDGEMENT The financial support of ASAN acknowledged.
Foundation
for this research is
REFERENCES 1. US NRC WASH-1400 (NUREG 75/014), October, 1975. 2. Colombo, A. G. and Jaarsma, R. J. A Powerful Numerical Method to Combine Random Variables, 1EEE Trans. Reliab., R-29(2) (1980), pp. 126 9. 3. Ahmed, S., Clark, R. E. and Metcalf, D. R. A Method for Propagating Uncertainty in PRA, Nucl. Technol., 59 (1982), pp. 238-45. 4. Ahmed, S., Metcalf, D. R. and Pergram, J. W. Uncertainty Propagation in PRA: A Comparative Study, Nucl. Engng Design, 68 (1981), pp. 1--3, 5. Apostolakis, G. and Lee, Y. T. Methods for the Estimation of Confidence Bounds for the Top Event Availability of Fault Tree, Nucl. Engng Design, 41 (1977), pp. 411-19. 6. Mazumdar, M. An Approximate Method for Computation of Probability Intervals for the Top Event Probability of Fault Trees, Nucl. Engng Design, 71 (1982), pp.45 50. 7. International Mathematical & Statistical Library Inc., The I MSL Library, 1981, 8. Cox, D. C., Baybutt, P. and Kurth, R. E. Comments on the Uncertainty in Accident Consequences Calculated by Large Code due to Uncertainty in Input, Nuel. Technol., 52 (1981), pp. 439-41. 9. Jackson, P. S., Hockenbury, R. W. and Yeater, M. g. Uncertainty Analysis of System Reliability and Availability Assessment, Nucl. Engng Design, 68 (1981), pp. 5 29.