CC BY
95
ESTIMATION OF RELIABILITY IN INTERFERENCE MODELS USING MONTE CARLO SIMULATION
A.N. Patowary, J. Hazarika and G. L. Sriwastav
Department of Statistics, Dibrugarh University, India e-mail: a patowary@yahoo.com, iitenhazarika@rediffmail.com, wastav2@yahoo.com
Abstract
This paper presents estimation of reliability R = P(X > Y) of a system, for the cases when its strength (X) and stress (Y) follow exponential, normal or gamma distributions, using Monte Carlo simulation (MCS). First the parameters of strength and / stress are estimated and substituting them in the reliability expressions, in different cases, the estimates of reliability are obtained.
A
Normal distribution is fitted to various sets of estimated reliability R, generated by MCS. The goodness of fit is tested using Kolmogorov-Smirnov one sample test.
Keywords: Stress-Strength; Monte-Carlo Simulation; Kolmogorov-Smirnov one sample test.
1. Introduction
In interference theory of reliability, reliability and other reliability characteristics of a system can be expressed as some functions of the parameters of the distributions of the random variables (r.v.'s), strength (X) and stress (Y) associated with the functioning of the system. We estimate these parameters and substitute these values in the expressions for reliability and other characteristics to get their estimates. The estimates of parameters used here are maximum likelihood estimators and as such from the invariance property of MLE's, the corresponding estimators of reliability are also MLE's. In absence of hard data the numerical values of the estimators can be obtained from simulation. There exists extensive literature for estimation of reliability analytically for single component systems e.g. Mazumder [12], Church and Harris [4] etc. But the reliability expressions for multi-component systems are not simple enough to facilitate analytical estimation of reliability and its other characteristics. Also due to lack of hard data, one way out is simulation, in particular Monte Carlo simulation.
With simulation technique it is possible to estimate reliability (or probability of failure) or other reliability characteristics without going into the analytical techniques. The availability of personal computer and software makes the process comparatively simple. In fact, to evaluate the accuracy of the sophisticated analytical techniques or to verify a new technique, simulation is routinely used to independently evaluate the underlying probability distributions.
1.1 Monte Carlo Simulation:
The Monte Carlo Simulation method is an artificial sampling method which may be used for solving complicated problems in analytic formulation and for simulating purely statistical problems. In the simplest form of simulation, each r.v. in a problem is sampled several times to
represent its real distribution. Each realization of r.v.'s in the problem produces a set of numbers that indicates one realization of the problem itself. Solving the problem deterministically for each realization is known as a simulation cycle, a trial or a run. Using many simulation cycles we get the overall probabilistic characteristics of the problem, particularly when the number of cycles N is sufficiently large. Simulation, using a computer, is an inexpensive way (compared to laboratory testing) to study the uncertainty in a problem.
The primary components of a Monte Carlo Simulation include the followings:
(i) Probability distribution function (or probability density function): The physical (or mathematical) system must be described by a set of probability distribution functions.
(ii) Random number generator: A source of random numbers uniformly distributed on the unit interval must be available.
(iii) Sampling rule: A prescription of sampling from specified distribution function, assuming the availability of random numbers on the unit interval.
(iv) Scoring (or tallying): The outcome must be accumulated into overall tallies or scores for the quantities of input.
In this paper except exponential distribution we have not used uniformly distributed random numbers, rather obtained random numbers following particular distribution directly from MATLAB.
In Section 2, we have estimated reliability of an n-standby system (n=1, 2, 3), through Monte Carlo Simulation technique. Simulation is performed for exponential stress-strength, normal stress-strength and gamma stress-strength. In Section-3 we have considered fitting of normal distribution to estimated reliability in each case, for different true values of the parameters. The goodness of fit is tested by K-S one sample test (Seigel [18]). Since we have taken a small sample, 20, only, when using x2test, the number of classes becomes too few, due to pooling. We have considered the fitting of normal distribution to check whether normal approximation is good enough for a small of 20.
Some literatures on the topic which we have come across are:
Kamat and Riley [8] presented MCS for a complex system for time to failure (TTF) models.
Some of the others studies of reliability estimation using MCS for TTF models includes Pulido et.al. [15], Goel [5], Hong and Lind [6], Landis et.al. [10], Tunak et.al. [23], Naess et.al. [13], Wu et.al. [24] etc.
Stancampiano [21] applied simulation to interference models. Manders et.al. [11], Aldrisi [1], Stumpf and Schwartz [22], Zhang et.al. [25] have simulated stress-strength. Paul and Borhanuddin [14], Rezaei et.al. [17] estimated reliability of stress-strength model, using MCS. Ahmad et.al.[2] obtain Bayes estimates of P (Y< X) using MCS. Borhanuddin et.al. [3] estimated reliability for multicomponent system using MCS. Rao et al [14] compared reliability estimates for multicomponent systems evaluated by different methods such as method of moments, modified ML method and Best Linear Unbiased Estimator through MCS technique.
Kakati and Sriwastav [7] and Sriwastav [20], used simple simulation by taking random exponential numbers to represent stress-strength. They considered very small samples. From these samples they first estimate the parameters and substituting these in the expressions of reliability they get estimated reliability.
2 Reliability Estimation through Monte Carlo Simulation:
Let us consider an n-standby system. Let X1, X2,.. .,Xn be the strengths of the n components in the system arranged in the order of activation. Let Y1, Y2,...,Yn be the stresses faced, respectively, by 1st, 2nd,...,nth component, when they are activated; Xi's and Yi's are all independent. For a detailed description of such a system one may refer (Sriwastav and Kakati, [19]).
The reliability Rn of an n-standby system for a single impact of stress is given by,
Rn = R(1) + R(2) + ... + R(n), (2.1)
where R(i) is the increment in the system reliability due to the ith component, defined as R (i )= P [Xx < Y,,X2 < Y2,..., X ; _! < Yi _ 1,Xi > Yi ] (2.2)
Here, we have assumed that all the components are having the same strength distributions and are working under the same environment (stress), i.e. all Xi's and Yi's are i.i.d. with probability density functions (pdf s) f(x) and g(y), respectively.
In this section, we use MCS to estimate reliability. The programs are developed in MATLAB, separately for exponential, normal and gamma. First a set of 5000 values of the particular r.v. viz. (exponential, normal or gamma) are generated for a particular value of the parameter(s). Using these values an estimate of the parameters involved is obtained. Substituting this estimate(s) in the expression of reliability we get an estimate of the reliability. This process is repeated j times to give j estimates of the parameter(s) and subsequently j estimates of reliability. The whole process is repeated for different true values of the parameters; for a particular true value of the parameter(s) j is the sample size.
2.1 (a) Exponential Stress-Strength:
Let us assume that the component's strength follows exponential distribution with mean X and the stress on it follows exponential distribution with mean strength unity, without loss of generality. Then the marginal reliability expression due to the n th component is, (ibid)
R(n) = R(n) = -\i-Vl (23)
1+ X^ 1+ X J
So, R1 = (2.4)
1+ X
R2 = R1 + (1 - R1) R1, (2.5)
R3 = R1 + (1 - R1) R1+ (1 - R1)2 R1. (2.6)
For MCS, let U be the uniform r.v. over (0, 1). Then by following inverse transformation we can generate exponential random variable with mean X as:
Let U = F(x) = 1 - exp (-x / X ) ^ X = - X log (1 - U) Now if U is uniform over (0, 1), (1 - U) is also uniform over (0, 1). So
X = - X log (U) (2.7)
From uniform r.v. U we can generate exponential r.v. with parameter X using the above transformation (2.7). We generate 5000 of U. Then from (2.7), for each U, - log U gives a value of the exponential r.v. X with mean unity. Thus we get 5000 values of X. Multiplying each of these 5000 values of X by X(= 0.5,2,3) we get 5000 values of exponential r.v. (say X1) with mean X. The mean of these 5000 values give an estimate of X for a particular true value of X. Substituting these estimates in (2.4), (2.5) and (2.6) we get an estimate of R1, R2 and R3, respectively. For each true value of X the whole process is repeated j times there by giving j estimates of X and R's for a particular true value ofX . Here, we have taken j = 20.
(b) Normal Stress-Strength: In case of normal stress-strength let X~ n(j,,o2 ) and stress Y ~N (0, 1) by without loss of generality. The reliability expressions are (ibid)
R (n) = [1 - O
^ 1 + n2 J
] • - 1 n J • (28)
R'=n!+? J, (29)
and R2 and R3 are given by (2.5) and (2.6), respectively.
2
For generating normal random numbers for particular true values of ^ and a we first
generate standardized normal numbers (Z) by the MATLAB. Next by the following transformation
2
we generate normal random numbers particular true values of ^ and a
X = + Z a where Z ~ N (0, 1).
A A
We estimate p and a2 from the sample of X of size 5,000. Substituting these estimates in Ri, R2 and R3, given above, we get estimates of the system reliabilities R's. This process is repeated j times for a particular set of set of true values of ^ and a . We have taken (^ , a) = (-1, .5), (0, .5), (1, .5), (2, .5), (-1,1), (0, 1), (1, 1), (2, 1), (-1, 2), (0, 2), (1, 2), (2, 2) and j = 20.
(c) Gamma Stress-Strength: We assume that the strength X ~r(1,m) and stress Y~r(1,k). Then if m and/ or k is an integer (ibid)
R(n) =
( m-1 r(m + k-i -1) ^
n-1
~0 rk(m - i -1)!2
m-1
r(m + k - i -1)
1 _V-1 Y . 3L ' (2.10)
£0 rk(m - i - 1)!2m+k-1-1
so, R1=z r(m+k-i-1> (2.11)
£0 rk(m - i - 1)!2m+k-1-1 v '
Here also R2 and R3 are given by (2.5) and (2.6). If m and k are not necessarily integers (Kapur and Lamberson, [9])
R(n) = [1 - Beta(m,k)Beta1/2(m,k)]n-1Beta(m,k)Beta1/2(m,k) (2.12)
where Beta (. , .) is Beta function and Beta1/2(. , .) is incomplete beta function, with parameters m and k.
So, R1 = Beta(m,k)Beta1/2(m,k) (2.13)
and then R2 and R3 are given by (2.5) and (2.6).
Here also, by MATLAB, we directly generate gamma random numbers X for strength population and Y for stress population, of sizes 5000 each with different true values of the parameters m and k. Substituting these estimates in the above expressions of R's we get a reliability estimate of the systems. This process is repeated j times for a particular set of true values of m and k. Then for each set of true values of m and k the above process is repeated. We have taken (m, k) = (1,1), (1,2), (2,1), (2,2) and j = 20.
3. Fitting of Normal Distribution to Systems Reliability:
From different expressions of reliability in Section-2 we have obtained the estimates of reliability substituting the estimated values of the parameters in respective cases. To these estimated values of reliability for different cases we have fitted normal distribution and tested the goodness of fit by one sample K-S test. The tabulated value of D for sample size 20 at 5% level of significance is 0.294 (see Seigel [18]).
Let us first consider the case of exponential stress-strength. For each X(= 0.5,2,3) j =20, the
AAA A
values of R1,R2 and R3 are obtained by substituting the corresponding values of X obtained in Section-2, in the expressions (2.4), (2.5) and (2.6). Then for each value of true X, mean and s.d. of
A A A
R1 ,R2 and R3 are calculated. In each case normal distribution is fitted and the goodness of fit is tested by K-S test. The values are tabulated in Table-3.1. True values of R1, R2, and R3 are also given in the same table for comparison.
For K-S test, if calculated value is of D < 0.294, the fit is good. From Table- 3.1, columns 5,
A A A
9, 13 it is clear that normal distribution gives good fit to the values of R1, R2 and R3.
Table -3.1: Exponential Stress-Strength
True True Ri Mean A Ri SD A Ri D for A Ri True R2 Mean A R2 SD A R2 D for A R2 True R3 Mean A R3 SD A R3 D for A R3
.5 .333 .333 .003 .090 .555 .556 .004 .067 .704 .704 .004 .116
1 .500 .500 .003 .072 .750 .750 .003 .072 .875 .875 .002 .075
2 .667 .667 .003 .076 .889 .889 .002 .063 .926 .963 .001 .063
3 .750 .750 .003 .071 .938 .948 .001 .088 .953 .984 .000 .023
N.B.: The entry .000 in the SD column indicates that the SD is very small. This is the situation for all the tables.
Next let us consider the case of normal stress-strength. The above procedure is repeated for different set of (a ) and their corresponding estimated values from Section-2 are used in expressions (2.9), (2.5) and (2.6). The results are tabulated in Table- 3.2. From values of D (see Seigel [18]) in column 6, 10, 14 we see that normal distribution gives good fits to the distributions
of Ri,R2 and R3.
Table -3.2: Normal Stress-Strength
True True True Mean SD D True Mean SD D True Mean SD D
a R1 A A for R2 A A for R3 A A for
R1 R1 A R2 R2 A R3 R3 A
R1 R2 R3
.5 -1 .212 .186 .001 .073 .379 .337 .002 .117 .510 .460 .003 .119
0 .500 .500 .003 .800 .750 .750 .003 .080 .875 .875 .002 .091
1 .788 .814 .002 .083 .955 .966 .000 .073 .991 .994 .000 .149
2 .945 .963 .001 .037 .997 .999 .000 .021 .999 .999 .000 .058
1 -1 .308 .240 .003 .141 .522 .424 .004 .115 .669 .561 .004 .205
0 .500 .501 .005 .064 .750 .751 .005 .064 .875 .876 .003 .077
1 .692 .760 .004 .047 .905 .942 .002 .043 .971 .986 .000 .026
2 .841 .922 .002 .123 .975 .994 .000 .105 .996 .999 .000 .129
2 -1 .421 .327 .004 .077 .665 .547 .006 .044 .805 .695 .006 .088
0 .500 .501 .005 .061 .750 .751 .005 .039 .875 .876 .004 .055
1 .579 .672 .004 .107 .823 .892 .003 .078 .926 .965 .001 .129
2 .655 .813 .005 .075 .899 .965 .002 .082 .977 .993 .001 .097
Here we would like to point out that for ^ = 0 and a = 0.5, the fit is not good.
Similarly for gamma stress-strength, for different sets of true values of stress-strength parameters (m, k) and taking their corresponding estimates from Section-2 and using this in expressions (2.11), (2.5) and (2.6) we obtain estimates of Ri, R2, and R3 in different situations and calculate D statistics in each case. All these values are tabulated in Table- 3.3. Comparing the values of D in columns 6, 10 and 14 with the tabulated values (ibid) we see that normal distribution gives good fit to reliabilities of systems for gamma stress-strength also.
Table-3.3: Gamma Stress-Strength (m and/ or k are Integer)
True True True Mean SD D True Mean SD D True Mean SD D
m k Ri A A for R2 A A for R3 A A for
Ri Ri A R 2 R 2 A R3 R3 A
R1 R2 R3
1 1 .500 .500 .007 .088 .750 .750 .007 .063 .875 .875 .005 .085
2 .250 .250 .003 .118 .438 .438 .005 .126 .579 .579 .006 .123
2 1 .750 .745 .012 .103 .938 .935 .006 .087 .984 .983 .002 .119
2 .500 .500 .004 .092 .750 .750 .004 .092 .875 .875 .003 .102
When neither m nor k is an integer then the corresponding estimates of R1, R2, and R3 are
A A
obtained by substituting the values of mand k from Section-2 in the expressions (2.13) etc. and the corresponding values are tabulated in Table-3.4.
Table- 3.4: Gamma Stress-Strength (m and K are not necessarily Integer)
True True True Mean SD D True Mean SD D True Mean SD D
m k R1 A A for R2 A A for R3 A A for
R1 R1 A R2 R2 A R3 R3 A
R1 R2 R3
1 1 .500 .499 .007 .090 .750 .749 .007 .090 .875 .874 .006 .076
2 .250 .250 .005 .092 .438 .438 .007 .005 .579 .579 .008 .065
2 1 .750 .750 .005 .069 .938 .938 .002 .060 .984 .984 .001 .071
2 .500 .500 .006 .095 .750 .750 .006 .070 .875 .875 .004 .096
Conclusion: In this paper, we have estimated the reliability through MCS. We have seen that
normal distribution is fitted to the data sets of estimated reliability obtained by MCS. Once we
know the distribution it is easy for us to obtain the other characteristics of the reliability data sets.
Acknowledgment: This paper was supported by UGC (India) under Major Research Project (Sc.).
References:
[1] Aldrisi, M. M., (1987): 'A Simulation Approach for Computing System Reliability', Microelectronics Reliability, Vol. 27, Issue: 3, pp. 463-467.
[2] Ahmad, K.E., Fakhry, M.E. and Jaheen, Z.F., (1997): 'Empirical Bayes Estimation of P(X>Y) and characterizations of Burr-type X model', Journal of Statistical Planning and Inference, Vol. 64, Issue: 2. pp. 297-308.
[3] Borhan Uddin, Md., Pandey, M., Ferdous, J. and Bhuiyan, M.R., (1993): 'Estimation of Reliability in a Multicomponent Stress-Strength Model', Microelectronics and Reliability, Vol.33, Issue: 13, pp.2043-2046.
[4] Church, J.D. and Harris, V., (1970): 'The Estimation of Reliability from Stress-Strength Relationship', Technometrics, Vol.12, No.1, pp.49-54.
[5] Goel, L., (2000): 'Monte Carlo Simulation Based Reliability Studies of a Distribution Test System', Electric Power System Research, Vol. 54, Issue: 1, pp.55-65.
[6] Hong, H.P. and Lind, N.C., (1996): 'Approximate Reliability Analysis using Normal Polynomial and Simulation Results', Vol.18, Issue:4, pp.329-339.
[7] Kakati, M.C. and Sriwastav, G.L., (1981): 'Accelerated Life-Testing under Stress-Strength Model', IAPQR Trans. Vol.11, No.1-2, pp.51-56.
[8] Kamat, S.J. and Riley, M.W., (1975): 'Determination of Reliability Using Event-Based Monte Carlo Simulation', IEEE Trans. on Reliab., Vol. R-24, No.1, pp.73-75.
[9] Kapur, K.C. and Lamberson, L.R., (1977): 'Reliability in Engineering Design', John Wiley and Sons, New York.
[10] Landis, C.M., Beyerlein, I.J. and McMeeking, R.M., (2000): 'Micromechanical Simulation of the Failure of Fiber Reinforced Composites', Journal of the Mechanics and Physics of Solids, Vol.48, Issue:3, pp.621-648.
[11] Manders, W.P., Bader, M.G. and Chou, T.W., (1982): 'Monte Carlo Simulation of the Strength of Composite Fiber Bundles', Fiber Science and Technology, Vol. 17, Issue:3, pp.183 -204.
[12] Mazumder, M.,(1970): 'Some Estimates of Reliability Using Interference Theory', NRLQ, Vol.17, pp.159-165.
[13] Naess, A., Leira, B.J. and Batsevych, O., (2009): 'System Reliability Analysis by Enhanced Monte Carlo Simulation', Structural Safety, Vol.31, Issue: 5, pp.349 - 355.
[14] Paul, R.M. and Borhar Uddin, MD., (1997): 'Estimation of Reliability of Stress-Strength Model with Non-identical Component Strengths', Microelectronics and Reliability, Vol. 37, Issue: 6, pp.923-927.
[15] Pulido, J.E., Jacobs, T.L. and Patters de Lima, E.C., (1992): 'Structural Reliability Using Monte- Carlo Simulation with Variance Reduction Techniques on Elastic-Plastic Structures', Computers and Structures, Vol. 43, Issue: 3, pp.419-430.
[16] Rao, G.S. and Kantam, RRL, (2010): 'Estimation of Reliability in Multicomponent Stress-Strength Model: Log Logistic Distribution' Eletron. J. App. Stat. Annal., Vol.3, Issue: 2 pp. 75-84.
[17] Rezaei, S., Tahmasbi, R. and Mahmoodi, M., (2010): 'Estimation of P [ Y < X ] for Generalized Pareto Distribution', Journal of Statistical Planning and Inference, Vol. 140, Issue: 2, pp. 480-494.
[18] Seigel, S. (1956): 'Non-Parametric Statistics for Behavioural Sciences', McGraw Hill Book. Co. Ltd., N.Y.
[19] Sriwastav, G.L. and Kakati, M.C., (1981): 'A Stress-Strength Model with Redundancy', IAPQR Trans. Vol.6, No.1, pp.21-27.
[20] Sriwastav, G.L., (1988): 'Reliability Estimation in Stress-Strength Model with Ancillary Information', Assam Statistical Review., Vol. 2, No.1-2, pp. 16-20.
[21] Stancampiano, P.A., (1978): 'Simulation Approaches to Probabilistic Structural Design at Component Level', Clear Engineering, Vol. 50, Issue:1, pp.123-129.
[22] Stumpf, H. and Schwartz, P., (1993): 'A Monte-Carlo Simulation of the Stress-Rupture of Seven-FiberMicrocomposites', Vol.49, Issue:3, pp.251-263.
[23] Tunak, M., Linka, A. and Volf, P., (2009): 'Load-Sharing and Monte Carlo Models of Defects in a Bundle fibres', Composites Science and Technology, Vol.69, Issue: 9, pp. 1417-1421.
[24] Wu, L., Yi, L., Xue-qin, G. and Deng-hai, W., (2009): 'Monte-Carlo Simulation for the Reliability Analysis of Multi-Status Network System Based on Breadth First Search', Information and Computing Science, pp.280-283.
[25] Zhang, H., Chandrangsu, T. and Ramussen, K.J.R., (2010): 'Probabilistic Study of the Strength of Steel Scaffold System', Structural Safety, Vol.32, Issue: 6, pp.393-401.
