V Kreinovich

C Jacob

D Dubois
email: dubois@irit.fr

J Cardoso
email: cardoso@isae.fr

M Ceberio
email: mceberio@utep.edu

I Batyrshin

FAILURE ANALYSIS OF A COMPLEX SYSTEM BASED ON PARTIAL INFORMATION ABOUT SUBSYSTEMS, WITH POTENTIAL APPLICATIONS TO AIRCRAFT MAINTENANCE

Keywords: complex system, probability of failure, interval uncertainty AMS Subject Classification: 90B25, 68T37, 65G40, 60K10, 62N05

In many real-life applications (e.g., in aircraft maintenance), we need to estimate the probability of failure of a complex system (such as an aircraft as a whole or one of its subsystems). Complex systems are usually built with redundancy allowing them to withstand the failure of a small number of components. In this paper, we assume that we know the structure of the system, and, as a result, for each possible set of failed components, we can tell whether this set will lead to a system failure. For each component A, we know the probability P (A) of its failure with some uncertainty: e.g., we know the lower and upper bounds P (A) and P (A) for this probability. Usually, it is assumed that failures of different components are independent events. Our objective is to use all this information to estimate the probability of failure of the entire the complex system. In this paper, we describe several methods for solving this problem, including a new efficient method for such estimation based on Cauchy deviates.

Formulation of the Problem

It is necessary to estimate the probability of failure for complex systems. In many practical applications, we need to estimate the probability of failure of a complex system. The need for such estimates comes from the fact that in practice, while it is desirable to minimize risk, it is not possible to completely eliminate it: no matter how many precautions we take, there are always some very low probability events that can potentially lead to a system's failure. All we can do is to make sure that the resulting probability of failure does not exceed the desired small value p 0 . For example, the probability of a catastrophic event is usually required to be at or below p 0 = 10 -9 .

In aircraft design and maintenance, we need to estimate the probability of a failure of an aircraft as a whole and of its subsystems. At the design stage, the purpose of this estimate is to make sure that this probability of failure does not exceed the allowed probability p 0 . At the maintenance stage, this estimate helps to decide whether a maintenance is needed: if the probability of failure exceeds p 0 , some maintenance is required to bring this probability down to the desired level p 0 (or below).

Information available for estimating system's probability of failure: general description. Complex systems consist of subsystems, which, in turn, consist of components (or maybe of subsubsystems which consist of components). So, to estimate the probability of failure of a complex system, we need to take into account when the failure of components and subsystems lead to the failure of the complex system as a whole, and how reliable are these components and subsystems.

From the failure of components and subsystems to the failure of the complex system as a whole. Complex systems are usually built with redundancy allowing them to withstand the failure of a small number of components. Usually, we know the structure of the system, and, as a result, for each possible set of failed components, we can tell whether this set will lead to a system failure. So, in this paper, we will assume that this information is available.

How reliable are components and subsystems? What do we know about the reliability of individual components? For each component A, there is a probability P (A) of its failure. When we have a sufficient statistics of failures of this type of components, we can estimate this probability as the relative frequency of cases when the component failed. Sometimes, we have a large number of such cases, and as a result, the frequency provides a good approximation to the desired probability -so that, in practice, we can safely assume that we know the actual values of these probabilities P (A).

If only a few failure cases are available, it is not possible to get an accurate estimate for P (A). In this case, the only information that we can extract from the observation is the interval P(A) = [P (A), P (A)] that contains the actual (unknown) value of this probability.

This situation is rather typical for aircraft design and maintenance, because aircrafts are usually built of highly reliable components -at least the important parts of the aircraft are built of such components -and there are thus very few observed cases of failure of these components.

Fuzzy case. In some cases, especially on the design stage, we do not yet have failure statistics, so we have to rely on expert estimates, estimates based on the expert's experience with similar components. These estimates are usually formulated by using words from natural language, like "about 1%". A natural way to describe such expert estimates is to use fuzzy techniques (see, e.g., [START_REF] Klir | Fuzzy Sets and Fuzzy Logic[END_REF][START_REF] Nguyen | A First Course in Fuzzy Logic[END_REF]), i.e., to describe each such estimate as a fuzzy number P(A). A fuzzy number means, crudely speaking, that for different degrees of uncertainty α, we have an interval P(A, α) that contains the actual (unknown) probability P (A) with this degree of uncertainty: e.g., the interval P(A, 0) contains P (A) with guarantee (uncertainty 0), while the interval P(A, 0.5) contains P (A) with uncertainty 0.5.

Are the component failures independent events or they are caused by a common cause? In many practical situations, failures of different components are caused by independent factors. For example, for an aircraft, possible failures of mechanical subsystems can be caused by the material fatigue, while possible failures of electronic systems can be caused by the interference of atmospheric electricity (e.g., when flying close to a thunderstorm) -and material fatigue and atmospheric electricity are independent events. As a result, usually, it is assumed that failures of different components are independent events. However, sometimes, we know that the failures of different components are caused by a common cause. In this case, the failures of different components are no longer independent events. In such situations, often, do not have any specific information about the correlation between these failures. This is a typical situation in aircraft design and maintenance -we have very few failures of each component, not enough to determine the exact probability of each such failure, and we have even fewer situations when two components failed, so an empirical determination of such correlations is out of question.

What we do in this paper. Our objective is to use all this information to estimate the probability of failure of the entire complex system. In this paper, we describe new methods for such estimation; some preliminary results first appeared in [START_REF] Jacob | Estimating probability of failure of a complex system based on partial information about subsystems and components, with potential applications to aircraft maintenance[END_REF][START_REF] Kreinovich | Estimating Probability of Failure of a Complex System Based on Inexact Information about Subsystems and Components, with Potential Applications to Aircraft Maintenance[END_REF].

