2.1. CTMC Models of Highly Reliable Systems

We further elaborate on the description of the highly reliable Markovian systems presented in Section 1. Let there be R types of components in the system with n_i components of type i. Some of the n_i components can act as spares for that type. Let N be the total number of components of all types in the system (i.e.,

). There are many repairman classes in the system, with a finite number of repairmen in each repairman class. Each component type is assigned a repairman class. The component types assigned to a repairman class are divided into priority classes and repaired in a preemptive or nonpreemptive manner. Any service discipline can be used for members of the same priority class.

Let {X(t) : t ≥ 0} be a CTMC model (assume left-continuous with right limits) of a highly reliable Markovian system. For the simplest such system, X(t) = (X₁(t),X₂(t),…,X_R(t)), where X_i(t) is the number of components of type i that are failed. For the general, highly reliable Markovian systems we consider, we have to add some more process descriptors in the state of the system. For example, we might have to add an ordered list of components waiting to be repaired in each repairman class. The following analysis is independent of which definition of X(t) we use.

Transitions in the CTMC correspond either to a component failing (which may cause the instantaneous failure of other components through failure propagation) or a component completing repair. We will refer to the first as a “failure transition” and to the second as a “repair transition.” We label the state with all components up as 1. We will also use 1 to denote the set of states that contains only the single state 1. Assume that X(0) = 1 unless stated otherwise. Let S be the set of states that are accessible from 1. For the purpose of simulation analysis, we will restrict {X(t) : t ≥ 0} to the set S. We partition S into two subsets: S = U ∪ F, where U is the set of up states and F ≠ [empty ] is the set of down states. Of course, 1 ∈ U and the state where all components are failed is in F. We will need the following assumptions:

Let Q = {q(x,y), x,y ∈ S} be the rate matrix (also called the infinitesimal generator matrix) of the CTMC. In this matrix, we arrange the states in the order of increasing number of components failed. Thus, the first state is one in which no components are failed (i.e., 1). This is followed by states in which exactly one component is failed and so on. Let q(x) = −q(x,x) denote the total rate out of state x and h(x) = 1/q(x) be the mean holding time in that state. From (A.2) and (A.3), q(x) > 0 for all x ∈ S. Let Φ denote the probability measure on the sample paths of this CTMC. For any E ⊂ S, let T_E = inf{t > 0 : X(t−) ∉ E,X(t) ∈ E}. Of particular interest are T₁ and T_F.

Let {Y_n : n ≥ 0} denote the embedded discrete-time Markov chain (DTMC) of {X(t) : t ≥ 0}; that is, {Y_n : n ≥ 0} has transition matrix P = {P(x,y) : x,y ∈ S}, where P(x,x) = 0 and P(x,y) = q(x,y)/q(x) for x ≠ y. One can simulate a CTMC by progressively generating the next state of the embedded DTMC and generating the random holding time in that state. For any E ∈ S, let τ_E = inf{n ≥ 1 : Y_n ∈ E}. Of particular interest are τ₁ and τ_F.

As is well known, positive-recurrent CTMCs are regenerative (see, e.g., Crane and Iglehart [8] for the definition) in nature. In the following study, we will be considering system regenerations that occur when the system enters state 1. In any regenerative cycle, let Z be the random variable denoting the holding time in state 1 and let W be the remaining time until either state 1 is reached (again) or the system fails. In mathematical terms, W = min(T₁,T_F) − Z, where T₁ and T_F are measured from the start of the cycle. We will let W (resp., W) be the random variable having the distribution of W given that a system failure (resp., no system failure) occurs in a regenerative cycle. Another quantity that will be used in the following analysis will be γ ≡ P{T_F < T₁} = P(τ_F < τ₁) (i.e., the probability of system failure occurring in a cycle). Assumption (A.3) ensures that γ is not trivially equal to zero. Assumptions (A.2) and (A.4) ensure that 1 − γ is not equal to zero. For any regenerative-cycle-based random variable, say W, we will use W_i to denote its value in the ith regenerative cycle.

2.2. Modeling Highly Reliable Components

