The estimation of the failure probability is an important and hard problem, arising in the optimal control of unreliable systems. The required performance measure is usually not analytically available. Thus, simulation becomes the most adequate technique to evaluate a rare event failure probability. We develop a variance reduction technique to estimate a small failure probability in a special unreliable system with repairs. A few variants of the importance sampling method are adapted to estimate the quantity of interest.
Currently, research in the eld of the technical system reliability is associated with rather complex models of degradation processes, for which it is di cult to obtain analytical or heuristic solutions. An important feature of assessing the parameters of the degradation process or the probability of system failure is the varying rate of degradation and failures with a natural or predetermined threshold level, which is usually a random variable [8,15]. The key point in the analysis of complex shock models is the transition to simulation methods and naive Monte Carlo method is still popular in a large number of modern works (see [5,6,11,16,17]), despite its widely known ine ciency for the study of highly reliable systems. Speed-up simulation methods and variance reduction techniques provide a more e cient way of investigating models associated with rare events than the naive Monte Carlo method.
In case of applying variance reduction techniques for the analysis of the degradation models the non-standard formulation of the rare events estimation problem to be considered, when the level of system failure is a random variable. In particular, earlier in [3] the variance reduction technique based on AsmussenKroese method [1] has been extended to estimate the failure probability that a random sum exceeds a random threshold. This approach provides a decreasing of the relative error for the heavy-tailed degradation stages.
In this work, we describe the di erent strategies of applying the importance sampling (IS) method for the considered problem of estimation of the rare-event failure probabilities related to a degradation process containing a few steps, in which a maintenance repair is used to prevent a failure. The rest of the paper is organized as follows. In the second section, we recall the basic de nitions of the degradation process. Then, in the next third section, the problem of rare-event estimation is formulated. Section 4 provides the description of the importance sampling (IS) estimators adapted to the considered failure probabilities. Finally, numerical results are presented in Section 5. 2
Consider the degradation process X := fX(t); t state F and state space 0g with a complete failure
E = f0; 1; : : : ; L; : : : ; M; : : : ; K; F g; that represents degradation stages of the system (see Fig. 1). The states (stages) L, K, M play a speci c role and assumed to be xed. With the help of this model, the dynamics of the state of the anti-corrosion coating was previously studied, and in some cases, an analytical solution for the characteristics of the process is known (see [4]). Nevertheless, this model allows us to develop and test a number of methods for variance reduction and further extend these techniques to more general degradation processes with independent increments. In what follows, we will assume that the two-level policy for managing preventive maintenance and restoration of the system is xed and the values of (K; L) are known, where K stands for preventive repair stage, after which the process returns to stage L. Consider the possible states of the system, illustrated in Fig. 1. The process starts in state X(0) = 0 and then passes K 1 intermediate degradation stages.
Two types of failures are possible in the system: gradual failure, which corresponds to the next stage transition, and instantaneous failure, when the process enters the nal state F . Tracking of instantaneous system failures is a key feature for highly reliable systems modeling.
Let Ti be a random length of the ith stage of the degradation, with a given probability density function (pdf) fTi . (Note that in general, distribution of Ti may depend on the stage number i.)
An instantaneous failure becomes possible starting from stage M , when the system does not have time to go to stage K and, accordingly, the time SM;K =
K 1 X Ti i=M of transition from stage M to stage K exceeds the time to failure V , where V is a random value (r.v.) with given distribution function FV . As a result, the process X jumps to the state F . Following the transition to the state F , the system again is ready for work from the state 0 after a random repair time UF with a known distribution. If, on the contrary, the process X enters the state K, then repair is performed during random time UK;L, and after that, the process returns to stage L.
Degradation process X is regenerative with two types of regeneration cycles. Denote by YF the length of the regeneration cycle with a failure, and let YNF be the length of the cycle without failure. Thus (unconditional) regeneration cycle length Y can be written as
Y = YF IfV SM;Kg + YNF IfSM;K<V g: where If g denotes indicator function and cycle lengths are determined as follows: YF = V + UF + S0;M
YNF = SM;K + UK;L + SL;M ; where r.v. V , UF , S0;M = PM 1
j=0 Tj are independent as well as r.v. SM;K , UK;L, SL;M = PM 1
j=L Tj .
Since for reliability models the greatest risk comes from the instantaneous failure, then the main target is to nd the probability of instantaneous failure within the regeneration cycle, that is l = P(SM;K
V ) = ESM;K E[I(SM;K
V )jSM;K ] = ESM;K [FV (SM;K )]; where it was taken into account independence of the r. v.'s V and SM;K .
Another critically important reliability descriptor of the system quality is a reliability function which is de ned as
R(t) = P[T > tjX(0) = 0]; t 0; where T stands for the lifetime of the system. In our model, the mean lifetime can be written as
E[T ] = E[YNF ](E[N ] 1) + E[V jV
SM;K ]; where N is the number of cycles completed before the total failure happens. It is easy to see that the mean number of such cycles equals E[N ] = 1=`.
For the exponential reliability model, the classical evaluation of the function R(t) based on the solution of the Kolmogorov di erential equations for the state probabilities, thus the results are obtained in terms of the Laplace transform. However, this approach is not acceptable in the rare failure case in general. Nevertheless, for rare failures, the reliability function asymptotic can be used (see [4]):
R(t)
e t E[YN`F ] as ` ! 0; and the mean lifetime becomes (approximately)
E[T ] =
R(t)dt Z 1 0
E[YNF ] : `
If ` is close to zero, then the naive Monte-Carlo fails, thus the variance reduction methods must be applied. 3
In this section, we describe a general problem arising in the estimation of a small probability (the rare event simulation problem). Suppose that > 0 is the so-called rarity parameter, namely ` := P(SM;K
V ) ! 0; as ! 1: The rarity parameter is associated with the parametric set of the distributions V , in particular, one can consider = EV .
Denote by Z an estimator of ` , that means EZ = ` . To estimate `
by Monte Carlo (MC) simulation, one has to generate the independent and
identically distributed (i.i.d) copies of the r. v. Z and calculate the sample
mean
`^ :=
1 XN Z(n):
N n=1
(
The typical measure of the goodness of the estimator is expressed by the relative error (RE) de ned as follows:
The simple example is indicator function of the target event, i.e.
RE h`^ i := r
Var h`^ i E h`^ i
: ZMC = I(SM;K
V );
(
1 p` N as ` ! 0; where a b means that a=b ! 1.
From the equation above it is clear that the RE of the naive MC estimator diverges for small values of the target probability and the sample size required to get a suitable RE is inversely proportional to p` .
There are many rare event simulation techniques [12,9] which aim at
modifying the estimator (
Importance sampling is a popular method for variance reduction aimed at changing the probability measure, so that the target rare event becomes more likely to occur.
We assume that Ti are i.i.d. with a common pdf fT , then the pdf of the random vector (TM ; :::; TK ) has the following product form f (x) = f (xM ; ::; xK 1) = and g(x) = g(xM ; :::; xK 1) be some proposal density function, then the target probability
Z ` = h (x) f (x) g(x) g(x)dx = Eg h (X) f (X) g(X)
; h (x) = FV Q( ) = E[e Ti ] where X1; :::; XN is the sample of the random vector with the density g. 4.1
In this subsection, we describe the classical algorithm based on the importance
sampling with a exponential change of measure. This approach can be applied
to a restricted class of problems when the moment generating function
exists, that is for the light-tailed stages Ti. The classical approach consists in
selecting an importance sampling distribution from a family of pdf's
(
The target probability can be expressed as
` = E [FV (SM;K )W (SM;K )] ;
where E is the expectation under pdf (
W (SM;K ) = exp n
SM;K + (K
M ) ln Q( )o: So, given the two sequences of the samples
fT (i) = (T M(i); :::; T K(i) 1); i = 1; :::; N g; fV (i); i = 1; :::; N g;
where elements of the random vector T (i) are independent realizations of the r.v.
with the twisted pdf fT de ned in (
N `^IS := 1 X FV (SM(i);K )W (SM(i);K );
N i=1 where SM(i);K are independent copies of SM;K .
To select the twisting parameter , adaptive IS is used. Let consider the second moment of the IS estimator:
E I(SM;K
V )W 2(SM;K ) (Q( ))2(K M)E[e 2 V ]
It is reasonable to choose which minimizes given above upper bound. For example, when Ti exp( ) and V = const, one can obtain
V K + M = :
V
In general, such a one-dimensional optimization problem can be solved numerically. 4.2
For the positive function h consider the following distribution
h (x)f (x) g (x) = R h (x)f (x)dx = h (x)f (x) `
It is quite straightforward to show that g is the optimal density which provides the zero variance of the estimator.
But in practice, it is not implementable because it requires knowledge of the target quantity `. Nevertheless, given above expression provides an insight into the general form of a good proposal distribution. Following [2] rewrite the intractable optimal proposal: where 1 ` 1 ` g (x) = h (x)f (x) =
exp( E(x)); E(x) =
log h (x) log f (x): The Laplace method gives a Gaussian approximation of g (x) based on a local expansion around a mode x = argmin E(x):
x In general, the mode can be calculated by numerical methods. Then a Taylor expansion up to second order around the mode gives
E(x)
E(x ) + (x
1
x )rE x=x + 2 (x
x )T H(x
x ):
where rE is the gradient vector and H is the Hessian matrix evaluated at the
mode. This leads to the following Gaussian approximation of the zero-variance
distribution:
g (x) =
e
Then, having the sample X1; :::; XN of the multivariate normal distribution, the
corresponding estimator `bLA is calculated by formula (
g( ) = N ( j x ; H 1): 4.3
The aim is to nd the proposal distribution g close to the desired zero-variance distribution g in the sense of the Kullback-Leibler divergence
D(g ; g) = Eg log g (x) g(x) =
Z g (x) log g (x)dx
Z g (x) log g(x)dx:
Additionally assume that the nominal pdf fT is parameterized by a nite dimensional vector u: fT (x) = fT (x; u). Let the proposal pdf is f ( ; v) for some parameter v (i. e. belongs to the same parametric class). The CE algorithm is based on nding an optimal parameter v : where v = argmin D(g ; f ( ; v))
v = argmax Ew [h (X)W (X; u; w) log f (X; v)] ; v
W (X; u; w) = f (X; u) f (X; w) :
The given above stochastic optimization problem is usually replaced by its stochastic counterpart: v = argmax v 1 N1
X h (Xk)W (Xk; u; w) log f (Xk; v);
N1 k=1
(
The main problem when dealing with the rare-event simulation is that most
values of the h (Xk) are zero. In this case the so-called multi-level CE procedure
[10,13,14] can be applied. This procedure constructs the sequence ( i; vi) of both
rarity and reference parameters, where vi is a solution of the problem (
In this section, we present the simulation analysis of the considered IS estimators.
We have used N = 10000 replications for all experiments. We compare proposed
IS estimators with the naive MC estimator in terms of the relative error (RE)
denoted by the formula (
` RER = jb
`j 100%; ` where ` is the true value (when available). 5.1
We assume that the stages Ti have exponential distribution with a rate parameter . Performance analysis of the IS with exponential twisting estimator is presented in Tables 1, 2. In the rst example, we assume that V is uniformly distributed on the interval ( ; + ), with > > 0, in which case is the rarity parameter. The following numerical parameters are used: K M = 10; = 1; = 5. The second example, where V is exponentially distributed with the rate parameter 1= , is used as a benchmark for which the closed-form expression is available [4] The same set of parameters was used as for the rst experiment.
Numerical results indicate a signi cant decrease of the RE. Moreover, simulation results are very close to the theoretical ones. However, as it was pointed out above, in general (non-exponential) case the sampling from the twisting distribution could be technically di cult.
RE(`bMC) RE(`bIS) RER(`bIS) In this subsection, two strategies of the approximation of zero-variance distribution are compared, namely the Laplace approximation (LA) and the crossentropy (CE) algorithm.
We assume that the stage T has Weibull distribution,
FT (t) = 1 exp( (t=b)a); t 0; where a > 0 is a shape parameter and b > 0 is a scale parameter. In the experiments we put b = 1.
Whereas the implementation of the Laplace approximation is quite
straightforward, the CE procedure requires additional speci cation of the parameters,
namely the reference parameter u of the nominal pdf f ( ; u) and the number
of replications N1 in the stochastic counterpart optimization problem (
u = (uM ; :::; uK 1) = (1; :::; 1):
The CE algorithm attempts to nd the new set of the scale parameters
v = (vM ; :::; vK 1)
which provides the optimum to (
The simulation experiments were performed with the following values of the parameters: K M = 10; = 5; N1 = 100. The numerical results presented in the Table 3 show that both LA and CE algorithms provide a signi cant decrease of the RE. It seems that the LA procedure slightly outperforms the CE for su ciently small probabilities but we also note that the RE is also estimated from generated samples, thus is a r. v. itself. Furthermore, increasing the number of samples N1 leads to further decreasing of the RE. 6
In this paper, we consider a variance reduction technique based on the importance sampling with the di erent strategies of choosing the proposal distribution, to estimate the failure probability of the degradation process. A few numerical examples were presented to demonstrate the properties and accuracy of the estimators. It is assumed in future research, to study the performance of the proposed estimator by both numerical and analytic methods and combine it with the advanced acceleration techniques based on the splitting of the trajectories. As a possible extension one can consider the so-called heterogeneous case when the stages are not identically distributed. 9. Kroese, D.P., Taimre, T., Botev, Z.I.: Handbook of Monte Carlo Methods. John
Wiley & Sons (2011) 10. Kroese, D.P., Rubinstein, R.Y., Glynn, P.W.: Chapter 2 - the cross-entropy method for estimation. In: Rao, C., Govindaraju, V. (eds.) Handbook of Statistics, Handbook of Statistics, vol. 31, pp. 19 { 34. Elsevier (2013). https://doi.org/https://doi.org/10.1016/B978-0-444-53859-8.00002-3 11. Lin, Y.H., Li, Y.F., Zio, E.: A comparison between monte carlo simulation and nite-volume scheme for reliability assessment of multi-state physics systems. Reliability Engineering & System Safety 174, 1{11 (2018) 12. Ross, S.M.: Simulation. Elsevier (2006) 13. Rubinstein, R.Y., Kroese, D.P.: The Cross-Entropy Method: A Uni ed Approach to Combinatorial Optimization, Monte Carlo Simulation and Machine Learning.
Springer-Verlag, New York (2004) 14. Rubinstein, R.Y., Kroese, D.P.: Simulation and the Monte Carlo method. John
Wiley & Sons, Inc., New Jersey (2017)
15. Shahraki, A., Yadav, O.P., Liao, H.: A review on degradation modelling and its
engineering applications. International Journal of Performability Engineering 13(