Comment. In most of this paper, we make a usual assumption that failures of different components are independent events. Sometimes, we know that the failures of different components are caused by a common cause; corresponding algorithms are described, e.g., in [START_REF] Ceberio | Interval-type and affine arithmetic-type techniques for handling uncertainty in expert systems[END_REF][START_REF] Chopra | Affine arithmetic-type techniques for handling uncertainty in expert systems[END_REF][START_REF] Chopra | Affine arithmetic-type techniques for handling uncertainty in expert systems[END_REF][START_REF] Jacob | Estimating probability of failure of a complex system based on partial information about subsystems and components, with potential applications to aircraft maintenance[END_REF].

2. Simplest Case: Component Failures are Independent and Failure Probabilities P (A) Are Exactly Known

Let us start our analysis with the simplest case when the component failures are independent and the failure probabilities P (A) for different components A are known exactly. As we mentioned, we assume that there exist efficient algorithms that, given a list of failed components, determines whether the whole system fails or not. In this case, it is always possible to efficiently estimate the probability P of the system's failure by using Monte-Carlo simulations. Specifically, we select the number of simulations N . Then, for each component A, we simulate a Boolean variable failing(A) which is true with probability P (A) and false with the remaining probability 1 -P (A). This can be done, e.g., if we take the result r of a standard random number generator that generates values uniformly distributed on the interval [0, 1] and select failing(A) to be true if r ≤ P (A) and false otherwise: then the probability of this variable to be true is exactly P (A).

Then, we apply the above-mentioned algorithm to the simulated values of the variables failing(A) and conclude whether for this simulation, the system fails or not. As an estimate for the probability of the system's failure, we then take the ratio p def = f /N , where f is the number of simulations on which the system failed. From statistics, it is known that the mean value of this ratio is indeed the desired probability, that the standard deviation can be estimated as

σ = p • (1 -p)/N ≤ 0.5/
√ N , and that for sufficiently large N (due to the Central Limit Theorem), the distribution of the difference P -p is close to normal. Thus, with probability 99.9%, the actual value P is within the three-sigma interval [p -3σ, p + 3σ].

This enables us to determine how many iterations we need to estimate the probability P with accuracy 10% (and certainty 99.9%): due to σ ≤ 0.5/ √ N , to guarantee that 3σ ≤ 0.1, it is sufficient to select N for which 3 • 0.5/ √ N ≤ 0.1, i.e., √ N ≥ (3 • 0.5)/0.1 = 15 and N ≥ 225. It is important to emphasize that this number of iterations is the same no matter how many components we have -and for complex systems, we usually have many thousands of components.

Similarly, to estimate this probability with accuracy 1%, we need N = 22, 500 iterations, etc. These numbers of iterations work for all possible values P . In practical applications, the desired probability P is small, so 1 -P ≈ 1, σ ≈ P/N and the number of iterations, as determined by the condition 3σ ≤ 0.1 or 3σ ≤ 0.01, is much smaller: N ≥ 900 • P for accuracy 10% and N ≥ 90, 000 • P for accuracy 1%.

Comment. In many cases, there are also efficient analytical algorithms for computing the desired probability of the system's failure; see, e.g., [START_REF] Dutuit | Approximate estimation of system reliability via fault trees[END_REF][START_REF] Flage | Handling of epistemic uncertainties in fault tree analysis: a comparison between probabilistic, possibilistic, and hybrid approaches[END_REF][START_REF] Guth | A probability foundation for vagueness and imprecision in fault tree analysis[END_REF][START_REF] Troffaes | On the use of the imprecise Dirichlet model with fault trees[END_REF].

Important Subcase of the Simplest Case: When Components are Very Reliable

In many practical applications (e.g., in important subsystems related to aircrafts), components are highly reliable, and their probabilities of failure P (A) are very small. In this case, the above Monte-Carlo technique for computing the probability P of the system's failure requires a large number of simulations, because otherwise, with high probability, in all simulations, all the components will be simulated as working properly.

For example, if the probability of a component's failure is P (A) = 10 -3 , then we need at least a thousand iteration to catch a case when this component fails; if P (A) = 10 -6 , we need at least a million iterations, etc.

In such situations, Monte-Carlo simulations may take a lot of computation time. In some applications, e.g., on the stage of an aircraft design, it may be OK, but in other cases, e.g., on the stage of routine aircraft maintenance, the airlines want fast turnaround, and any speed up is highly welcome.

To speed up such simulations, we can use a re-scaling idea; see, e.g., [START_REF] Jacob | Estimating probability of failure of a complex system based on partial information about subsystems and components, with potential applications to aircraft maintenance[END_REF][START_REF] Jaksurat | Probabilistic approach to trust: ideas, algorithms, and simulations[END_REF]. Specifically, instead of using the original values P (A), we use re-scaled (larger) values λ • P (A) for some λ 1. The value λ is chosen in such a way that the resulting probabilities are larger and thus, require fewer simulations to come up with cases when some components fail. As a result of applying the above Monte-Carlo simulations to these new probabilities λ • P (A), we get a probability of failure P (λ).

In this case, one can show that while the resulting probabilities λ • P (A) are still small, the probability P (λ) depends on λ as P (λ) ≈ λ k • P for some positive integer k.

Thus, to find the desired value P , we repeat this procedure for two different values λ 1 = λ 2 , get the two values P (λ 1 ) and P (λ 2 ), and then find both unknown k and P from the resulting system of two equations with two unknowns: P (λ 1 ) ≈ λ k 1 • P and P (λ 2 ) ≈ λ k 2 • P . To solve this system, we first divide the first equation by the second one, getting an equation P (λ 1 )/P (λ 2 ) ≈ (λ 1 /λ 2 ) k with one unknown k, and find k ≈ ln(P (λ 1 )/P (λ 2 ))/(λ 1 /λ 2 ). Then, once we know k, we can find P as P ≈ P (λ 1 )/λ k 1 .

4.

Cases When We Know the Probabilities P (A) with Uncertainty: Fuzzy Uncertainty Can Be Reduced to Interval One

In many cases, we do not know the exact probabilities P (A) of the component's failure, we only know the intervals P(A) that contain these probabilities, or, even more generally, a fuzzy number P(A) that describes this probability.