In mathematical models of highly reliable Markovian systems, the failure rate of any component, say component i, is assumed to be of the form λ_iε^r_i, where ε is a small positive parameter called the rarity parameter; r_i and λ_i are positive constants (Shahabuddin [36]; see also Gertsbakh [10]). Let r₀ = min{r₁,…,r_N} > 0. Because the repair rates are large compared to failure rates, the repair rate of component i is represented by a constant μ_i > 0. The failure-propagation probabilities are either assumed to be constants or of the same form as the failure rates (i.e., constants multiplied by ε raised to positive powers). Once we introduce this ε-parameterization, then the performance measures become functions of ε. However, for simplicity, we do not specify this dependence in the notation. For example, we continue to use γ for what should ideally be denoted by γ(ε). In highly reliable Markovian systems, the aim is to study the performance measures and the variance of their estimators for small ε.

This particular ε-parameterization guarantees that if Q(x,y) > 0 (resp., P(x,y) > 0) for some ε = ε₀ > 0, then Q(x,y) > 0 (resp., P(x,y) > 0) for all 0 < ε < ε₀. Given the arrangement of states in Q that we mentioned earlier, it can easily be seen that all the elements above the diagonal are O(ε) and all elements below the diagonal are O(1). (A function f (ε) is defined to be O(ε^d) (resp., O(ε^d)), d ≥ 0, if there exists a constant K such that | f (ε)| ≤ Kε^d (resp., ≥ Kε^d) for all sufficiently small ε. A function is said to be Ω(ε^d), d ≥ 0, if it is both O(ε^d) and O(ε^d); that is, a function is Ω(ε^d) if it is exactly of order ε^d. A function f (ε) is defined to be o(ε^d), d ≥ 0, if | f (ε)|/ε^d → 0 as ε → 0.) This structure played an important role in the steady-state-simulation analysis of highly reliable Markovian systems in Shahabuddin [36].

One key property of highly reliable Markovian systems that will be used later is the special structure of the regenerative cycle. Because state 1 has no repair transitions, q(1) is Ω(ε^r₀) since it is the sum of only component failure rates. An important consequence of this is that E(Z) is Ω(ε^−r₀). By Assumption (A.2) and the fact that μ_i's are positive constants, we have that q(x) = Ω(1) for all x ∈ S − 1. Finally, as mentioned earlier, repair rates are very large compared to failure rates. Using these three facts, we see that most of the regenerative cycles consist of a long time interval in which the system is in state 1, after which there is a single failure transition (which might correspond to more than one component failing because of failure propagation), followed by a short time interval in which the failed components are repaired (and thus the cycle completes). This is the intuition behind the fact shown in Shahabuddin [36] that E(min(T_F,T₁)) = E(Z + W) is of the same order as E(Z) (i.e., Ω(ε^−r₀)). Using similar methods, we can show that the expected regenerative-cycle time, E(T₁), is also Ω(ε^−r₀) and that E(W), E(W), and E(W²) are O(1). For completeness, the formal proofs are given in the Appendix. It was also shown in [36] that γ = Ω(ε^r), where r is some nonnegative number depending on the structure of the system. Using the fact that the mean time to system failure (MTTF) for regenerative systems can be expressed as E(min(T_F,T₁))/γ (see, e.g., Goyal et al. [15]), we see that the MTTF is Ω(ε^−(r₀+r)).

2.3. Importance Sampling for Highly Reliable Systems

Consider the problem of estimating the probability, α, of a rare event {X ∈ A}, where X is a random sample path with probability measure Φ and A is a rare set. The “naive” way of estimation is to generate samples of X, say X₁,X₂,…,X_n, from Φ and then form the sample mean

where I(·) is the indicator function of the event inside the bracket. The variance of this estimator is (α − α²)/n ≈ α/n, as α is small. The expected half-width of the 100(1 − δ)% confidence interval, HW, is proportional to

. Thus, the relative error (RE), which is defined by RE = HW/α, is proportional to

. For fixed n, as α → 0, RE → ∞. This is the problem with naive estimation of the probability of rare events.

In importance sampling, we use a change of measure Φ_new, such that for each sample path x ∈ A, Φ_new(x) > 0 if Φ(x) > 0. Then, we can express α as

