Methodology for the Assessment of Imprecise Multi-State System Availability

Abstract

Most existing studies of a system’s availability in the presence of epistemic uncertainties assume that the system is binary. In this paper, a new methodology for the estimation of the availability of multi-state systems is developed, taking into consideration epistemic uncertainties. This paper formulates a combined approach, based on continuous Markov chains and interval contraction methods, to address the problem of computing the availability of multi-state systems with imprecise failure and repair rates. The interval constraint propagation method, which we refer to as the forward–backward propagation (FBP) contraction method, allows us to contract the probability intervals, keeping all the values that may be consistent with the set of constraints. This methodology is guaranteed, and several numerical examples of systems with complex architectures are studied.

1. Introduction

In traditional binary-state systems studies, we suppose that system components can only have two states: total functioning or complete failure. However, in many industrial systems, this assumption may not be sufficient. Some systems can have more than two states of levels of performance, ranging from total functioning to total failure. In other words, some systems can have many intermediate states between the two extremum cases. These systems are considered multi-state systems (MSS). The first work related to the theory of MSS was that of Barlow and Wu [1] and El-Neveihi et al. [2]. A complete introduction of MSS reliability theory and their applications are represented in the first book related to the reliability analysis of MSSs [3]. An interesting recent review of the literature was proposed in [4,5]. The main known methods for MSS reliability are the Monte Carlo simulation method [6,7], the Markov approach [8,9], the structure function method [10,11], and the universal generating function (UGF) method [12,13,14].

In the structure function method, the classical Boolean models are extended to multi-valued models. It was the first technique introduced in MSS reliability assessment. However, this method and the Markov method generate all the possible states of the system. The Monte Carlo simulations method can reduce the number of simulations, but it leads to an approximative result dependent on the chosen number of simulations. Consequently, these methods are extremely time consuming when considering large MSSs due to the important number of the system states.

On the other hand, the conventional MSS reliability assessment considers two major assumptions [15]:

• The probabilities related to the states of MSS components should be characterized by probability distributions or measures.
• The performance levels of MSS components should be precisely determined.

However, when studying MSS reliabilities, many types of uncertainties concerning the state probabilities and the performance levels of components [16] have to be modeled. The most known categorization is to divide the uncertainty into two types: epistemic and aleatory uncertainties. The latter (also called irreducible uncertainty) arises from the intrinsic variability of a phenomenon across space, through time, or within a population. The former (also called reducible uncertainty) arises from the incompleteness of data/knowledge [17,18,19]. Note that the distinction between these two types of uncertainty is useful when considering computing the reliability of systems [20]. If the data are sufficient, the classical probabilistic approach is the more adequate theory to be used in reliability and risk assessments [20,21,22].

In this paper, we study the problem of representing epistemic uncertainty in the evaluation of multi-state systems/components availabilities, in contrast with works that have considered systems with multi-state performance and only components with binary states [23,24]). To our knowledge, there are only a few works considering the representation of epistemic uncertainty in the evaluation of multi-state systems/components availabilities:

• In [25], Huang et al. used p-boxes and intervals to represent unknown probabilities of component states when considering UGF functions. In [26], the authors proposed a combined approach based on the universal generating function (UGF) and the Dempster–Shafer evidence theory when considering the epistemic uncertainty, transmission loss, and performance storage simultaneously in the availability/reliability assessment of multi-state phased mission systems. Another study considering the same methodology (the combined UGF and Dempster–Shafer theories) was proposed in [27] when considering common cause failures, and was compared to the interval UGF methods in [28]. In [29], the traditional UGF approach for multi-state system (MSS) availability and reliability assessment was extended by using fuzzy random variables and by considering simultaneously the probabilistic convolution and the fuzzy extension principle. However, the UGF method used in these works can only be applied to series/parallel and parallel/series configurations. This differs from our paper in which the proposed approach can be applied to all configurations.
• Li et al. [30] proposed an approach based on the use of interval arithmetics with interval-valued probability masses as a model of unknown probability distributions [31], which are different from belief functions. This approach also uses the UGF method.
• Some works [32,33,34] represented unknown probabilities by fuzzy sets, and propagated uncertainties using the extension principle of Zadeh and fuzzy arithmetics (see, for example, [35]). More specifically, the propagation of uncertainties is based on the use of the extension principle of Zadeh to perform algebraic operations on fuzzy numbers. The extension principle provides a general method for extending crisp mathematical concepts to address fuzzy quantities, such as real algebraic operations on fuzzy numbers, which are a generalization of real numbers in the sense that they do not refer to single values, but rather to a connected set of possible values, where each possible value has its own weight between 0 and 1. In these works, the used fuzzy sets were defined over probabilities and not over component states. They correspond to some categories of hierarchical models [36], i.e., uncertainties were defined over uncertainty models. These approaches are in contrast to our approach defined in the current paper, where the uncertainty representation is characterized by interval values of failure and repair rates [37].
• In [38], the authors studied the reliability analysis of complex MSS with epistemic uncertainty and common cause failures using the Bayesian network (BN) and Dempster–Shafer (DS) theory. The idea here is to add to the Bayesian system network an uncertain state (a virtual state) in order to express the epistemic uncertainty.

To our knowledge, the combination of a Markovian approach with interval contraction methods has not been investigated. Note that a preliminary work concerning our approach was briefly presented in our short conference papers [39,40]. Apart from the novelty of our proposed work, we can note many reasons to justify why such an extension is appealing:

• The theory of imprecise probabilities is completely based on the classical probability theory and thus, it can be considered its generalization. Consequently, conventional MSS models can be considered a special case of the proposed imprecise models.
• Interval contraction methods are guaranteed. For example, when using interval contractions, the result of an interval computation is an interval, which guarantees that the exact value is in the obtained interval.
• The MSS availability results can have a large imprecision (the length of availability intervals). This fact reflects the incompleteness of initial failure/repair data and can help the decision maker to collect more data if it is possible.

Interval analysis was first developed to quantify the error in numerical computations [41]. It is considered a very rigorous theoretical framework for dealing with uncertainties. In order to efficiently characterize the domain of the variables involved in an equation, interval methods use operators called contractors and separators. The interval operators are used to remove some parts of a box that do not satisfy the equation’s constraints. The result of an interval computation is an interval, also known as a pair of numbers or an upper and a lower bound. This pair of numbers guarantees that the exact value is contained in the interval between the two numbers. In our case, we consider that uncertain transition rates are bounded by intervals, $\lambda =[\underline{\lambda },\overline{\lambda }]$ and $\mu =[\underline{\mu },\overline{\mu }]$.

To model the complex structure of the MSS, we propose using a Markovian model and calculating its asymptotic availability by first solving the following system of equations:

$$\Pi .Q=0$$

(1)

where $\Pi =\left[\left[{\pi }_1\right]\dots \left[{\pi }_n\right]\right]$ is the state’s probability vector of the system, and its elements are $\left[{\pi }_i\right]$, the probability interval of the system to be in the ith state. Q is the interval transition matrix. The availability is then equal to the sum of the probabilities related to the functioning states.

In this paper, we propose a new method that can calculate the imprecise availability in a simple way. We propose dealing with imprecise data as intervals and applying interval analysis on imprecise Markovian approaches, making it possible for us to calculate the imprecise availability of a system. To calculate the interval form of the availability, we use “the technique of contractors” as defined in [42]. In particular, we use the contractor of the forward–backward propagation (FBP) in order to contract the intervals of steady probabilities and to efficiently assess the imprecise availability of the studied system.

The paper is organized as follows: In Section 1, we detail the imprecise Markovian model in terms of dependability. In Section 2, we present the original study, which deals with the use of the imprecise Markovian model in terms of dependability, and we show the negative aspects of the method. In Section 3, we briefly describe the interval contraction methods. In Section 4, we present our contribution for the purpose of estimating the imprecise availability of an MSS. This method takes all the constraints into account with the aim of computing the availability of a complex MSS. To understand all the steps of the methodology, numerical examples are presented and compared to other results obtained by different existing methods (IUGF [30], BUGF [43], and exact method). We conclude the paper with a discussion and a summary statement.

2. Imprecise Markovian Models in Dependability Studies

2.1. Basics of Imprecise Markovian Models

In Markovian models, the transition rates (for example failure/repair rates) can be usually influenced by many factors. Thus, the evaluation of the transition rates can be difficult due to several reasons: the use of new components, the lack of data, etc. We are convinced that, in this case, we should incorporate these imprecisions into the model instead of ignoring them. The development of some models integrating imprecise probabilities makes it possible. In [44], the authors introduced the use of imprecise probabilities theory in continuous time Markov chain. According to the procedure used in [44], a discretization of imprecise continuous-time Markov chains used lower and upper transition operators [45] based on some expectation operators.

Thus, in our proposition, we consider that the transition rates are given in the form of intervals. Moreover, the imprecision may exist in the transition matrix or/and the initial state probabilities. Formally, the considered imprecise model is written as follows:

$$\left\{\begin{array}{c}P(X_0=i_0)=[\underline{p_i},\overline{p_i}]\hfill \\ P(X_n=j=i_n|X_{n-1}=i=i_{n-1},\dots ,X_0=i_0)\hfill \\ =P(X_n=j|X_{n-1}=i)\hfill \end{array}\right.$$

(2)

where $x=\left\{1,\dots ,N\right\}$, i and j are two elements of x, and $X_0$, $X_{n-1}$, and $X_n$ are random variables that belong to x and satisfy Equation (2).

We are interested in finding the probability that a system is in a working state at infinity (asymptotic availability). Under the assumptions in Equation (2), the stationarity of the system is thus determined in the form of a vector of intervals $[\Pi ]=\left[\left[{\pi }_0\right]\dots \left[{\pi }_i\right]\dots \left[{\pi }_n\right]\right]$ with $i\in \left\{1,\dots ,n\right\}$, where $\left[{\pi }_i\right]$ is the probability interval of being in a state $e_i$. When considering precise data (i.e., precise transition matrix), we compute the asymptotic availability by solving Equation (1). However, when considering imprecise data, the bounds of the asymptotic availability intervals cannot be computed by simply considering the bounds of the transition matrix [44]. This important remark is discussed in detail in the next subsection. More specifically, some methods were proposed to solve this problem. The “exact method” consists of finding all the possible transition matrices of a system. Then, one has to solve Equation (1) and find a vector of stationary probabilities in order to obtain the vector of intervals that contains all of the possible vectors, and then we simply take the minimum and maximum of this vector. Obviously, this method leads to an accurate result, but its complexity, which can be very important, depends on the system’s size since it takes into account all the possibilities for the different transition matrices. The belief universal generated function (BUGF) [43] and the interval universal generated function (IUGF) [30] proposed two original techniques combining the classical UGF method and intervals. They are, therefore, used in the case of interval-modeled imprecision. These two methods are efficient and give good results. The BUGF is considered to be more efficient than the IUGF, but their use is limited to cases where the system does not have a complex structure (e.g., series-parallel and parallel-series configurations) because the UGF method can be only used for series-parallel and parallel-series configurations.

2.2. Troffaes Method

To the best of our knowledge, only one study has been done on dependability for the availability computing of an MSS based on imprecise Markov models. We now focus on this study and show that part of it is wrong. This paper, published by Troffaes et al. [44], used imprecise continuous-time Markovian models to assess the reliability of power networks. In this study, the authors assumed that the failure and repair rates are bounded by intervals. However, this is not the case for the transition matrix $\left[Q\right]$. Instead, they compute two matrices, $\underline{Q}$, the lower bound of $\left[Q\right]$, and $\overline{Q}$, the upper bound of $\left[Q\right]$. The transition rates $\left[q_{ij}\right]$ depend on the intervals of the failure and repair rates. $\underline{Q}$ is obtained by finding its elements $\underline{q_{ij}}$, where $\underline{q_{ij}}$ are the lower values of $\left[q_{ij}\right]$. $\overline{Q}$ is obtained by finding its elements $\overline{q_{ij}}$, where $\overline{q_{ij}}$ are the higher values of $\left[q_{ij}\right]$. By obtaining two matrices, the problem is turned into two precise cases of Markovian models. By solving $\Pi .Q=0$, the two vectors of $[\Pi ]$ are found, as is the probability vector of the system, to be in two specific states: the lower bound $\underline{\Pi }$ and the upper bound $\overline{\Pi }$. In this way, the interval of the availability of the system is found, and it is formed from two bounds: the lower bound is the sum over all the elements of $\underline{\Pi }$ corresponding to working states, and the upper bound is the sum over all the elements of $\overline{\Pi }$ corresponding to working states.

To more clearly illustrate this methodology, let us analyze the example proposed by the authors. Let us consider a network system composed of two power lines, A and B. The states related to the functioning of this system are represented by the set $E=\left\{AB,A,B,\varnothing \right\}$. The labels of these states represent the non-faulty components (i.e., both A and B are non-faulty in AB, whereas both are faulty in ∅). The common cause failures are modeled by considering the following three events:

• CAB: common cause failure of both A and B.
• BI: independent failure of B.
• AI: independent failure of A.

To simplify, any interval $\left[x\right]$ is written as x. Moreover, $q_1^A$ represents the rate of AI, $q_1^B$ represents the rate of BI, and $q_2$ represents the rate of CAB. Similarly, let $r_B$ be the repair rate of B and $r_A$ be the repair rate of A. For the sake of simplicity, the authors excluded all simultaneous repairs. The transition rate matrix is then

$$\left[Q\right]=\left[\begin{array}{cccc}-q_1^A-q_1^B-q_2& q_1^B& q_1^A& q_2\\ r_B& -q_1^A-q_2-r_B& 0& q_1^A+q_2\\ r_A& 0& -q_1^B-q_2-r_A& q_1^B+q_2\\ 0& r_A& r_B& -r_A-r_B\end{array}\right]$$

(3)

The Markov chain is represented in Figure 1. The interval numerical values of the data [46] are

$$\begin{array}{c}\hfill q_1^A=[0.32,0.37]\\ \hfill q_1^B=[0.32,0,37]\\ \hfill q_2=[0.19,0,24]\\ \hfill r_A=[730,1460]\\ \hfill r_B=[730,1460]\end{array}$$

(4)

expressed per year.

The transition matrix Q is the interval $Q=\left[\underline{Q},\overline{Q}\right]$. The lower and upper bounds are defined by

$$\underline{Q}=\left[\begin{array}{cccc}-\overline{q_1^A}-\overline{q_1^B}-\overline{q_2}& \underline{q_1^B}& \underline{q_1^A}& \underline{q_2}\\ \underline{r_B}& -\overline{q_1^A}-\overline{q_2}-\overline{r_B}& 0& \underline{q_1^A}+\underline{q_2}\\ \underline{r_A}& 0& -\overline{q_1^B}-\overline{q_2}-\overline{r_A}& \underline{q_1^B}+\underline{q_2}\\ 0& \underline{r_A}& \underline{r_B}& -\overline{r_A}-\overline{r_B}\end{array}\right]$$

(5)

$$\overline{Q}=\left[\begin{array}{cccc}-\underline{q_1^A}-\underline{q_1^B}-\underline{q_2}& \overline{q_1^B}& \overline{q_1^A}& \overline{q_2}\\ \overline{r_B}& -\underline{q_1^A}-\underline{q_2}-\underline{r_B}& 0& \overline{q_1^A}+\overline{q_2}\\ \overline{r_A}& 0& -\underline{q_1^B}-\underline{q_2}-\underline{r_A}& \overline{q_1^B}+\overline{q_2}\\ 0& \overline{r_A}& \overline{r_B}& -\underline{r_A}-\underline{r_B}\end{array}\right]$$

(6)

By replacing each term by its value, we obtain

$$\underline{Q}=\left[\begin{array}{cccc}-0.98& 0.32& 0.32& 0.19\\ 730& -1460.61& 0& 0.51\\ 730& 0& -1460.61& 0.51\\ 0& 730& 730& -2920\end{array}\right]$$

(7)

$$\overline{Q}=\left[\begin{array}{cccc}-0.83& 0.37& 0.37& 0.24\\ 1460& -730.51& 0& 0.61\\ 1460& 0& -730.51& 0.61\\ 0& 1460& 1460& -1460\end{array}\right]$$

(8)

The authors evaluated the lower and upper stationary distributions using the above two transition matrices and solving $\underline{\Pi }.\underline{Q}=0$ and $\overline{\Pi }.\overline{Q}=0$, without taking into account the fact that the sum of the steady probabilities is equal to 1 ($\sum {\pi }_i=1$). For the stationary distribution, the authors found

$${\underline{\Pi }}_1=\left[\begin{array}{c}9.985\times {10}^{-1}\\ 2.623\times {10}^{-4}\\ 2.623\times {10}^{-4}\\ 6.513\times {10}^{-5}\end{array}\right]$$

(9)

$${\overline{\Pi }}_1=\left[\begin{array}{c}9.994\times {10}^{-1}\\ 7.252\times {10}^{-4}\\ 7.252\times {10}^{-4}\\ 1.647\times {10}^{-4}\end{array}\right]$$

(10)

such that $[{\underline{\Pi }}_1,{\overline{\Pi }}_1$] is the interval of $\Pi$ obtained using the method presented in the paper. The interval vector $[\Pi ]$ is found by regrouping the two bounds ($\underline{\Pi }$ and $\overline{\Pi }$) together. The authors were able to find the availability of the system from the interval vector $[\Pi ]$. In this example, by calculating the sum over ${\pi }_1,{\pi }_2$ and ${\pi }_3$, we find that the availability is $\left[A_1\right]=[0.9990246,1.00085]$.

However, the authors made a mistake in this study. The lower transition matrix $\underline{Q}$, where all its elements are the lower bounds of $q_{ij}$, is surely not the matrix that gives the lower stationary vector $\underline{\Pi }$. The same remark applies to the upper transition matrix $\overline{Q}$, where all its elements are the upper bounds of $q_{ij}$, which is surely not the matrix that gives the upper stationary vector $\overline{\Pi }$. In other words, there are some combinations of the transition rate intervals that can lead to lower or higher matrices that are different from the ones chosen by the authors to find the solutions (vector of probabilities) and that are outside of the range of values proposed by the authors.

In fact, by applying the exact method, we can demonstrate that their assumption is not always correct. The main idea of the exact method is to obtain all possible combinations of the lower and upper bounds of $q_{ij}$ for each combination for which we have a transition matrix. With the transition matrix, we find the vector $\Pi$ and then the availability of the system. We compare all the values of the availability so that we can choose the matrix $Q_l$ that belongs to $[\underline{Q},\overline{Q}]$, and that gives the lower bound of $[\Pi ]$, which corresponds to the lowest vector among all the $\Pi$s and the lower bound of $\left[A\right]$. $Q_l$ is given as follows:

$$Q_l=\left[\begin{array}{cccc}-0.83& 0.32& 0.32& 0.19\\ 730& -1460.61& 0& 0.51\\ 730& 0& -1460.61& 0.51\\ 0& 730& 730& -2920\end{array}\right]$$

(11)

such that $[{\underline{\Pi }}_2,{\overline{\Pi }}_2$] is the interval of $\Pi$ obtained by the exact method, with

$${\underline{\Pi }}_2=\left[\begin{array}{c}9.981\times {10}^{-1}\\ 5.674\times {10}^{-4}\\ 5.674\times {10}^{-4}\\ 6.978\times {10}^{-4}\end{array}\right]$$

(12)

Another matrix $Q_u$, which also belongs to $[\underline{Q},\overline{Q}]$, gives the upper bound of $[\Pi ]$, which corresponds to the highest vector among all the $\Pi$s and the upper bound of $\left[A\right]$. $Q_u$ is given as follows:

$$Q_u=\left[\begin{array}{cccc}-0.83& 0.32& 0.32& 0.19\\ 1460& -1460.61& 0& 0.51\\ 1460& 0& -730.51& 0.51\\ 0& 730& 730& -2920\end{array}\right]$$

(13)

with

$${\overline{\Pi }}_2=\left[\begin{array}{c}9.994\times {10}^{-1}\\ 1.894\times {10}^{-4}\\ 3.787\times {10}^{-4}\\ 5.91\times {10}^{-5}\end{array}\right]$$

(14)

By calculating the sum over the probabilities of the working states (in this example, the sum of $\left[{\pi }_1\right]$, $\left[{\pi }_2\right]$, and $\left[{\pi }_3\right]$), the availability is $\left[A_2\right]=[0.999302,1.00005]$.

Table 1. The imprecise availability obtained by the two methods.

The article’s method	$\left[A_1\right]=[0.99902,1.00085]$	$\cap [0,1]=[0.99902,1]$
The Exact method	$\left[A_2\right]=[0.99930,1.00005]$	$\cap [0,1]=[0.99930,1]$

shows the availability obtained by the method proposed in the article and the availability obtained by the exact method. To show the efficiency of the exact method, we keep the obtained interval as is. However, we should perform the intersection of the obtained interval with $[0,1]$ since the probability must belong to this interval, or we can also normalize the probabilities of being in each state. Note that in this example, and in the rest of the manuscript, the results are obtained after rounding off calculations.

$Q_l\in [\underline{Q},\overline{Q}]$ gives the lower stationary vector $\underline{\Pi }$, which gives the lower bound of the system’s availability interval. Thus, $Q_u\in [\underline{Q},\overline{Q}]$ gives the upper bound $\overline{\Pi }$, which gives the upper bound of the system’s availability interval. Therefore, it is not as claimed in the article by Troffaes that the lower and upper bounds of $[\Pi ]$ are obtained using $\underline{Q}$, all of whose elements are the lower bounds of the intervals, and $\overline{Q}$, all of whose elements are the upper bounds of the intervals. Note that in the article, the authors do not consider the condition that the sum of the probabilities is equal to 1. Later, in Section 2.3, we test our proposed method against the exact method to solve the problem by considering the lower and the upper bounds of $q_{ij}$.

2.3. Exact Method

In this subsection, we propose a method that we call the exact method since it is the most reliable method for providing a narrow interval of the system availability, and because it contains the real exact value of the availability, which is otherwise impossible to find. We need to take all the possible lower bound or upper bound values of the elements $\left[q_{ij}\right]$ of the transition matrix (combinations of the bounds of $\left[\lambda \right]$, the failure rates intervals, and $\left[\mu \right]$, the repair rate intervals of the system components based on their position in the matrix), which means all the possible combinations k. For each combination k, we find the corresponding transition matrix $Q_k$, and we turn the problem into the precise case in order to solve Equation (1) and to find the corresponding stationary vector ${\Pi }_k$. Using the vector ${\Pi }_k$, we can calculate the corresponding availability $A_k$ of the system for the combination k. The number of possible combinations is $2^l$, where l is the number of $\left[q_{ij}\right]$ elements of the transition matrix $\left[Q\right]$, and $k\in \left\{1,\dots ,2^l\right\}$. The bounds of the system availability are obtained by taking the maximal and minimal availabilities $A_k$ of all the combinations.

Example: Consider a series system formed by two components, A and B. Each of these components has two possible states: perfect functioning and total failure. Each component has the following failure and repair rates per hour:

$\left\{\begin{array}{c}{\lambda }_A=[2\times {10}^{-3},4\times {10}^{-3}]\hfill \\ {\mu }_A=[1.8\times {10}^{-2},3.5\times {10}^{-2}]\hfill \\ {\lambda }_B=[3\times {10}^{-3},4.5\times {10}^{-3}]\hfill \\ {\mu }_B=[1.5\times {10}^{-2},3\times {10}^{-2}]\hfill \end{array}\right.$

The system only works if both of the components work. The system has four possible states:

• State 1: Components A and B work.
• State 2: Component A works and B fails.
• State 3: Component B works and A fails.
• State 4: Components A and B fail.

Since the system is a series system, which means that the system works as long as all the components work, we can regroup the system’s states into two groups: working state (state 1) and failure state (states 2, 3, and 4).

The transition matrix elements $q_{ij}$ can take either the lower or the upper bound of the interval $\left[q_{ij}\right]$ (the combination of the failure and the repair rates of A and B). For example, the first combination $k=1$ is when all the $q_{ij}$ take the lower bounds of their intervals as values. The second combination $k=2$ is when $q_{11}$ takes the upper bound of its intervals as values and the rest take the lower bounds as values, and so on. In this case, we have $l=16$ and $2^l=2^{16}$ possible precise transition matrices $Q_k$, where $k\in \left\{1,\dots ,2^{16}\right\}$ is in the following form:

$$Q_k=\left[\begin{array}{cccc}-({\lambda }_A+{\lambda }_B)& {\lambda }_A& {\lambda }_B& 0\\ {\mu }_A& -({\mu }_A+{\lambda }_B)& 0& {\lambda }_B\\ {\mu }_B& 0& -({\lambda }_A+{\mu }_B)& {\lambda }_A\\ 0& {\mu }_B& {\mu }_A& -({\mu }_A+{\mu }_B)\end{array}\right]$$

(15)

For each $Q_k$, we solve ${\Pi }_k·Q_k=0$ to find the vector ${\Pi }_k$ and the availability $A_k$. When comparing all the k availabilities, we reform the interval of the availability of the system, where the lower bound has the smallest value from all the $A_k$s and the upper bound has the highest value from all of the $A_k$s.

In this example, the availability of the system is

$A=[0.5207,0.9427]$.

where the corresponding $Q_l$ is the transition matrix that gives $\underline{A}$, and $Q_u$ is the transition matrix that gives $\overline{A}$ ($Q_l$ and $Q_u$ are not the lower and upper transition matrices of Q).

$Q_l=\left[\begin{array}{cccc}-0.005& 0.002& 0.003& 1\\ 0.018& -0.038& 0& 1\\ 0.015& 0& -0.032& 1\\ 0& 0.015& 0.018& 1\end{array}\right]$.

$Q_u=\left[\begin{array}{cccc}-0.005& 0.004& 0.0045& 1\\ 0.035& -0.038& 0& 1\\ 0.03& 0& -0.032& 1\\ 0& 0.015& 0.018& 1\end{array}\right]$.

In conclusion, the exact method is time consuming and computation of the availability at each time is not straightforward.

3. Interval Constraint Propagation

The propagation of interval constraints is done by combining interval arithmetic computation [41] and constraint propagation [42]. The constraints are then used in order to eliminate all the inconsistent values defined in the variables domains.

3.1. Definitions, Notations, and Basic Operations

A real interval $\left[x\right]$ is defined as a connected and closed subset of $\mathbb{R}$. It is expressed as follows:

$$\left[x\right]=[\underline{x},\overline{x}]=\left\{x\in \mathbb{R}|\underline{x}\le x\le \overline{x}\right\},$$

(16)

where $\underline{x}$ and $\overline{x}$ represents the minimal and maximal bounds of $\left[x\right]$, respectively.

A box is a vector [x] in ${\mathbb{R}}^n$. It can be expressed as the Cartesian product of n intervals:

$$\left[x\right]=\left[x_1\right]\times ...\times \left[x_n\right].$$

(17)

The classical set operation (intersection and union) can be used in interval arithmetics. The interval union of the intervals $\left[x\right],\left[y\right]\subset \mathbb{R}$ is expressed as

$$\left[x\right]\bigcup \left[y\right]=\left[\right[x]\cup [y\left]\right]$$

(18)

where ∪ refers to the set union operation. The symbol $[.]$ denotes that for any set S in $\mathbb{R}$, the interval null operator will return the smallest interval enclosing S. The classical binary operations $\left\{+,-,\times ,÷\right\}$ are also extended to intervals. The resulting interval $\left[z\right]=\left[x\right]\diamond \left[y\right]$ is expressed by

$$\left[z\right]=\left[x\right]\diamond \left[y\right]=\left[x\diamond y|x\in [x],y\in [y]\right].$$

(19)

The detailed extensions of arithmetic operations are given by

$$\begin{array}{c}\hfill \left[x\right]+\left[y\right]=[\underline{x}+\underline{y},\overline{x}+\overline{y}],\end{array}$$

(20)

$$\begin{array}{c}\hfill \left[x\right]-\left[y\right]=[\underline{x}-\overline{y},\overline{x}-\underline{y}],\end{array}$$

(21)

$$\begin{array}{c}\hfill \left[x\right]\times \left[y\right]=[min(\underline{x}\underline{y},\underline{x}\overline{y},\overline{x}\underline{y},\overline{x}\overline{y}),max(\underline{x}\underline{y},\underline{x}\overline{y},\overline{x}\underline{y},\overline{x}\overline{y})],\end{array}$$

(22)

$$\begin{array}{c}\hfill \left[x\right]÷\left[y\right]=\left[x\right]\times [1/\overline{y},1/\underline{y}],0\notin \left[y\right]\end{array}$$

(23)

3.2. Contractors

The technique of contractors was first introduced in some applications of artificial intelligence in the 1970s and was considered as a constraints satisfaction [42]. The constraints satisfaction is the process of finding a solution to a set of constraints that impose conditions that the variables must satisfy. Then, the technique of contractors becomes an important parts of the interval analysis field [42]. This technique helps to contract a big interval $\left[x\right]$ into a smaller one $\left[x^{\prime }\right]$.

Let us consider $n_x$ variables: $x_i\in \mathbb{R}$ where $i\in 1,..,n_x$. These variables are linked by $n_f$ constraints, which have the following form:

$$f_j(x_1,..,x_{n_x})=0,j\in 1,..,n_f.$$

(24)

The variables $x_i$ belong to a domain ${\mathbb{X}}_i$ (these domains are represented by intervals denoted by $\left[x_i\right]$). A vector x is defined as follows:

$$x={(x_1,..,x_{n_x})}^T$$

(25)

The prior domain of x is represented by the box:

$$\left[x\right]=\left[x_1\right]\times ..\times \left[x_{n_x}\right]$$

(26)

Let us consider a function f whose coordinates are $f_j$s. The first equation has the form $f\left(x\right)=0$. This form represents a constraint satisfaction problem (CSP) H [42] that is expressed as follows:

$$H:\left(f\right(x)=0,x\in [x\left]\right)$$

(27)

The solution set of H is defined as

$$\mathbb{S}=x\in \left[x\right]\left|f\right(x)=0$$

(28)

The contraction of H consists in replacing [x] by a smaller domain [x’] such that the solution set remains unchanged, i.e., $\mathbb{S}\subset \left[x^{\prime }\right]\subset \left[x\right]$. It was proven that there exists an optimal contraction of H that replaces [x] by the smallest box containing $\mathbb{S}$. A contractor of H is any operator used for the contraction of H.

3.3. Forward–Backward Propagation

Many contractors were proposed in the literature. Each contractor can be efficient for one specific CSPs and has its own algorithm of functioning. The most popular contractor is the forward–backward propagation (FBP) contractor. We chose this contractor because it is known for its ease and simplicity, and many researchers consider that it is more general than other contractors since it works for all types of systems. It also leads to guaranteed results (i.e., when performing the contractions, we always obtain an interval that belongs to the initial interval). Compared to other techniques, the FBP offers great accuracy. For all these reasons, we chose to use the FBP technique to help us to contract the intervals with the aim of finding the availability when t tends to inf, as given in Equation (1) for the case of imprecise data.

The forward–backward contractor $C_{\downarrow \uparrow }$ contracts the CSP domains H by considering each constraint of the constraints $n_f$ separately. It is based on two steps: forward and backward steps. In the forward step, interval arithmetic operations are applied to all the operators of the function $y=f\left(x\right)$ (from the variable’s domain ([x]) up to the function’s domain ([y])). In other words, it is applied to the direct forms of equations. In the backward step, it sets the interval associated with the new function’s domain [y] to $[0,0]$. Then, it applies backward arithmetic from the function’s domain to the variable’s domain. In other words, it uses the inverse of the functions that appear in the equations f(x).

The following example explains the procedure of the FBP technique:

Let us consider the constraint $y=-5x_1+2x_2=0$ and the initial box-domain $\left[x\right]=\left[1,4\right]\times \left[-3,7\right]$.

This constraint can be decomposed, as shown in Equation (29), into three primitive constraints (i.e., constraints involving a single operator, such as ($+,-,\ast ,/$), or a single function) by introducing two intermediate variables, $z_1$ and $z_2$, defined as $z_1=-5·x_1$ and $z_2=2·x_2$, respectively. The initial domains for these variables are determined as follows:

$$\left\{\begin{array}{c}\left[z_1\right]:=-5\left[x_1\right]=-5\times [1,4]=[-20,-5]\hfill \\ \left[z_2\right]:=2\left[x_2\right]=2\times [-3,7]=[-6,14]\hfill \\ \left[y\right]:=\left[z_1\right]+\left[z_2\right]=[-20,-5]+[-6,14]=[-26,9]\hfill \end{array}\right.$$

(29)

We aim to contract H while respecting the constraint $f\left(x\right)=5x_1+2x_2=0$. Thus, we contract the primitive constraints presented in Equation (29) until the contractions become inefficient (i.e., when the solution converges).

In this case, the domain for y is equal to 0 because $f\left(x\right)=0$. Then, we add the following step:

$$\left[y\right]:=\left[y\right]\cap 0$$

(30)

If $\left[y\right]$ is empty, we can say that the CSP does not have any solutions. Otherwise, $\left[y\right]$ is replaced by the value 0. Then, we apply the backward propagation from the function’s domain ($z_1$ and $z_2$) to the variable’s domain ($x_1$ and $x_2$). Formally, all the domains associated to the variables should be updated as follows:

$$\left\{\begin{array}{c}\left[z_1\right]:=(\left[y\right]-\left[z_2\right])\cap \left[z_1\right]=([0,0]-[-6,14])\cap [-20,-5]⟹\left[z_1\right]=[-14,-5]\hfill \\ \left[z_2\right]:=(\left[y\right]-\left[z_1\right])\cap \left[a_2\right]=([0,0]-[-20,-5])\cap [-6,14]⟹\left[z_2\right]=[5,14]\hfill \\ \left[x_2\right]:=(\left[z_1\right]/-5)\cap \left[x_1\right]=([-14,-5]/-5)\cap [1,4]⟹\left[x_1\right]=[1,14/5]\hfill \\ \left[x_2\right]:=(\left[z_2\right]/2)\cap \left[x_2\right]=([5,14]/2)\cap [-3,7]⟹\left[x_2\right]=[5/2,7].\hfill \end{array}\right\}$$

(31)

The obtained box after the first contraction is:

$$\left[x\right]\left(1\right)=[1,14/15]\times [5/2,7]$$

(32)

We then repeat the obvious procedure until the boxes $\left[x\right]\left(k\right)$ converge to the smallest possible domain.

Table 2. Steps of $C_{\downarrow \uparrow }$.

Step 1: Forward	$\left[z_1\right]:=-5\left[x_1\right]$
Step 2: Forward	$\left[z_2\right]:=2\left[x_2\right]$
Step 3: Forward	$\left[y\right]:=\left[z_1\right]+\left[z_2\right]$
Step 4	$\left[y\right]:=\left[y\right]\cap 0$
Step 5: Backward	$\left[z_1\right]:=(\left[y\right]-\left[z_2\right])\cap \left[z_1\right]$
Step 6: Backward	$\left[z_2\right]:=(\left[y\right]-\left[z_1\right])\cap \left[z_2\right]$
Step 7: Backward	$\left[x_1\right]:=(\left[z_1\right]/-5)\cap \left[x_1\right]$
Step 8: Backward	$\left[x_2\right]:=(\left[z_2\right]/2)\cap \left[x_2\right]$

shows the corresponding steps of the FBP contractor applied to this example.

4. The Proposed Methodology

Our main purpose is to calculate the availability of an MSS when t tends to ∞. The uncertainties of the data (failure and repair rates) must be taken into account. The first method, presented in Section 2, cannot be used since the authors found the lower bound of the stationary vector by taking the lower bound of the transition matrix and the upper bound of the stationary vector by taking the upper bound of the same matrix. By detailing the calculation, we show that another matrix may exist that gives the upper and the lower bounds of $[\Pi ]$.

The exact method is the most efficient method for finding the imprecise availability of an MSS. However, when the system becomes much larger, the number of transitions becomes enormous. In this case, we have a combinatorial explosion of the calculation. To avoid this, we need a guaranteed method that helps us to efficiently calculate the imprecise availability of the MSS in order to guarantee that our solution includes the true value of the availability. Our search is aimed at simple calculation procedures with a minimal time cost.

We choose to model the uncertainties of the data in terms of intervals, and we use imprecise Markov models and introduce the contracting technique (cf. Section 3.3) to compute the imprecise availability of an MSS.

In this section, we give the detailed steps of our methodology. To better demonstrate this methodology, we present the different steps through simple examples.

4.1. The Steps of the Methodology

In order to compute the MSS availability, we should follow and apply the steps presented in Figure 2, summarized below:

1.

Step 1: Define the Markov graph to present the states and all transitions between the different system states.

2.

Step 2: Construct the interval transition matrix of the system.

3.

Step 3: Apply the contracting technique to the system of equations $\Pi ·Q=0$.

4.

Step 4: Find $\Pi$ and the availability of the system.

Before we start with the steps of the method, we need to list all the information about the system we are considering. What type of system do we have and what is its structure? How many components does it have and how many states does each component have? To answer these questions, we determine the system/components states, and we define the corresponding interval failure rates $\lambda$ and the interval repair rates $\mu$ of each component. We then determine the type of connection between one component’s state and the other states. In this way, we can easily define the possible states of the system and understand the types of connections between them.

We should mention that when the components are independent, we can take each component separately and apply the steps of the methodology to compute its availability before computing the availability of the entire system. The interval obtained in this way is more conservative than the availability interval obtained when applying the methodology to the entire system. The following example illustrates the different steps.

Imagine a parallel system of two components, A and B. A parallel system works when at least one of the two components works. For each component, we have two possible states (binary components): state 1 (perfect functioning) and state 2 (total failure). Thus, the whole system has four possible states, where the first three states represent the functioning state of the parallel system and the last state is the failure state. The components have the following failure and repair rates per hour:

$\left\{\begin{array}{c}{\lambda }_A=[3\times {10}^{-3},4\times {10}^{-3}];{\mu }_A=[2.8\times {10}^{-2},3.5\times {10}^{-2}]\hfill \\ {\lambda }_B=[3.5\times {10}^{-3},4.5\times {10}^{-3}];{\mu }_B=[2\times {10}^{-2},3\times {10}^{-2}]\hfill \end{array}\right.$

4.1.1. Step 1

Given the system structure and the states of each component, we can easily find the states of the whole system. In the above example, we have two binary components. Thus, the system has four states: total functioning (state 1) and degraded states (states 2 and 3), which are all considered functioning states in the parallel system, and the last state, total failure (state 4). Figure 3 shows the Markov graph of this system.

4.1.2. Step 2

Once we have the Markov graph, either for the system or for each component, we move on to the second step. As in the precise case, we can determine the transition matrix from the Markov graph, the only difference being that we are dealing with uncertainties in terms of intervals. We define the interval transition matrix $\left[Q\right]$, where each term of the matrix $q_{ij}=[\underline{q_{ij}},\overline{q_{ij}}]$ is an interval of failure rates ${\lambda }_i=[\underline{{\lambda }_i},\overline{{\lambda }_i}]$ and an interval of repair rates ${\mu }_i=[\underline{{\mu }_i},\overline{{\mu }_i}]$ of dimension $(n\times n)$.

We should mention that each term of the transition matrix represents the transition rate from one state to another. The transition matrix of our example is

$$\left[Q\right]=[\underline{Q},\overline{Q}]=\left[\begin{array}{cccc}[\underline{q_{11}},\overline{q_{11}}]& [\underline{q_{12}},\overline{q_{12}}]& [\underline{q_{13}},\overline{q_{13}}]& [\underline{q_{14}},\overline{q_{14}}]\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ \left[\underline{q_{41}},\overline{q_{41}}\right]& [\underline{q_{42}},\overline{q_{42}}]& [\underline{q_{43}},\overline{q_{43}}]& [\underline{q_{44}},\overline{q_{44}}]\end{array}\right]$$

(33)

In our example, $\left[Q\right]$ is expressed as

$\left[\begin{array}{cccc}[-(\overline{{\lambda }_A}+\overline{{\lambda }_B}),-(\underline{{\lambda }_A}+\underline{{\lambda }_B})]& [\underline{{\lambda }_A},\overline{{\lambda }_A}]& [\underline{{\lambda }_B},\overline{{\lambda }_B}]& 0\\ [\underline{{\mu }_A},\overline{{\mu }_A}]& [-(\overline{{\mu }_A}+\overline{{\lambda }_B}),-(\underline{{\mu }_A}+\underline{{\lambda }_B})]& 0& [\underline{{\lambda }_B},\overline{{\lambda }_B}]\\ [\underline{{\mu }_B},\overline{{\mu }_B}]& 0& [-(\overline{{\lambda }_A}+\overline{{\mu }_B}),-(\underline{{\lambda }_A}+\underline{{\mu }_B})]& [\underline{{\lambda }_A},\overline{{\lambda }_A}]\\ 0& [\underline{{\mu }_B},\overline{{\mu }_B}]& [\underline{{\mu }_A},\overline{{\mu }_A}]& [-(\overline{{\mu }_A}+\overline{{\mu }_B}),-(\underline{{\mu }_A}+\underline{{\mu }_B})]\end{array}\right]$

By replacing the values of failure and repair rates, we obtain

$\left[\begin{array}{cccc}[-8.5\times {10}^{-3},-6.5\times {10}^{-3}]& [3\times {10}^{-3},4\times {10}^{-3}]& [3.5\times {10}^{-3},4.5\times {10}^{-3}]& 0\\ [2.8\times {10}^{-2},3.5\times {10}^{-2}]& [-3.95\times {10}^{-2},-3.15\times {10}^{-2}]& 0& [3.5\times {10}^{-3},4.5\times {10}^{-3}]\\ [2\times {10}^{-2},3\times {10}^{-2}]& 0& [-3.4\times {10}^{-2},-2.3\times {10}^{-2}]& [3\times {10}^{-3},4\times {10}^{-3}]\\ 0& [2\times {10}^{-2},3\times {10}^{-2}]& [2.8\times {10}^{-2},3.5\times {10}^{-2}]& [-6.5\times {10}^{-2},-4.8\times {10}^{-2}]\end{array}\right]$

4.1.3. Step 3

After we have our transition matrix $\left[Q\right]$, we move on to the next step. To determine the availability of the system, we must have the probabilities of being in each state of the system. Since we want to find the availability when the time t tends to ∞, we need to find the steady probabilities of being in each state of the system:

$[\Pi ]=[\underline{\Pi },\overline{\Pi }]=[[\underline{{\pi }_1},\overline{{\pi }_1}]\dots [\underline{{\pi }_n},\overline{{\pi }_n}]]$.

To do this, and since this vector is unknown to us, we assume an initial vector $[\Pi ]$ expressed as

$$[\Pi ]=\left[\right[0,1\left]\right[0,1]\dots [0,1\left]\right]$$

(34)

where $[\Pi ]$ presents all the initial intervals of steady probability of the system to be in each state. To compute the intervals of this vector, we use the contracting technique, which is why we apply Equation (1). We then obtain a system of n equations. The last equation is related to the following constraint:

$$\sum _{i=1}^n{\pi }_i=1,$$

(35)

which means that the sum of all probabilities is equal to 1. By performing the operations of the FBP contractor, we are able to contract our steady probabilities into smaller intervals that contain the unknown but real exact values.

In the previous example, our system of equations is as follows:

$$\left\{\begin{array}{c}\left[q_{11}\right]\left[{\pi }_1\right]+\left[q_{21}\right]\left[{\pi }_2\right]+\left[q_{31}\right]\left[{\pi }_3\right]=0\Rightarrow \left(1\right)\hfill \\ [q_{12}]\left[{\pi }_1\right]+\left[q_{22}\right]\left[{\pi }_2\right]+\left[q_{41}\right]\left[{\pi }_4\right]=0\Rightarrow \left(2\right)\hfill \\ [q_{13}]\left[{\pi }_1\right]+\left[q_{33}\right]\left[{\pi }_3\right]+\left[q_{43}\right]\left[{\pi }_4\right]=0\Rightarrow \left(3\right)\hfill \\ [q_{24}]\left[{\pi }_2\right]+\left[q_{34}\right]\left[{\pi }_3\right]+\left[q_{44}\right]\left[{\pi }_4\right]=0\Rightarrow \left(4\right)\hfill \\ [{\pi }_1]+\left[{\pi }_2\right]+\left[{\pi }_3\right]+\left[{\pi }_4\right]=1\Rightarrow \left(5\right)\hfill \end{array}\right.$$

(36)

We apply FBP on the intervals $\left[{\pi }_i\right]=\left[\underline{{\pi }_i},\overline{{\pi }_i}\right]$ to reduce each $\left[{\pi }_i\right]$ as much as possible, which means reducing the elements of the vector $\Pi$.

We start with the first contraction. Each constraint (equation) in Equation (36) has the following form:

$$f_i({\pi }_1,\dots ,{\pi }_{n_k})=0$$

(37)

where $f_i$ is broken down into a sequence of elementary operations (i.e., operations with only elementary operators), such as ($+,-,\ast ,/$). We break down this constraint into constraints [47] involving only single operators. For example, the first constraint in Equation (36) can be decomposed into the following primitive constraints:

$$\left\{\begin{array}{c}\left[a_1\right]=\left[q_{11}\right]\left[{\pi }_1\right]=[-8.5\times {10}^{-3},-6.5\times {10}^{-3}]\ast [0,1]=[-8.5\times {10}^{-3},0],\hfill \\ \left[a_2\right]=\left[q_{21}\right]\left[{\pi }_2\right]=[2.8\times {10}^{-2},3.5\times {10}^{-2}]\ast [0,1]=[0,3.5\times {10}^{-2}],\hfill \\ \left[a_3\right]=\left[q_{31}\right]\left[{\pi }_3\right]=[2\times {10}^{-2},3\times {10}^{-2}]\ast [0,1]=[0,3\times {10}^{-2}],\hfill \\ \left[y\right]=\left[a_1\right]+\left[a_2\right]+\left[a_3\right]=[-8.5\times {10}^{-3},6.5\times {10}^{-2}]\hfill \end{array}\right.$$

(38)

After we obtain the new interval [y], we find the intersection with 0 because the main constraint is $y=0$:

$\left[y\right]\cap 0=[-8.5\times {10}^{-3},6.5\times {10}^{-2}]\cap [0,0]=[0,0]$.

With this final value of [y], we perform the backward propagation by re-evaluating the three intervals $\left[a_1\right],\left[a_2\right]$, and $\left[a_3\right]$ such that, in the end, we can reconstruct the required intervals $\left[{\pi }_1\right],\left[{\pi }_2\right]$, and $\left[{\pi }_3\right]$, as shown in Equation (39):

$$\left\{\begin{array}{c}\left[a_1\right]=(\left[y\right]-\left[a_2\right]-\left[a_3\right])\cap \left[a_1\right]\hfill \\ \left[a_2\right]=(\left[y\right]-\left[a_1\right]-\left[a_3\right])\cap \left[a_2\right]\hfill \\ \left[a_3\right]=(\left[y\right]-\left[a_1\right]-\left[a_2\right])\cap \left[a_3\right]\hfill \\ \left[{\pi }_1\right]=(\left[a_1\right]/q_{11})\cap \left[{\pi }_1\right]\hfill \\ \left[{\pi }_2\right]=(\left[a_2\right]/q_{21})\cap \left[{\pi }_2\right]\hfill \\ \left[{\pi }_3\right]=(\left[a_3\right]/q_{31})\cap \left[{\pi }_3\right]\hfill \end{array}\right.$$

(39)

The intervals obtained are

$$\left\{\begin{array}{c}\left[a_1\right]=[-6.5\times {10}^{-2},0]\cap [-8.5\times {10}^{-3},0]=[-8.5\times {10}^{-3},0]\hfill \\ \left[a_2\right]=[-3\times {10}^{-2},8.5\times {10}^{-3}]\cap [0,3.5\times {10}^{-2}]=[0,8.5\times {10}^{-3}]\hfill \\ \left[a_3\right]=[-3.5\times {10}^{-2},8.5\times {10}^{-3}]\cap [0,3\times {10}^{-2}]=[0,8.5\times {10}^{-3}]\hfill \\ [{\pi }_1]=[0,1.30769]\cap [0,1]=[0,1]\hfill \\ [{\pi }_2]=[0,0.3035714]\cap [0,1]=[0,0.3035714]\hfill \\ [{\pi }_3]=[0,0.425]\cap [0,1]=[0,0.425]\hfill \end{array}\right.$$

(40)

After applying constraint number 1, we obtain the vector $\Pi$ in the following form:

$$\left[{\Pi }_1\right]=\left[\begin{array}{cccc}[0,1]& [0,0.3035714]& [0,0.425]& [0,1]\end{array}\right]$$

(41)

We continue the technique of contraction with constraint number 2. By following the same steps, we have the primitive constraints

$$\left\{\begin{array}{c}\left[b_1\right]=\left[q_{12}\right]\left[{\pi }_1\right]=[3\times {10}^{-3},4\times {10}^{-3}]\ast [0,1]=[0,4\times {10}^{-3}],\hfill \\ \left[b_2\right]=\left[q_{22}\right]\left[{\pi }_2\right]=[-3.95\times {10}^{-2},-3.15\times {10}^{-2}]\ast [0,0.3035714]=[-1.199107\times {10}^{-2},0],\hfill \\ \left[b_4\right]=\left[q_{42}\right]\left[{\pi }_4\right]=[2\times {10}^{-2},3\times {10}^{-2}]\ast [0,1]=[0,3\times {10}^{-2}],\hfill \\ \left[y\right]=\left[b_1\right]+\left[b_2\right]+\left[b_4\right]=[-1.199107\times {10}^{-2},7\times {10}^{-2}]\hfill \end{array}\right.$$

(42)

The intersection of [y] with $[0,0]$ yields $[0,0]$. We now repeat the same operations backwards to re-evaluate $\left[b_1\right],\left[b_2\right]$, and $\left[b_4\right]$, and the newest version of $\left[{\pi }_1\right],\left[{\pi }_2\right]$, and $\left[{\pi }_4\right]$. We find

$$\left\{\begin{array}{c}\left[b_1\right]=(\left[y\right]-\left[b_2\right]-\left[b_4\right])\cap \left[b_1\right]\hfill \\ \left[b_2\right]=(\left[y\right]-\left[b_1\right]-\left[b_4\right])\cap \left[b_2\right]\hfill \\ \left[b_4\right]=(\left[y\right]-\left[b_1\right]-\left[b_2\right])\cap \left[b_4\right]\hfill \\ \left[{\pi }_1\right]=(\left[b_1\right]/q_{12})\cap \left[{\pi }_1\right]\hfill \\ \left[{\pi }_2\right]=(\left[b_2\right]/q_{22})\cap \left[{\pi }_2\right]\hfill \\ \left[{\pi }_4\right]=(\left[b_4\right]/q_{42})\cap \left[{\pi }_4\right]\hfill \end{array}\right.$$

(43)

where

$$\left\{\begin{array}{c}\left[b_1\right]=[-3\times {10}^{-2},1.199107\times {10}^{-2}]\cap [0,4\times {10}^{-3}]=[0,4\times {10}^{-3}]\hfill \\ [b_2=[-3.4\times {10}^{-2},0]\cap [-1.199107\times {10}^{-2},0]=[-1.199107\times {10}^{-2},0]\hfill \\ \left[b_4\right]=[-4\times {10}^{-3},1.199107\times {10}^{-2}]\cap [0,3\times {10}^{-2}]=[0,1.199107\times {10}^{-2}]\hfill \\ [{\pi }_1]=[0,1.33333]\cap [0,1]=[0,1]\hfill \\ [{\pi }_2]=[0,0.3807]\cap [0,0.3035714]=[0,0.3035714]\hfill \\ [{\pi }_4]=[0,0.05996]\cap [0,1]=[0,0.05996]\hfill \end{array}\right.$$

(44)

The newest vector of probabilities is, therefore,

$$\left[{\Pi }_2\right]=\left[\begin{array}{cccc}[0,1]& [0,0.3035714]& [0,0.425]& [0,0.05996]\end{array}\right]$$

(45)

We then take the third, fourth, and fifth (or final) constraints, applying the same steps as before.

Table 3. The probability vector obtained after each constraint.

Constraint Number	Probability Vector $\Pi$
1	$\left[\begin{array}{cccc}[0,1]& [0,0.303]& [0,0.425]& [0,1]\end{array}\right]$
2	$\left[\begin{array}{cccc}[0,1]& [0,0.303]& [0,0.425]& [0,0.059]\end{array}\right]$
3	$\left[\begin{array}{cccc}[0,1]& [0,0.303]& [0,0.032]& [0,0.599]\end{array}\right]$
4	$\left[\begin{array}{cccc}[0,1]& [0,0.303]& [0,0.425]& [0,0.006]\end{array}\right]$
5	$\left[\begin{array}{cccc}[0.207,1]& [0,0.303]& [0,0.425]& [0,0.063]\end{array}\right]$

represents the probability vector $\Pi$ that is obtained after using the five constraints and applying the contracting technique only one time. We continue repeating the contractions on all the equations until $\pi$, the elements of the vector $\Pi$, converge. This allows us to obtain the final version of the interval vector $\Pi$.

Table 4. The probability vector obtained after k contractions.

Number of Contractions k	Probability Vector $\Pi$
1	$\left[\begin{array}{cccc}[0.207,1]& [0,0.303]& [0,0.425]& [0,0.063]\end{array}\right]$
2	$\left[\begin{array}{cccc}[0.496,0.959]& [0.016,0.182]& [0.022,0.280]& [0.001,0.040]\end{array}\right]$
3	$\left[\begin{array}{cccc}[0.571,0.901]& [0.039,0.155]& [0.054,0.238]& [0.004,0.034]\end{array}\right]$
6	$\left[\begin{array}{cccc}[0.625,0.871]& [0.051,0.135]& [0.070,0.208]& [0.006,0.030]\end{array}\right]$
10	$\left[\begin{array}{cccc}[0.628,0.869]& [0.052,0.134]& [0.071,0.206]& [0.006,0.029]\end{array}\right]$

presents the vector $\Pi$ after k contractions.

Figure 4 shows how the probability interval for each state changes with each contraction.

4.1.4. Step 4

The system availability $A_s$ is given by:

$$A_s=\sum _{i,\mathrm{over}\mathrm{the}\mathrm{working}\mathrm{states}}\left[{\pi }_i\right]$$

(46)

In the last example, we had a parallel system of two components. ${\pi }_1,{\pi }_2$, and ${\pi }_3$ are the probabilities that correspond to the working states. ${\pi }_4$ is the probability that corresponds to the failure state of the system.

On the basis of the last version of $[\Pi ]$ obtained after ten contractions (when each $\left[{\pi }_i\right]$ converges to a certain interval), we compute the availability of the system:

$$\left\{\begin{array}{c}\left[{\pi }_1\right]+\left[{\pi }_2\right]+\left[{\pi }_3\right]=[0.7529529,1.21107924]\hfill \\ A_s=[0.7529529,1]\hfill \end{array}\right.$$

(47)

Note that we can find the availability of the system when its components are independent by taking each component separately and applying the same steps mentioned in this section. By performing the steps for each component, we find the corresponding $\left[{\Pi }_j\right]$ of the component j, which allows us to compute the availability of each component $\left[A_j\right]$ separately. After determining the availability of all the components, we determine the relationship between all the $\left[A_j\right]$, and we are able to calculate the imprecise availability of the system $\left[A_s\right]$. In this example, we have a parallel system, which means that the relationship between the availability of the components is

$$A_s=1-(1-A_1)(1-A_2)$$

(48)

To verify if the result obtained by the contracting technique is accurate, we apply the exact method and compare the calculated availability with the answer obtained from the contracting technique. Recall that the main idea of the exact method is to take all the possible transition matrices by combining the lower and upper bounds of $\left[q_{ij}\right]$ (i.e., taking the lower and upper bounds of the intervals of the failure and repair rates).

When we look at the previous example using the exact method, we obtain the lower bound of the steady probabilities:

${\underline{\Pi }}_E=\left[\begin{array}{c}0.666963\hfill \\ 0.098112\hfill \\ 0.146102\hfill \\ 0.088821\hfill \end{array}\right]$

and the upper bound of the steady probabilities:

${\overline{\Pi }}_E=\left[\begin{array}{c}0.837987\hfill \\ 0.060367\hfill \\ 0.111134\hfill \\ 0.009490\hfill \end{array}\right]$

The availability determined by the exact method is

$\left[A_E\right]=[0.911,1]$

We also consider the “precise method” by considering the center of each transition matrix $\left[q_{ij}\right]$. In other words, we do not take into account the imprecision of the transition rates by simply considering the center value of the transition rates intervals $\frac{q_{ij}^{min}+q_{ij}^{max}}{2}$. When we consider the previous example using the precise method, we obtain $A_P=0.986$.

Table 5. The system’s availability obtained by different methods.

Contracting Technique	Exact Method	Precise Method
$A_s=[0.752,1]$	$\left[A_E\right]=[0.911,1]$	$A_P=0.986$

presents the results obtained by the three methods. As we can see, our proposed method is more conservative than the “exact method”.

5. Metrics for Validation of the Methodology

5.1. Guaranteed Index (GI)

In general, the results provided by the FBP contraction method are guaranteed, which means that on the basis of correct assumptions, the interval contraction provides correct conclusions.

More specifically, in the availability analysis, we consider that a guaranteed availability result is obtained when we are sure that the “real” availability belongs to the availability interval obtained. In order to quantify and validate the efficiency of our method, we consider that a guaranteed result is an interval that contains the exact interval availability. Recall that the exact interval availability is obtained when considering all the combinations of the transition matrix elements $q_{ij}$ (upper and lower failure and repair rates). We can compute the optimal interval from all the obtained intervals.

Thus, we introduce the guaranteed index (GI), which is a binary variable that is equal to 1 if the obtained availability interval contains the exact interval availability, and equal to 0 if it is not the case.

5.2. Conservatism Rate (CR)

We define the conservatism rate (CR) as the ratio of the difference between the width of the interval of the availability obtained by our proposed method and the width of the availability obtained by the exact method.

Formally, the CR is given by

$$CR=\frac{w\left(\left[A_{obtained}\right]\right)-w\left(A_E\right)}{w\left(A_{obtained}\right)}$$

(49)

where $w\left(\left[A\right]\right)=\overline{A}-\underline{A}$.

6. Case Studies

To understand and validate the proposed methodology, we propose two MSS case studies and apply the methodology to calculate the imprecise availability. After calculating the imprecise availability, we compare the results to the ones obtained using two well-known methodologies (IUGF [30] and BUGF [43]). The first case study shows the accuracy of our method, compared to other methods. The second one verifies that the method can efficiently handle cases of complex systems that have a large number of components and states.

6.1. Case Study 1

In this section, we use the example presented by Soroudi et al. [48]. Our aim is to compute the system availability described in [48] and presented in Figure 5. It is a flow transmission system from left to right and composed of three pipes. The pipes’ performance is computed using their transmission capacities (tons per minute).

The components 1 and 2 have 3 states: full capacity, partial failure (a partial capacity), and total failure (a capacity of 0). The component 3 has two states: full capacity and total failure. All the performance levels of the component’s states are precise. The system availability is computed using Markov models. We also compare this availability to those obtained by the IUGF method proposed in [30] and by the BUGF method proposed in [43].

In this case study, the repair and failure rates are in form of intervals (cf.

Table 6. Transition rates and states for each component.

	${\mathit{G}}_1$	${\mathit{G}}_2$	${\mathit{G}}_3$
State 1	$g_3^1=1.5$	$g_3^2=2$	$g_2^3=4$
	(Completely successful)	(Completely successful)	(Completely successful)
State 2	$g_2^1=1$	$g_2^2=1.5$	$g_1^3=0$
	(Degraded successful)	(Degraded successful)	(Total failure)
State 3	$g_1^1=0$	$g_1^2=0$	−
	(Total failure)	(Total failure)
Failure rates	${\lambda }_{3,1}^1=[{10}^{-5},3\times {10}^{-4}]$$h^{-1}$	${\lambda }_{3,1}^2=[2\times {10}^{-5},6\times {10}^{-4}]$$h^{-1}$	${\lambda }_{2,1}^3=[{10}^{-5},4\times {10}^{-4}]$$h^{-1}$
	${\lambda }_{2,1}^1=[4\times {10}^{-5},5\times {10}^{-4}]$$h^{-1}$	${\lambda }_{2,1}^2=[3\times {10}^{-5},4\times {10}^{-4}]$$h^{-1}$
Repair rates	${\mu }_{2,1}^1=[2\times {10}^{-2},5\times {10}^{-2}]$$h^{-1}$	${\mu }_{2,1}^2=[3\times {10}^{-2},6\times {10}^{-2}]$$h^{-1}$	${\mu }_{1,1}^3=[5\times {10}^{-2},9\times {10}^{-2}]$$h^{-1}$
	${\mu }_{1,1}^1=[4\times {10}^{-2},8\times {10}^{-2}]$$h^{-1}$	${\mu }_{1,1}^2=[3\times {10}^{-2},7\times {10}^{-2}]$$h^{-1}$

). The system has at most 18 $(3\times 3\times 2)$ possible states, ranging from the total functional state to the failure state (cf. Figure 6).

The notation $0_j$ indicates that component j is in the total functional state, the notation $1_j$ indicates that component j is in the partial working state, and the notation $2_j$ indicates that component j is in the state of total failure. The performance level g of a state is computed by considering the performance level of components. Thus, the state 1, which is denoted by $0_1{0}_2{0}_3$, refers to the state where all 3 components are totally functional. The components 1 and 2 are in parallel so that the obtained performance level is equal to the sum of the two performance instances: $2+1.5=3.5$. The component 3 is in series with the parallel configuration of 1 and 2. Thus, the system performance level g is equal to the minimum, i.e., $g=min(3.5,4)=3.5$. We compute the system availability when considering a demand level $w=1.5$. Consequently, we consider the system working states when the total performance level is greater than or equal to w. In Figure 6, the colored states represent working states.

The system availabilities calculated using the IUGF method and BUGF method are shown in

Table 7. Availability when t tends to ∞ for each method.

Exact Method	FBP Contraction Technique	BUGF	IUGF
$[0.9809,1.0323]$	$[0.9551,1.0360]$	[0.9550, 1.0361]	[0.9137, 1.0900]

. To obtain the availability using our proposed methodology, we first use the FBP technique to solve the system of equations by contracting the intervals. The system of Equation (1) is equal to zero (18 equations) plus one additional equation representing the fact that the sum of each steady availability term $\pi$ will be equal to 1. Following the steps of the methodology, we find our imprecise availability. The obtained results are represented in Table 7. The system availability of the exact method is also given in Figure 7 and Table 7.

From Table 7, we can see that the accuracy of the FBP method is close to that obtained when we use the exact method. Moreover, the interval availability obtained by BUGF and IUGF methods are also less conservative than the FBP method since the interval of the contracting technique contains the interval obtained by the exact method, is smaller than the intervals obtained by the IUGF method, and is very close to the one obtained by the BUGF method.

In this case, the GI is equal to 1 since the interval obtained using the FBP contraction technique contains the interval obtained using the exact method. Moreover, the CR of this example is

$$C{R}_1=\frac{0.0809-0.0514}{0.0809}=0.364$$

(50)

We conclude that the FBP method offers, in this case, accurate results regarding the results of the BUGF and IUGF methods (the BUGF and IUGF availability intervals contain the FBP availability interval).

6.2. Case Study 2

We have studied the accuracy of the FBP method regarding the exact, the IUG, and the BUGF methods. We now show that the proposed methodology also remains efficient when the system has a large number of states. To illustrate this, we apply our methodology to the following case study. We aim to compute the availability of the system presented in Figure 8. The system is composed of eight components: A, B, C, D, E, F, G, and H. The transmission flow of the system is considered to be from left to right. The four components B, C, E, and H have states of working (total working and total failure). The three components A, F, and G have three states of working: full working, partial failure, and total failure. The last component D has four possible states: total failure, degraded functioning of type 2, degraded functioning of type 1, and full working. We deduce that the number of system states is $4\times 3\times 3\times 3\times 2\times 2\times 2\times 2=1728$. The values of the failure and repair rates of each component are presented in

Table 8. Failure and repair rates from one state to another for each component.

Component	Failure Rate	Repair Rate
A	${\lambda }_{1,1}=[3\times {10}^{-5},4.5\times {10}^{-5}]$$h^{-1}$	${\mu }_{1,1}=[4.1\times {10}^{-1},7\times {10}^{-1}]$$h^{-1}$
	${\lambda }_{1,2}=[2.6\times {10}^{-5},7\times {10}^{-5}]$$h^{-1}$	${\mu }_{1,2}=[3\times {10}^{-1},6.2\times {10}^{-1}]$$h^{-1}$
B	${\lambda }_2=[2.6\times {10}^{-5},3.1\times {10}^{-5}]$$h^{-1}$	${\mu }_2=[1.6\times {10}^{-1},4\times {10}^{-1}]$$h^{-1}$
C	${\lambda }_3=[3.5\times {10}^{-5},4.6\times {10}^{-5}]$$h^{-1}$	${\mu }_3=[2.1\times {10}^{-1},3\times {10}^{-1}]$$h^{-1}$
D	${\lambda }_{4,1}=[2.6\times {10}^{-5},3.5\times {10}^{-5}]$$h^{-1}$	${\mu }_{4,1}=[3.1\times {10}^{-1},3.5\times {10}^{-1}]$$h^{-1}$
	${\lambda }_{4,2}=[2.9\times {10}^{-5},3.7\times {10}^{-5}]$$h^{-1}$	${\mu }_{4,2}=[2.9\times {10}^{-1},3.3\times {10}^{-1}]$$h^{-1}$
	${\lambda }_{4,3}=[2.7\times {10}^{-5},3.4\times {10}^{-5}]$$h^{-1}$	${\mu }_{4,3}=[3\times {10}^{-1},3.5\times {10}^{-1}]$$h^{-1}$
E	${\lambda }_5=[3\times {10}^{-5},4.1\times {10}^{-5}]$$h^{-1}$	${\mu }_5=[2.5\times {10}^{-1},3\times {10}^{-1}]$$h^{-1}$
F	${\lambda }_{6,1}=[3.1\times {10}^{-5},5.6\times {10}^{-5}]$$h^{-1}$	${\mu }_{6,1}=[2\times {10}^{-1},2.4\times {10}^{-1}]$$h^{-1}$
	${\lambda }_{6,2}=[2.9\times {10}^{-5},4.6\times {10}^{-5}]$$h^{-1}$	${\mu }_{6,2}=[2.6\times {10}^{-1},3\times {10}^{-1}]$$h^{-1}$
G	${\lambda }_{7,1}=[4.1\times {10}^{-5},4.7\times {10}^{-5}]$$h^{-1}$	${\mu }_{7,1}=[2.3\times {10}^{-1},2.9\times {10}^{-1}]$$h^{-1}$
	${\lambda }_{7,2}=[2.9\times {10}^{-5},3.2\times {10}^{-5}]$$h^{-1}$	${\mu }_{7,2}=[2.7\times {10}^{-1},3.1\times {10}^{-1}]$$h^{-1}$
H	${\lambda }_8=[4.1\times {10}^{-5},5\times {10}^{-5}]$$h^{-1}$	${\mu }_8=[2.7\times {10}^{-1},3\times {10}^{-1}]$$h^{-1}$

(the parameters ${\mu }_{(k,i)}$ and ${\lambda }_{(k,i)}$ are used to represent the repair rate i and failure rate i, respectively, of component k).

Since we assume that the components are independent, we can construct the Markov graph for each component (Figure 9) in order to determine its interval transition matrix. We then apply Equation (1) to find ${\Pi }_{\left(k\right)}$ for each component, where $k=A,\dots ,H$. Finally, we compute the availability of each component by computing the sum of ${\Pi }_{(k,i)}$ over all the working states. In this example, the availability of the entire system $A_{system}$ is presented in the following equation:

$$A_{system}=1-(1-A_A{A}_B{A}_C{A}_D)(1-A_E{A}_F{A}_G{A}_H)$$

(51)

We apply the FBP method to compute ${\Pi }_{\left(k\right)}$ of components. Then, we compare the system availability obtained using this method with the availabilities obtained using the exact method and the precise method. To simplify the comparison, we elected to present the results in terms of the unavailability ($1-availability$). All the results are presented in

Table 9. The unavailability of the system using different methods.

Exact Method	FBP Contraction Method	Precise Method
$[3.318\times {10}^{-8},2.987\times {10}^{-7}]$	$[8.508\times {10}^{-9},5.310\times {10}^{-7}]$	$8.346\times {10}^{-8}$

and in Figure 10.

In this example, the obtained interval of system availability using the FBP method contains the interval obtained using the exact method. The GI metric is also 1. When we calculate the CR for this case, we have

$$C{R}_2=\frac{5.224\times {10}^{-7}-2.655\times {10}^{-7}}{5.224\times {10}^{-7}}=0.491$$

(52)

Table 9 shows that the interval obtained using the FBP contraction technique contains the interval obtained using the exact method, and the precise availability belongs to the two availability intervals. As we can see, the FBP method is an efficient methods since it leads to accurate results, even when dealing with larger systems. However, $C{R}_2>C{R}_1$ since the second system is larger than the first and it has a larger number of components and a larger number of states.

7. Conclusions

In this work, we proposed an original methodology for computing the imprecise availability of multi-state systems (MSS) based on the Markovian model and interval contraction techniques. The imprecision of data was represented by intervals of repair/failure rates. The FBP algorithm was chosen for interval contractions. Regarding the previous work in the literature, this methodology provides many advantages:

• It can be considered a generalization of classical availability theory for MSS, and thus the results obtained when considering only the center of repair/failure rates intervals are special cases of our proposed method.
• The availability results are guaranteed in the sense that we are sure that the obtained availability intervals contain the exact value.
• The imprecision of the availability results (the length of availability intervals) reflects the incompleteness of initial transitions rates and can be used to assist the decision maker by collecting, for example, more data if it is possible.
• It is less time consuming comparing to the exact method, which consists of considering all the possibilities (upper and lower bounds) of transition rates. On the other hand, it is always more conservative than the exact method.
• It is straightforward regarding the BUGF and IUGF methods in the sense that we do not have to transform the repair/failure rates into probabilities.
• It can be applied to all types of configurations, whereas the BUGF and IUGF methods can only be used in series parallel and parallel series configurations.

Moreover, two case studies were proposed. The first one helped us in illustrating the use of the methodology. It also indicates that for small systems, our proposed method gives results very similar to those obtained by the exact method. The second one was proposed to show the computational efficiency of our proposed method for large systems. Moreover, some other techniques can be used to further improve the time computation of our method. For example, the reduction in the state space of components can be performed by finding a smaller component with equivalent behavior and substituting the original component by the smaller one in the Markov model. Another idea for the model reduction is to aggregate states into super states to obtain a Markov model on the aggregated state space.