In the case of intervals, different combinations of values P (A) ∈ P(A) lead, in general, to different values P . Let us denote the dependence of P on the values P (A) by f : P = f (P (A), P (B), . . .). In this case, we want to know the range of possible values of the desired probability P : The problem of computing such an interval is a particular case of the general problem of interval computations, i.e., a problem of computing the range of a given function f (x 1 , . . . , x n ) when each of the variables x i takes value from a given interval x i ; see, e.g., [7,[START_REF] Jaulin | Applied Interval Analysis: With Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF][START_REF] Moore | Introduction to Interval Analysis[END_REF] and references therein.

When each probability is described by a fuzzy number P(A), i.e., by a membership function µ A (P ) that assigns, to every real number P , a degree to which this number is possible as a value of P (A), we want to find the fuzzy number P that describes f (P (A), P (B), . . .). A natural way to define the corresponding membership function µ(P ) leads to Zadeh's extension principle:

µ(P ) = sup{min(µ A (P (A)), µ B (P (B)), . . .) : f (P (A), P (B), . . .) = P }.
It is known that from the computational viewpoint, the application of this formula can be reduced to interval computations.

Specifically, for each fuzzy set with a membership function µ(x) and for each α ∈ (0, 1], we can define this set's α-cut as X (α)

def = {x : µ(x) ≥ α}.
Vice versa, if we know the α-cuts for all α, we, for each x, can reconstruct the value µ(x) as the largest value α for which x ∈ X (α). Thus, to describe a fuzzy number, it is sufficient to find all its α-cuts.

It is known that when the inputs µ i (x i ) are fuzzy numbers, and the function y = f (x 1 , . . . , x n ) is continuous, then for each α, the α-cut Y(α) of y is equal to the range of possible values of f (x 1 , . . . , x n ) when x i ∈ X i (α) for all i:

Y(α) = f (X 1 (α), . . . , X n (α));
see, e.g., [?, 12, ?, 17]. This is how processing fuzzy data is often done -by reducing to interval computations. P = [P , P ] = f (P(A), P(B), . . .). In particular, for our problem, once know the α-cuts P(A, α) corresponding to different components A, we can find the α-cuts P(α) corresponding to the desired probability P as the corresponding interval range:

P(α) = f (P(A, α), P(B, α), . . .).
So, if we know how to solve our problem under interval uncertainty, we can also solve it under fuzzy uncertainty -e.g., by repeating the above interval computations for α = 0, 0.1, . . . , 0.9, 1.0.

In view of this reduction, in the following text, we will only consider the case of interval uncertainty.

5. Component Failures are Independent, Failure Probabilities P (A) Are Known with Interval Uncertainty: Monotonicity Case

In view of the above analysis, let us now consider the case when the probabilities P (A) are only known with interval uncertainty, i.e., for each component A, we only know the interval [P (A), P (A)] that contains the actual (unknown) value P (A). We still assume that failures of different components are independent events.

Let us start with the simplest subcase when the dependence of the system's failure is monotonic with respect to the failure of components. To be precise, we assume that if for a certain list of failed components, the system fails, it will still fail if we add one more components to the list of failed ones. In this case, the smaller the probability of failure P (A) for each component A, the smaller the probability P that the system as a whole will fail. Similarly, the larger the probability of failure P (A) for each component A, the larger the probability P that the system as a whole will fail.

Thus, to compute the smallest possible value P of the failure probability, it is sufficient to consider the values P (A). Similarly, to compute the largest possible value P of the failure probability, it is sufficient to consider the values P (A). Thus, in the monotonic case, to compute the range [P , P ] of possible values of overall failure probability under interval uncertainty, it is sufficient to solve two problems in each of which we know probabilities with certainty:

• to compute P , we assume that for each component A, the failure probability is equal to P (A); • to compute P , we assume that for each component A, the failure probability is equal to P (A).

In Practice, the Dependence is Sometimes Non-Monotonic

In some practically reasonable situations, the dependence of the system's failure on the failure of components is non-monotonic; see, e.g., [START_REF] Jacob | Estimating probability of failure of a complex system based on partial information about subsystems and components, with potential applications to aircraft maintenance[END_REF]. This may sound counter-intuitive at first glance: adding one more failing component to the list of failed ones suddenly makes the previously failing system recover, but here is an example when exactly this seemingly counter-intuitive behavior makes perfect sense. Please note that this example is over-simplified: its only purpose is to explain, in intuitive terms, the need to consider non-monotonic case.

To increase reliability, systems include duplication: for many important functions, there is a duplicate subsystem ready to take charge if the main subsystem fails. How do we detect that the main system failed? Usually, a subsystem contains several sensors; sensors sometimes fail, as a result of which their signal no longer reflect the actual value of the quantity they are supposed to measure. For example, a temperature sensor which is supposed to generate a signal proportional to the temperature, if failed, produces no signal at all, which the system will naturally interpret as a 0 temperature. To detect the sensor failure, subsystems often use statistical criteria. For example, for each sensor i, we usually know the mean m i and the standard deviation σ i of the corresponding quantity. When these quantities are independent and approximately normally distributed, then, for the measurement values x i , the sum

X 2 def = n i=1 (x i -m i ) 2 σ 2 i
is the sum of n standard normal distributions and thus, follows the chi-square distributed with n degrees of freedom. So, if the actual value of this sum exceeds the threshold corresponding to confidence level p = 0.05, this means that we can confidently conclude that some of the sensors are malfunctioning. If the number n of sensors is large, then one malfunctioning sensor may not increase the sum X 2 too high, and so, its malfunctioning will not be detected, and the system will fail.

On the other hand, if all n sensors fail, e.g., show 0 instead of the correct temperature, each term in the sum will be large, the sum will exceed the threshold -and the system will detect the malfunctioning. In this case, the second redundant subsystem will be activated, and the system as a whole will thus continue to function normally. This is exactly the case of non-monotonicity: when only one sensor fails, the system as a whole fails; however, if, in addition to the originally failed sensor, many other sensors fail, the system as a whole becomes functioning well. Other examples of non-monotonicity can be due to the fact that some components may be in more than two states [START_REF] Jacob | Uncertainty handling in quantitative BDD-based fault-tree analysis by interval computation[END_REF].