where L(x) = dΦ(x)/dΦ_new(x) when dΦ_new(x) > 0 and 0 otherwise. The function L is the Radon–Nikodym derivative; it is also called the likelihood ratio. The subscript in the second expectation operator denotes the new probability measure under which the expectation is taken. Equation (1) forms the basis of the technique of importance sampling. The last expectation term suggests that we use the probability measure Φ_new(·) and generate the samples (X_i,L(X_i)). Then, a new unbiased estimator is given by

and its variance is

The main problem in importance sampling is to find an easily implementable Φ_new so that E(I(X ∈ A)L(X)) << α, and so the variance of the new estimator is significantly less than that of the naive one.

The most common importance-sampling technique used for highly reliable systems is failure biasing, proposed by Lewis and Böhm [25]. The basic idea behind failure biasing is to make component-failure transitions in the embedded DTMC happen with a probability that is much higher than in the original system. In state 1, there are no repair transitions. Therefore, we do not need to failure bias in this state. However, in states that have both failure and repair transitions, the total probability of repair transitions ≈ 1 and the total probability of failure transitions ≈ 0. In such states, the total probability of failure transitions is increased to θ, where θ is some constant (i.e., it is independent of ε) between 0 and 1 that is significantly larger than the failure transition probabilities (in practice, θ is typically taken to be 0.5). Therefore, the total probability of repair transitions is decreased to 1 − θ. In a version of failure biasing called balanced failure biasing (Goyal et al. [15] and Shahabuddin [36]), the probability of each failure transition (that had positive probabilities in the original system), conditioned on the event that the transition that occurs is a failure transition, is made the same (this is also done in state 1). The probability of each repair transition, conditioned on the event that the transition that occurs is a repair transition, is left unchanged. Let P′ be the new transition-probability matrix corresponding to balanced failure biasing. Note that P′(x,y) > 0 if and only if P(x,y) > 0.

We will now review the order of magnitude results for the variance associated with importance sampling in the estimation of γ. To obtain a sample of I(T_F < T₁), we need only simulate one regenerative cycle of the CTMC. Thus, γ is called a “regenerative-cycle-based measure.” Other examples of regenerative-cycle-based measures are E(W), E(WI(T_F < T₁)), E(WI(T_F > T₁)), E(T₁), and E(min{T_F,T₁}). Let Φ′′ be a change of measure on the sample paths of the CTMC where we use P′ until time min(T_F,T₁) and then P after that. For any stopping time τ of the embedded DTMC, define

Let τ_min = min(τ_F,τ₁). Then, in the same spirit as Eq. (1), we can write

The variance of the importance sampling estimator, σ_γ²(Φ′′) ≡ Var_Φ′′(I(T_F < T₁)L_τmin), is given by σ_γ²(Φ′′) = E_Φ′′(I(T_F < T₁)L_τmin²) − γ². Let σ_γ²(Φ) be the variance of the naive estimator.

The following theorem was proved in Shahabuddin [36] under the following strengthened version of Assumption (A.4):

Theorem 1: Both γ and σ_γ²(Φ) are Ω(ε^r), where r is a positive constant depending on the structure of the system, the r_i's, and the failure propagation probabilities. Also, E_Φ′′(I(T_F < T₁)L_τmin²) is Ω(ε^2r) and, so, σ_γ²(Φ′′) is O(ε^2r).

This theorem implies that we get a bounded RE in the estimation of γ. One can prove similar bounded RE results for other “rare” regenerative-cycle-based measures like E(WI(T_F < T₁)). Moreover, observe that Assumption (A.4′) is not at all restrictive. If it does not hold, then γ is Ω(1) and we do not have to use importance sampling to estimate it. Unless otherwise stated, in the following sections we will assume that Assumptions (A.1), (A.2), (A.3), and (A.4′) hold.

For transient measures like the unreliability and expected interval unavailability, in addition to failure biasing, we have to use another technique called forcing (Lewis and Böhm [25]). If the time horizon is orders of magnitude less than E(Z), then the system will fail very rarely in [0,t], even though we failure bias. To avoid this, the time of the first event is sampled from the distribution of Z conditioned on the fact that it is less than t. Because Z is exponentially distributed with rate q(1), the time of the first event that we use in the simulation is sampled from the distribution function given by

where 0 ≤ s ≤ t.

3.1. The Small-Time-Horizon Case

We can express the unreliability as U(t) = E(I(A)), where A is the event {T_F < t}. Let Φ′ be the new measure on the sample paths of {X(t) : t ≥ 0}, in which in each replication we use P′ until the system fails and P from then on. For the unreliability, for each sample, we only need to simulate the CTMC until either the system fails or the time horizon is exceeded. Let Φ_Forcing′ be the measure corresponding to both balanced failure biasing and forcing. The variances, without and with balanced failure biasing, are denoted by σ_U(t)²(Φ) and σ_U(t)²(Φ′), respectively. The variance with balanced failure biasing and forcing is denoted by σ_U(t)²(Φ_Forcing′). The following theorem provides the orders of magnitude of the unreliability and the variances of its estimators when the time horizon is small.

Theorem 2: Consider the case where t = Ω(ε^r_t) with r_t = 0. Then, both U(t) and σ_U(t)²(Φ) are Ω(ε^r+r₀). Also, σ_U(t)²(Φ′) = O(ε^2r+r₀) and σ_U(t)²(Φ_Forcing′) = O(ε^2(r+r₀)) (r is the same as in the order of magnitude expression for γ in Theorem 1).

From Theorem 2, we get the following corollary:

Corollary 1: The RE using Φ_Forcing′ (corresponding to a fixed 100(1 − δ)% level of confidence), for a fixed number n of replications, remains bounded as ε → 0.

Lemmas 1 and 2 give bounds on U(t) that are used to prove the order of magnitude result for U(t) in Theorem 2. The result for σ_U(t)²(Φ) is a consequence of the fact that σ_U(t)²(Φ) = U(t) − U²(t). To simplify notation, let q ≡ q(1). Let K be the random variable denoting the number of times the Markov chain is in state 1 (including the first time) before hitting a state in F. Clearly, K has a geometric distribution with parameter (i.e., success probability) γ.

Lemma 1: Let U(t) = 1 − e^−qtγ. Then, for all ε and t,

Proof: Let W_tot be the total amount of time the CTMC spends in states other than 1 before hitting F. Let f_Wtot|k(·) be the density of W_tot given that K = k and let Erlang(q,k) denote an Erlang random variable with rate q and shape parameter k. For a real-valued a, let (a)⁺ = max(a,0). Then,

Lemma 2: Let q_min = min(q(x) : x ∈ U − 1). Then, for all ε and t,

where k is some constant.

Proof: From Shahabuddin [36], it can be shown that there exists a sequence of transitions that start from state 1 and reach a state in F without reentering state 1 such that the product of their probabilities is Ω(ε^r). All other paths have probability of order O(ε^r). In highly reliable systems terminology (see, e.g., Gertsbakh [10]), this is one of the “most likely paths” to system failure. Let (x₀,x₁,x₂,…x_k) be the sequence of states visited in one such path, where x₀ = 1 and x_k ∈ F and x_i ∈ U − 1 for 1 ≤ i < k. Then,

where the inequality in the fifth line above follows from the fact that the region of integration in the fifth line is contained in that of the previous line. █

In the same spirit as Eq. (1), we can write U(t) = E(I(T_F < t)) = E_Φ′(I(T_F < t)L_τF) because we terminate a replication once the failed state is reached. Hence, now we can use Φ′ and obtain samples of (I(T_F < t),L_τF). Consequently, if we use balanced failure biasing, the new variance is given by σ_U(t)²(Φ′) ≡ Var_Φ′(I(T_F < t)L_τF) = E_Φ′(I(T_F < t)L_τF²) − U²(t).

Lemma 3, which follows, gives an upper bound on E_Φ′(I(T_F < t)L_τF²) which will enable us to evaluate the order of magnitude of σ_U(t)²(Φ′). Let D ≡ E_Φ′′(I(T_F < T₁)L_τmin²) = E(I(T_F < T₁)L_τmin) and B ≡ E_Φ′′(I(T_F > T₁)L_τmin²) = E(I(T_F > T₁)L_τmin), where Φ′′ has been defined in Section 2.3. The order of magnitude of D (resp., B − 1) is given in Theorem 1 (resp., Proposition 1).

The key to this result lies in the fact that we can decompose the likelihood ratio L_τF into a product of likelihood ratios over individual cycles. For any sample path with K = k, we can write