In the following text, we will consider the non-monotonic case, in which a simple algorithm (given above) is not applicable.

What If We Have Few Components, With Respect to Which the Dependence is Non-Monotonic

Let us start with the simplest case when there are only few components with respect to which the dependence is non-monotonic, and with respect to all other components, the dependence is still monotonic. In this case, for all monotonic components B, as before, we can take P (B) = P (B) when we are computing P , and take P (B) = P (B) when we are computing P .

For non-monotonic components A, for computing each of the values P and P , we need to to take into account all possible values P (A) ∈ [P (A), P (A)]. For each such component, by using the formula of full probability, we can represent the probability P of the system's failure as follows:

P = P (A) • P (F |A) + (1 -P (A)) • P (F |¬A),
where P (F |A) is the conditional probability that the system fails under the condition that the component A fails, and P (F |¬A) is the conditional probability that the system fails under the condition that the component A does not fail. The conditional probabilities P (F |A) and P (F |¬A) do not depend on P (A), so the resulting dependence of P on P (A) is linear. A linear function attains it minimum and maximum at the endpoints. Thus, to find P and P , it is not necessary to consider all possible values P (A) ∈ [P (A), P (A)], it is sufficient to only consider two values: P (A) = P (A) and P (A) = P (A).

For each of these two values, for another non-monotonic component A , we have two possible options P (A ) = P (A ) and P (A ) = P (A ); thus, in this case, we need to consider 2 × 2 = 4 possible combinations of values P (A) and P (A ).

In general, when we have k non-monotonic components A 1 , . . . , A k , it is sufficient to consider 2 k possible combinations of values P (A i ) and P (A i ) corresponding to each of these components. When k is small, this is doable -and, as our preliminary experiments show, works very well.

This procedure requires times which grows as 2 k . As we mentioned earlier, when k is large, the needed computation time becomes unrealistically large.

General Case: the Problem is NP-Hard, and Even Checking Whether a

Component Is Monotonic Is NP-hard Natural question. The fact that the above algorithm requires unrealistic exponential time raises a natural question: is it because our algorithm is inefficient or is it because the problem itself is difficult?

The problem is NP-hard. In the general case, when no assumption is made about monotonicity, the problem is as follows:

• we have a propositional formula F with n variables A i -each variable A i is true if the corresponding component fails, and the formula F is true if the system as whole fails; • for each component i, we know the interval [P (A i ), P (A i )] that contains the actual (unknown) P (A i ) that this component fails; • we assume that the failures of different components are independent events. Different values P (A i ) ∈ [P (A i ), P (A i )] lead, in general, to different values of the probability P that F is true (i.e., that the system failed). Our objective is to compute the range [P , P ] of possible values of this probability.

One can easily show that, in general, this problem is NP-hard (for precise definitions, see, e.g., [START_REF] Kreinovich | Computational Complexity and Feasibility of Data Processing and Interval Computations[END_REF][START_REF] Papadimitriou | Computational Complexity[END_REF]). Indeed, it is well known that the following propositional satisfiability problem SAT is NP-hard: given a propositional formula F , check whether this formula is satisfiable, i.e., whether there exist values A i that make it true. We will prove that our problem is NP-hard by reducing SAT to it: for every particular case F of SAT, there is a particular case of our problem whose solution leads to a solution to the F . Indeed, for every formula F , let us take [P (A i ), P (A i )] = [0, 1] for all variables i.

If the formula F is not satisfiable, then F is always false, so the probability of its being true is 0. In this case, the range of possible values of P consists of a single value 0: [P , P ] = [0, 0].

On the other hand, if F is satisfiable, e.g., if F is true for A 1 true, A 2 false, etc., then we can take P (A 1 ) = 1, P (A 2 ) = 0, etc., and conclude that under these probabilities, F is always true, i.e., P = 1. Thus, in this case, P = 1.

So, by computing P , we will get either P = 0 or P = 1. In the first case, F is satisfiable, in the second case, it is not. The reduction proves that our problem is NP-hard.

Even checking monotonicity is NP-hard. Let us show that even checking monotonicity is NPhard. Intuitively, monotonicity with respect to a component A means that if the system was in a failing state, and we change the state of A to failing, the system remains failing. In precise terms, monotonicity means that if the formula F was true, and we change the value of the variable A from false to true, then F remains true.

Let us prove that the problem of checking monotonicity is NP-hard by reducing SAT to this problem. Indeed, for every propositional formula F (A 1 , . . . , A n ), we can form a new formula

F def = F (A 1 , . . . , A n ) &¬A 0 .
When the original formula F is not satisfiable, then F is always false and thus, the new formula F is also always false. In this case, F is, by definition, monotonic with respect to A 0 .

When F is satisfiable, this means that F is equal to "true" for some values of A 1 , . . . , A n . In this case, for A 0 ="false", the new formula F is true, but if we make A 0 ="true", F becomes false. Thus, in this case, F is not monotonic with respect to A 0 .

Summarizing: the new formula F is monotonic with respect to A 0 if and only if the original formula F was satisfiable. This reduction proves that checking monotonicity is indeed NP-hard.

A Practically Important Case when Dependence May Be Non-Monotonic but

Intervals Are Narrow: Towards a New Algorithm General non-monotonic case: a possible algorithm. For each component A, by using the formula of full probability, we can represent the probability P of the system's failure as follows:

P = P (A) • P (F |A) + (1 -P (A)) • P (F |¬A),
where P (F |A) is the conditional probability that the system fails under the condition that the component A fails, and P (F |¬A) is the conditional probability that the system fails under the condition that the component A does not fail. The conditional probabilities P (F |A) and P (F |¬A) do not depend on P (A), so the resulting dependence of P on P (A) is linear. A linear function attains it minimum and maximum at the endpoints. Thus, to find P and P , it is not necessary to consider all possible values P (A) ∈ [P (A), P (A)], it is sufficient to only consider two values: P (A) = P (A) and P (A) = P (A).