where L⁽ⁱ⁾ represents the likelihood ratio over the ith regenerative cycle in the sample path. Given K = k, the L⁽ⁱ⁾ are mutually (conditionally) independent, with L⁽¹⁾,…,L^(k−1) having the distribution of L_τmin given that {τ₁ < τ_F} and L^(k) has the distribution of L_τmin given that {τ₁ > τ_F}.

Lemma 3: For all ε and t,

Proof: In the same spirit as Eq. (1),

Given K = k, (Z₁ + Z₂ + ··· + Z_k) is (conditionally) independent of L_τF. Therefore, we can write

Substituting this in Eq. (4), we obtain

Carrying out the necessary algebra, we get the result of the lemma. █

We will now describe the contribution of forcing to the reduction in variance. Recall that Φ_Forcing′ denotes the new measure on the sample paths of {X(t) : t ≥ 0} when with balanced failure biasing (i.e., Φ′) we also use forcing. Then, we have the following lemma.

Lemma 4: σ_U(t)²(Φ_Forcing′) = (1 − e^−qt)E_Φ′(I(T_F ≤ t)L_τF²) − U²(t).

Proof: Let L_Forcing be the likelihood ratio incurred due to forcing. Note that L_Forcing = 1 − e^−qt. Then,

Proof of Theorem 2: From Lemma 1, we get that U(t) is O(ε^r+r₀). Because q is Ω(ε^r₀), we have that (1 − e^−qt/k) is Ω(ε^r₀). Moreover, since q_min is Ω(1), we have that (1 − e^{−q_min t/k}) is also Ω(1). It then follows that U(t) is O(ε^r+r₀) and we get the first part of the theorem.

From Theorem 1, we see that D = Ω(ε^2r). Using Proposition 1 and the fact that e^x = 1 + x + o(x), we get that e^(B−1)qt − 1 is Ω(ε^r₀). Therefore, by Lemma 3, we get that E_Φ′(I(T_F < t)L_τF²) is O(ε^2r+r₀). Hence, σ_U(t)²(Φ′) = E_Φ′(I(T_F < t)L²) − U²(t) is O(ε^2r+r₀). Then, using Lemma 4, we get that

3.2. The Large-Time-Horizon Case

Consider the following generalization of Theorem 2, which gives the orders of magnitude of the variances of our estimators for large time horizons.

Theorem 3: Consider the case of large time horizons where t = Ω(ε^−r_t) with 0 ≤ r_t ≤ r₀. Then, both U(t) and σ_U(t)²(Φ) are Ω(ε^r+r₀−r_t) and σ_U(t)²(Φ_Forcing′) = O(ε^{2(r+r₀−r_t)}).

The proof of this theorem follows from Lemmas 1–4 and Proposition 1 by substituting t = Ω(ε^r_t) in all of the expressions involving t and using the same method as in the proof of Theorem 2.

We saw in the previous section that for small t, for the bounded RE property to hold, the variance reduction using importance sampling had to be of the same order as the unreliability. From Theorem 3, we see that for the case of large time horizons, we also get a variance reduction that is of the same order as the unreliability as long as r_t ≤ r₀.

Let

We can easily show that for 0 < r_t ≤ r₀, E_Φ′(I(T_F < t)L_τF²)/V(t) → 1 as ε → 0. We conjecture that this holds even for r₀ < r_t < r₀ + r_t. If this conjecture is true, then we have an approximate expression for the variance for the case when ε is small. Then, it is easy to see why simulation using importance sampling becomes very inefficient for the case where r_t > r₀. This is also explained by the results in Glynn [11]; the variance of the likelihood ratio increases (roughly) exponentially fast in the number of transitions of the Markov chain. In experiments in Nicola, Nakayama, Heidelberger, and Goyal [30], it is shown that even for the case where r_t = r₀, one has to tune the value of the failure-biasing parameter θ by trial and error, which can be computationally expensive. Hence, for the case where r_t ≥ r₀, it is best to use the “bounding approach,” as mentioned in Section 1.

Motivated by this, we now investigate bounds on the unreliability. As mentioned in Section 1, because the system is regenerative, the time to failure is roughly a geometric sum of i.i.d. random variables. There is some literature on bounds on the distribution function of such random variables, which, in our case, corresponds to the unreliability. We investigate the bounds given in [4]. For completeness, we first present the bound in its original form. We then adapt it to the reliability model described in this article. Let V_i (generically denoted by V) be i.i.d. nonnegative random variables and let V′ be another nonnegative random variable independent of the V_i's. Let

, where N₀ is a geometric random variable, independent of the V_i's and V′, with probability of “success” p (i.e., P(N₀ = i) = (1 − p)ⁱ⁻¹p for i ≥ 1). It is well known in the literature (e.g., Keilson [24] and Solovyev [40]) that as p → 0, T converges in distribution to an exponentially distributed random variable with mean E(T). Here, E(T) = (1 − p)E(V)/p + E(V′). Theorem 2.2 of Brown [4] gives the following bound:

In our setting, a “success” in the geometric random variable definition corresponds to system failure happening in a regenerative cycle. More specifically, if we define N₀ ≡ K, V_i ≡ Z_i + W_i for 0 ≤ i ≤ K − 1, V′ ≡ Z_K + W_K (recall that W_i is the random variable W_i conditioned on no failure event occurring in that cycle and that W_i is the random variable W_i conditioned on a failure event occurring in that cycle), and p ≡ γ, then T ≡ T_F. An important point to note is that in our setting, even though W_i is not independent of K, the W_i, 0 ≤ i ≤ K − 1, and W_K are independent of K. Also,

Define U_b(t) = 1 − e^−t/E(T_F). Let I ≡ I(T_F < T₁) (resp., I ≡ I(T_F > T₁)) and define

Then, using Eq. (5), we see that upper and lower bounds on U(t) = P(T_F ≤ t) are given by U_b(t) = U_b(t) + U_err,b(t) and U_b(t) = U_b(t) − U_err,b(t), respectively. The bounds are in terms of regenerative-cycle-based measures.

The quality of any upper bound and lower bound combination depends on how close they are to each other. Define the relative error of the bounds (this should not be confused with the relative error, RE, in the simulation context, that was defined earlier), REB, to be the difference between the upper bound and the lower bound divided by twice the measure of interest. (The “twice” is motivated by the fact that if we use the arithmetic mean of the upper and lower bound as an approximation for our measure of interest, then the relative error between the approximation and the actual value is always less than REB.) In this case, the REB is given by REB_b(t) ≡ (U_b(t) − U_b(t))/2U(t) = U_err,b(t)/U(t).

Theorem 4: If r_t > r₀, then U_err,b(t)/U_b(t) → 0 as ε → 0; if r_t = r₀, then U_err,b(t)/U_b(t) is Ω(1); if r_t < r₀, then U_err,b(t)/U_b(t) → ∞ as ε → 0. Hence, REB_b(t) → 0, as ε → 0 if and only if r_t > r₀.

Proof: We use the representation given in Eq. (6). Assumptions (A.2) and (A.4′) imply that 1 − γ is Ω(1). Using the orders of γ, 1 − γ, E(W), E(W), E(W²), and E(Z) (see Sect. 2.2), we get that E(Z + W) = Ω(ε^−r₀), E(Z + W) = Ω(ε^−r₀), and E((Z + W)²) = Ω(ε^−2r₀). Consequently, U_err,b(t) = Ω(ε^r). Moreover, since E(T_F) = Ω(ε^−r−r₀), U_b(t) = Ω(ε^r+r₀−r_t) for r_t < r + r₀ and U_b(t) = Ω(1) for r_t ≥ r + r₀. The result for all three cases follows from these facts. The last part of the theorem follows from the fact that

Therefore, REB_b(t) → 0 if and only if U_err,b(t)/U_b(t) → 0. █

We also tried using the bounds in Kalashnikov [23], but there seems to be some error in the bounds. Moreover, as mentioned earlier, the bounds do not make use of the special structure of the regenerative cycles; that is, they do not make use of the fact that min{T₁,T_F} = Z + W, where Z is exponentially distributed and E(Z) >> E(W). The following theorem provides bounds that make use of this special structure.

Theorem 5:

(i) Let U(t) = 1 − e^−qγt. Define

and let U(t) ≡ U(t) − U_err(t), where

Remark 1: These bounds are in terms of regenerative-cycle-based measures.

Remark 2: These bounds converge for a much wider range of the time horizon as compared to the range in Theorem 4.