For each of these two values, for another component A , we have two possible options P (A ) = P (A ) and P (A ) = P (A ); thus, in this case, we need to consider 2×2 = 4 possible combinations of values P (A) and P (A ).

In general, when we have k components A 1 , . . . , A k , it is sufficient to consider 2 k possible combinations of values P (A i ) and P (A i ) corresponding to each of these components. This procedure requires times which grows as 2 k . As we mentioned earlier, when k is large, the needed computation time becomes unrealistically large. Natural question. The fact that the above algorithm requires unrealistic exponential time raises a natural question: is it because our algorithm is inefficient or is it because the problem itself is difficult?

The problem is NP-hard. In the general case, when no assumption is made about monotonicity, the problem is as follows:

• Let F be a propositional formula with n variables A i • for each variable A i , we know the interval [P (A i ), P (A i )] that contains the actual (unknown) P (A i ) that this variable is true; • we assume that the Boolean variables are independent. Different values P (A i ) ∈ [P (A i ), P (A i )] lead, in general, to different values of the probability P that F is true (e.g., that the system fails). Our objective is to compute the range [P , P ] of possible values of this probability.

In [START_REF] Jacob | Estimating probability of failure of a complex system based on partial information about subsystems and components, with potential applications to aircraft maintenance[END_REF], we have proven that, in general, the problem of computing the desired range [P , P ] is NP-hard. From the practical viewpoint, this means, that (unless P=NP, which most computer scientists believe to be not true), there is no hope to avoid non-feasible exponential time. Since we cannot have a feasible algorithm that is applicable to all possible cases of the general problem, we therefore need to restrict ourselves to practically important cases -and try to design efficient algorithms that work for these cases. This is what we do in this paper. A practically important case of narrow intervals. When there is enough information, the intervals [P (A), P (A)] are narrow. If we represent them in the form [ P (A) -∆(A), P (A) + ∆(A)], with P (A) = P (A) + P (A) 2 and ∆(A) = P (A) -P (A) 2 , then values ∆(A) are small, so we can safely ignore terms which are quadratic or of higher order in terms of ∆P (A).

Linearization: analysis of the problem. In the case of narrow intervals, the difference ∆P (A) def = P (A)-P (A) is bounded by ∆(A) and thus, also small: |∆P (A)| ≤ ∆(A). Hence, we can expand the dependence of the desired system failure probability P = P (P (A), . . .) = P ( P (A) + ∆P (A), . . .) into Taylor series and keep only terms which are linear in ∆P (A):

P ≈ P + A c A • ∆P (A)

.).

For those A for which c A ≥ 0, the largest value of the sum

A c A • ∆P (A) (when ∆P (A) ∈ [-∆(A), ∆(A)]
) is attained when ∆P (A) attains its largest possible value ∆(A). Similarly, when c A < 0, the largest possible values of the sum is attained when ∆P (A) = -∆(A). In both cases, the largest possible value of the term

c A • ∆P (A) is |c A | • ∆(A).
Thus, the largest possible value of P is equal to P + ∆, where

∆ def = A |c A | • ∆(A).
Similarly, one can show that the smallest possible value of P is equal to P -∆, so the range of possible values of the failure probability P is [ P -∆, P + ∆].

We already know how to compute P -e.g., we can use the Monte-Carlo approach. How can we compute ∆? How to compute ∆: numerical differentiation and its limitations. A natural idea is to compute all the partial derivatives c A and to use the above formula for ∆. By definition, c A is the derivative, i.e., c A = lim h→0 P ( P (A) + h, P (B), P (C), . . .) -P ( P (A), P (B), P (C), . . .) h .

By definition of the limit, this means that to get a good approximation for c A , we can take a small h and compute c A = P ( P (A) + h, P (B), P (C), . . .) -P ( P (A), P (B), P (C), . . .) h .

This approach to computing derivatives is called numerical differentiation.

The problem with this approach is that each computation of the value P ( P (A) + h, P (B), P (C), . . .)

by Monte-Carlo techniques requires a lot of simulations, and we need to repeat these simulations again and again as many times as there are components. For an aircraft, with thousands of components, the resulting increase in computation time is huge. Moreover, since we are interested in the difference P ( P (A) + h, . . .) -P ( P (A), . . .) between the two probabilities, we need to compute each of these probabilities with a high accuracy, so that this difference would be visible in comparison with the approximation error ∼ 1/ √ N of the Monte-Carlo estimates. This requires that we further increase the number of iterations N in each Monte-Carlo simulation and thus, even further increase the computation time.

Cauchy deviate techniques: reminder. In order to compute the value

A |c A | • ∆(A)
faster, one may use a technique based on Cauchy distributions (e.g., [START_REF] Kreinovich | A new Cauchy-based black-box technique for uncertainty in risk analysis[END_REF][START_REF] Trejo | Error estimations for indirect measurements: randomized vs. deterministic algorithms for 'black-box' programs[END_REF]) , i.e., probability distributions with probability density of the form ρ

(z) = ∆ π • (z 2 + ∆ 2 )
; the value ∆ is called the scale parameter of this distribution, or simply a parameter, for short. Cauchy distribution has the following property: if z A corresponding to different A are independent random variables, and each z A is distributed according to the Cauchy law with parameter ∆(A), then their linear combination z = A c A • z A is also distributed according to a Cauchy law, with a scale parameter ∆ =

A |c A | • ∆(A).
Therefore, using Cauchy distributed random variables δ A with parameters ∆(A), the difference c def = P ( P (A) + δ A , P (B) + δ B , . . .) -P ( P (A), P (B), . . .)

= A c A • δ A
is Cauchy distributed with the desired parameter ∆. So, repeating this experiment N c times, we get N c values c (1) , . . . , c (Nc) which are Cauchy distributed with the unknown parameter, and from them we can estimate ∆. The bigger N c , the better estimates we get. Comment. To avoid confusion, we should emphasize that the use of Cauchy distributions is a computational technique, not an assumption about the actual distribution: indeed, we know that the actual value of ∆P (A) is bounded by ∆(A), but for a Cauchy distribution, there is a positive probability that the simulated value is larger than ∆(A). Cauchy techniques: towards implementation. In order to implement the above idea, we need to answer the following two questions:

• how to simulate the Cauchy distribution;

• how to estimate the parameter ∆ of this distribution from a finite sample. Simulation can be based on the functional transformation of uniformly distributed sample values:

δ A = ∆(A) • tan(π • (r A -0.5))
, where r A is uniformly distributed on the interval [0, 1].

In order to estimate ∆, we can apply the Maximum Likelihood Method

ρ(c (1) ) • ρ(c (2) ) • . . . • ρ(c (Nc) ) → max,
where ρ(z) is a Cauchy distribution density with the unknown ∆. When we substitute the above-given formula for ρ(z) and equate the derivative of the product with respect to ∆ to 0 (since it is a maximum), we get an equation

1 1 + c (1) ∆ 2 + . . . + 1 1 + c (Nc) ∆ 2 = N c 2 .
Its left-hand side is an increasing function that is equal to 0(< N c /2) for ∆ = 0 and > N c /2 for ∆ = max c (k) ; therefore the solution to this equation can be found by applying a bisection method to the interval 0, max c (k) .

It is important to mention that we assumed that the function P is reasonably linear when the values δ A are small: |δ A | ≤ ∆(A). However, the simulated values δ A may be larger than ∆(A). When we get such values, we do not use the original function P for them, we use a normalized function that is equal to P within the given intervals, and that is extended linearly for all other values; we will see, in the description of an algorithm, how this is done. Cauchy deviate technique: main algorithm.

• Apply P to the values P (A) and compute P = P ( P (A), P (B), . . .).

• For k = 1, 2, . . . , N c , repeat the following:

• use the standard random number generator to compute n numbers r

A that are uniformly distributed on the interval [0, 1];

• compute Cauchy distributed values c (k) A = tan(π • (r (k) A -0.5)); • compute the largest value of |c (k)
A | so that we will be able to normalize the simulated measurement errors and apply P to the values that are within the box of possible values:

K = max A |c (k) A |; • compute the simulated measurement errors δ (k) A := ∆(A) • c (k) A /K; • compute the simulated probabilities P (k) (A) = P (A) + δ (k) A ;
• estimate P (P (k) (A), P (k) (B), . . .) and then compute c (k) = K • (P (P (k) (A), P (k) (B), . . .) -P );

• Compute ∆ by applying the bisection method to solve the corresponding equation.

Resulting gain and remaining limitation. By using the Monte-Carlo techniques, we make sure that the number of iterations N c depends only on the accuracy with which we want to find the result and not on the number of components. Thus, when we have a large number of components, this method is faster than numerical differentiation.

The computation time of the new algorithm is smaller, but it is still not very fast. The reason is that the Cauchy method was originally was designed for situations in which we can compute the exact value of P (P (k) (A), P (k) (B), . . .). In our problem, these values have to be computed by using Monte-Carlo techniques, and computed accurately -and each such computation requires a lot of iterations. Instead of running the maximum likelihood, we can also just estimate ∆ by means of the sample interquartile range instead of solving the non-linear equation. But this method will be less accurate. Final idea to further decrease the needed number of simulations. (see, e.g., Section 5.4 of [START_REF] Trejo | Error estimations for indirect measurements: randomized vs. deterministic algorithms for 'black-box' programs[END_REF]) For each combination of values δ A , the corresponding Monte-Carlo simulation produces not the actual probability P ( P (A) + δ A , P (B) + δ B , . . .), but an approximate value P ( P (A) + δ A , P (B) + δ B , . . .) = P ( P (A) + δ A , P (B) + δ B , . . .) + c n that differs from the desired probability by a random variable c n which is normally distributed with mean 0 and variance σ The components c c and c n are independent. Thus, for c = c c + c n , for the characteristic function χ(ω)

def = E[exp(i • ω • c)], we have E[exp(i • ω • c)] = E[exp(i • ω • c c ) • exp(i • ω • c n )] = E[exp(i • ω • c c )] • E[exp(i • ω • c n )],
i.e., χ(ω) = χ c (ω) • χ n (ω), where χ c (ω) and χ n (ω) are characteristic functions of c c and c n . For Cauchy distribution and for the normal distribution, the characteristic functions are known: χ c (ω) = exp(-|ω| • ∆) and χ n (ω) = exp(-ω 2 • σ 2 ). So, we conclude that χ(ω) = exp(-|ω| • ∆ω 2 • σ 2 ). Hence, to determine ∆, we can estimate χ(ω), compute its negative logarithm, and then compute ∆ (see the formula below).

Since the value χ(ω) is real, it is sufficient to consider only the real part cos(. . .) of the complex exponent exp(i • . . .). Thus, we arrive at the following algorithm: Algorithm. First, we use a lengthy Monte-Carlo simulation to compute the value P = P ( P (A), P (B), . . .). Then, for k = 1, 2, . . . , N , we repeat the following:

• use a random number generator to compute n numbers r (k) A , that are uniformly distributed on the interval [0, 1]; • for a real number ω > 0, compute

• compute δ (k) A = ∆ i • tan(π • (r
χ(ω) = 1 N • N k=1 cos ω • c (k) ; • compute ∆ = - ln(χ(ω)) ω -σ 2 • ω 2 .
Comment. Of course, we also need, as before, to "reduce" the simulated values δ A to the given bounds ∆(A).

What If We Do Not Assume Independence

Exact methods. If we do not assume independence, then, in principle, we can have different probabilities P of failure even when the failure probabilities P (A i ) of all components are known exactly. For example, if the system consists of two components and it fails if both components fail, i.e., if In general, to describe probabilities of all possible combinations of statements A i , it is sufficient to describe 2 n probabilities of atomic statements A ε 1 1 & . . . & A εn n , where ε i ∈ {-, +}, A + means A, and A -means ¬A. These probabilities satisfy the condition that their sum is 1; each given probability P (A i ) and the desired probability P can be described as a sum of some such probabilities, so the problem of finding the range of P becomes the particular case of linear programming problems, when we need to find the minimum and maximum of a linear function under linear constraints; see, e.g., [START_REF] Walley | Statistical reasoning with imprecise probabilities[END_REF].

F = A 1 & A 2 ,
For example, in the above case n = 2, we have four non-negative unknowns

P ++ = P (A 1 & A 2 ), P +-= P (A 1 & ¬A 2 ), P -+ = P (¬A 1 & A 2 )
, and P --= P (A 1 & A 2 ) that satisfy the constraints P ++ + P +-+ P -+ + P --= 1, P ++ + P +-= P (A 1 ), and P ++ + P -+ = P (A 2 ). Here, P = P ++ , so, e.g., to find P , we need to minimize P ++ under these constraints.

Limitations of the exact methods. This works well if n is small, but for large n, this method requires an unrealistically long time.

Heuristic approximate methods. In principle, we can use technique similar to straightforward interval computations. Indeed, for simple formulas F like ¬A 1 , A 1 ∨ A 2 , or A 1 & A 2 , we have explicit formulas for the range of the probability P (F ). So, to estimate the range of the probability P for an arbitrary formula F , we can do the following:

• we parse the expression F , i.e., represent is as a sequence of simple boolean operations -the same sequence that a computer computing F would follow, and • replace each computation step with corresponding operations with probability ranges.

In this algorithm, at each moment of time, we keep the bounds on the probabilities P (A i ) and on the probabilities P (F j ) of the corresponding intermediate formulas.

The problem with this approach is that at each step, we ignore the dependence between the intermediate results F j ; hence intervals grow too wide. For example, the estimate for P (A ∨ ¬A) computed this way is not 1, but an interval containing the correct value 1.

A more accurate algorithm was proposed in [START_REF] Ceberio | Interval-type and affine arithmetic-type techniques for handling uncertainty in expert systems[END_REF][START_REF] Chopra | Affine arithmetic-type techniques for handling uncertainty in expert systems[END_REF][START_REF] Chopra | Affine arithmetic-type techniques for handling uncertainty in expert systems[END_REF]. In this algorithm, at each stage, besides the bounds on the values P (F j ) (including the original bounds on P (A i ), or -if available -exact values of P (A i )), we also compute the bounds for the probabilities P (F j & F k ), P (F j & ¬F k ), P (¬F j & F k ), and P (¬F j & ¬F k ).

On each computation step, when we add a new intermediate result F a (e.g., F a = F b & F c ), we add bounds for the probabilities of the new statement F a and of all possible combinations that include F a , i.e., on the probabilities P (F a & F k ), P (F a & ¬F k ), P (¬F a & F k ), and P (¬F a & ¬F k ). To compute each of these probabilities, we use the known bounds on the probabilities of combinations of F b and F k , F c and F k , F b and F c , and get the desired bounds on the combinations of F a and F k by solving the corresponding linear programming problem. In this linear programming problem, we consider, as variables, 2 3 = 8 probabilities of atomic statements

F ε b b & F εc c & F ε k k .
At the end of the process, we get an interval P . Similarly to interval computations, we can prove, by induction, that every possible value P of the system's failure is contained in this interval, i.e., that the interval P is an enclosure for the desired range [P , P ]: [P , P ] ⊆ P . This algorithm requires more computation time that the above straightforward algorithmsince for s intermediate steps, we now need s 2 estimations instead of s -but as a result, we get more accurate accurate, with smaller "excess width" for the resulting enclosure. For example, the probability P (A ∨ ¬A) is estimated as 1.

To get an even better accuracy, instead of probabilities of all possible combinations of two intermediate results, we can compute, on each step, probabilities of all possible combinations of three such results: F i , F j , and F k . In this case, computation time grows as s 3 but the resulting enclosure is even more accurate. Similarly, we can have combinations of fours, fives, etc.

Conclusion

In this paper we considered the problem of estimating the probability of failure P of a complex system such as an aircraft, assuming we only know upper and lower bounds of probabilities of elementary events such as component failures. The assumptions in this paper is that failures of different components are independent events, and that there is enough information to ensure narrow probability intervals. The problem of finding the resulting range [P , P ] of possible values of P is computationally difficult (NP-hard). In this paper, we describe several efficient algorithms for computing the range. In particular, for the practically important case of narrow intervals [P (A), P (A)], we describe an efficient method that uses Cauchy deviates to estimate the desired range [P , P ]. Future works concern the estimation of intervals [P (A), P (A)] from the imprecise knowledge of failure rates. Moreover, it is interesting to study what can be done in practice when the independence assumption on component failures no longer holds. 12. Acknowledgment V. Kreinovich was supported by the National Science Foundation (NSF) grants HRD-0734825 and DUE-0926721 and by Grant 1 T36 GM078000-01 from the National Institutes of Health (NIH). C. Jacob was supported by a grant from @MOST Prototype, a joint project of Airbus, Institut de Recherche en Informatique deToulouse (IRIT), LAAS-CNRS, ONERA, and ISAE.

We are thankful to the Airbus colleagues for providing us examples of the practical problems, general guidance about which assumptions are reasonable and which are not for aircraft maintenance problems, and for testing our techniques.

We are also thankful to all the participants of the International Workshop on Soft Computing Applications and Knowledge Discovery SCAKD'2011 (Moscow, Russia, June 25, 2011) and the 10th Mexican International Conference on Artificial Intelligence MICAI'2011 (Puebla, Mexico, November 26 -December 4, 2011) for valuable discussions.

P

  = [P , P ] = f (P(A), P(B), . . .) def = {f (P (A), P (B), . . .) : P (A) ∈ P(A), P (B) ∈ P(B), . . .}.

  , where P def = P ( P (A), . . .) and c A def = ∂ ∂P (A) P ( P (A), . .

  2 = P • (1 -P ) N . As a result, the difference c def = P ( P (A)+δ A , P (B)+ δ B , . . .) -P between the two observed probabilities can be represented as c = c c + c n , where c c def = P ( P (A) + δ A , P (B) + δ B , . . .) -P is, as we have mentioned, Cauchy distributed with parameter ∆, while c n = P ( P (A) + δ A , P (B) + δ B , . . .) -P ( P (A) + δ A , P (B) + δ B , . . .) is normally distributed with mean 0 and known standard deviation σ.

A

  -0.5)); • use Monte-Carlo simulations to find the frequency (probability estimate) P ( P (A) + δ (k) A , P (B) + δ (k) B , . . .) and then c (k) = P ( P (A) + δ (k) A , P (B) + δ (k) B , . . .) -P ;

  then possible values of P take an interval [P , P ] = [max(P (A 1 ) + P (A 2 ) -1, 0), min(P (A 1 ), P (A 2 )].

  Vladik Kreinovich -received his M.Sc. in Mathematics and Computer Science from St. Petersburg University, Russia, in 1974, and Ph.D. from the Institute of Mathematics, Soviet Academy of Sciences, Novosibirsk, in 1979. In 1975-80, he worked with the Soviet Academy of Sciences, in particular, in 1978-80, with the Special Astrophysical Observatory (representation and processing of uncertainty in radioastronomy). In 1982-89, he worked on error estimation and intelligent information processing for the National Institute for Electrical Measuring Instruments, Russia. In 1989, he was a Visiting Scholar at Stanford University. Since 1990, he is with the Department of Computer Science, University of Texas at El Paso. Also, served as an invited professor in Paris (University of Paris VI), Hong Kong, St. Petersburg, Russia, and Brazil. Main interests: representation and processing of uncertainty, especially interval computations and intelligent control. Published 3 books, 6 edited books, and more than 800 papers. Member of the editorial board of the international journal Reliable Computing (formerly, Interval Computations), and several other journals. Co-maintainer of the international website on interval computations http://www.cs.utep.edu/interval-comp Christelle Jacob -A graduate student from Université Paul-Sabatier (France) in 2006. She obtained a Master degree in Applied Mathematics from Université Paul-Sabatier (France) in 2008, and a Specialized Master in Embedded Systems from Institut Supérieur de l'Aéronautique et de l'Espace (ISAE, France) in 2009. Presently, she is preparing a Ph.D degree on uncertainty management for preventive maintenance of aircrafts in the project @MOST of Airbus company. Martine Ceberio is an assistant professor in Computer Science at the University of Texas at El Paso (UTEP). She received her PhD in Computer Science in 2003 from the University of Nantes, France and joined UTEP in fall 2003. Her research is in numerical constraint solving, global optimization, multi-criteria decision making, and their applications. She has been engaged in activities related to the participation of women in computing since 2008. She has implemented projects with an all-girls high school of El Paso and she coordinates the local affiliate of a national competition of the National Center for Women in Information and Technology (NCWIT) that recognizes high-school girls with interest and experience in computer science. Ildar Batyrshin -received Diploma of Moscow Physical-Technical Institute in 1975; Ph.D. from Moscow Power Engineering Institute in 1983; Dr.Sci. (Habilitation) from Higher Attestation Committee of Russian Federation in 1996. Since 1975 he has been with the Department of Informatics and Applied Mathematics of Kazan State Technological University, Russia, (Department head in 1997-2003). Since 1999, he has been also with the Institute of Problems of Informatics of Academy of Sciences of the Republic of Tatarstan, Russia as a leading researcher. He is now with the Research Program of Applied Mathematics and Computations of Mexican Petroleum Institute as a leading researcher and head of the project. He is also a visiting professor of the Center of Computing Research of the National Polytechnic Institute of Mexico. He is a Past President of Russian Association for Fuzzy Systems and Soft Computing, member of the International Fuzzy Systems Association (IFSA) Council, member of the Council of the Mexican Society of Artificial Intelligence, IEEE member, Honorary Professor of Budapest Tech and Honorary Researcher of the Republic of Tatarstan, Russia. He is a coauthor of 5 books. His areas of research activity are Fuzzy Logic and Fuzzy Reasoning, Hardware Implementation of Fuzzy Systems, Computational Intelligence, Expert Systems, Decision Making, Cluster Analysis, Time Series Data Mining and Data Visualization.

APPL. COMPUT. MATH., V.XX, N.XX, 20XX Didier Dubois is a Research Advisor at IRIT, the Computer Science Department of Paul Sabatier University in Toulouse, France and belongs to the French National Centre for Scientific Resarch (CNRS). He holds a Doctorate in Engineering from ENSAE, Toulouse (1977), a Doctorat d'Etat from Grenoble University (1983) and an Honorary Doctorate from the Facult Polytechnique de Mons, Belgium (1997). He is the coauthor, with Henri Prade, of two books on fuzzy sets and possibility theory, and 15 edited volumes on uncertain reasoning and fuzzy sets. Also with Henri Prade, he coordinated the HANDBOOK of FUZZY SETS series published by Kluwer (7 volumes, 1998-2000) including the book Fundamentals of Fuzzy Sets. He has contributed about 200 technical journal papers on uncertainty theories and applications. He is Editor-in -Chief of the journal Fuzzy Sets and Systems, an Advisory Editor of the IEEE Transactions on Fuzzy Systems, and a member of the Editorial Board of several technical journals, such as the International Journals on Approximate Reasoning, General Systems, Applied Logic, and Information Sciences among others. He is a former president of the International Fuzzy Systems Association (1995-1997). He received the 2002 Pioneer Award of the IEEE Neural Network Society, and the 2005 IEEE TFS Outstanding Paper Award. His topics of interest range from Artificial Intelligence to Operations Research and Decision Sciences, with emphasis on the modelling, representation and processing of imprecise and uncertain information in reasoning and problem-solving tasks. Janette Cardoso -graduated from the Universidade Federal de Santa Catarina, Brazil. She received her PhD degree in Computer Science from the University of Toulouse, France, in 1990 and her Habilitation degree in 2007. Presently, she is Professor at ISAE -Institut Supérieur de l'Aéronautique et de l'Espace, Toulouse, France. Her research interests are focused on Possibility theory, Petri nets and formal methods for embedded systems.