Model-Based Fault Diagnosis System Verification Using Reachability Analysis

In model-based fault detection and isolation (FDI) systems, fault indicating signals (FISs) such as residuals and fault estimates are corrupted by various noises, uncertainties and variations. It becomes challenging to verify whether an FDI system still works or not in real life applications. It is also challenging to select a threshold so that false alarm rate and missed detection rate are kept low depending on real operation conditions. This paper proposes solutions to the aforementioned problems by quantitatively analyzing the effect of uncertainties on FIS. The problems are formulated into reachability analysis problem for uncertain systems. The reachable sets of FIS are calculated under normal and selected faulty cases, respectively. From these reachable sets, the effectiveness of an FDI system can be qualitatively verified under described uncertainties. A dedicated threshold can be further chosen to be robust to all possible described uncertainties. As a by-product, the minimum detectable fault can also be quantitatively determined by checking the intersection of the computed reachable sets. The proposed approach is demonstrated by evaluating an FDI algorithm of a motor in the presence of parameter uncertainties, unknown load, and sensor noises, where a fault estimation-based approach is adopted to diagnose amplifier, velocity, and current sensor faults.

approach [6], [7], a fault is indicated when fault estimates deviate from a predefined threshold.
The most important and challenging issue in model-based FDI is robustness [4].Since there always exist some mismatches between the real plant and the mathematical model used for FDI observer design in real life applications, such as system parameter uncertainties (e.g., resistance in a motor), external disturbances (e.g., unknown load) and sensor noises [8], [9].As a result, it brings many challenges to the problem of FDI algorithm selection and verification, i.e., proving whether an FDI algorithm works or not in the presence of all kinds of system uncertainties and variations.Any product or function must go through verification and validation processes before being applied in real engineering to see whether or not certain performance specifications are satisfied [10].Besides, the FIS is also inevitably corrupted by all kinds of uncertainties, i.e., being nonzero even without any faults [11], [12], which brings challenges to the problem of threshold selection.If the threshold is selected too small, false alarm rate is high while missed detection rate is high if the threshold is selected too large.To this end, a lot of algorithms have been proposed to handle the robustness issue.According to in what stage uncertainties are considered, they can be categorized into active and passive approaches [12], [13].The active ones achieve robustness in the stage of observer design, while the passive ones achieve robustness in the stage of threshold selection (or residual evaluation, decision making).
Typical active approaches reduce or even decouple the effect of uncertainties on FIS by designing dedicated observers, such as unknown input observer [8], robust filter (e.g., H ∞ observer) [14].There is another important active approach, i.e., interval/set-membership-based diagnosis approach [12], [13], [15]- [17].In this approach, the conventional FDI observer is replaced by new interval/setmembership-based observers which can capture the possible upper and lower bounds of system states in the presence of interval uncertainties.Then a fault is alarmed by consistency check between the calculated state interval based on uncertain model and measurements from sensors.This approach has recently been applied to FDI problem of wind turbine system in [12] and [17], where uncertain set is represented by zonotopes due to its computational efficacy.There is also some interesting research on how to represent uncertainty using constrained zonotopes [16] such that a better tradeoff between computation accuracy and efficiency can be made.
This work is licensed under a Creative Commons Attribution 3.0 License.For more information, see http://creativecommons.org/licenses/by/3.0/Although FDI observer design, algorithm verification, and threshold selection are equally important for a successful fault diagnosis system in real applications, less attention has been paid to algorithm verification and threshold selection in the existing literature.The commonly used approach in the community for algorithm verification is stochastic simulation [18].This approach obtains the envelope of states or outputs based on a finite (possibly large) number of different linear models selected from a continuum of models corresponding to each possible value of the system parameters.However, as pointed out in [19] and [20], the number of required simulations grows exponentially with the number of state, input, and parametric variables due to a necessary gridding of the multidimensional set bounding all variables.Besides, even a large number of different models are evaluated on fairly fine grids of uncertain parameters or Monte Carlo sampling, it is still possible to miss the model corresponding to the most critical parameter combinations [10].As a result, more efficient and reliable verification approaches are needed for FDI systems, which motivates the research in this paper.
On the other hand, a threshold of the FIS should be selected such that FDI system is robust to all possible model uncertainties, unknown inputs, and faults of no interest (see [21, p. 8]).The existing robust threshold selection methods are mainly based on uncertainty amplification/inflation [21] (see [22, p. 289].In this approach, a conventional FDI observer is first designed, and then in the stage of FIS evaluation the system uncertainties are modeled (or approximated) as unknown input vector with limited bound or bounded by a known function and their effect on FIS is amplified based on norm inequality [21], integral inequality [22].This method may result in conservative or even useless threshold, especially in the presence of multiple parameter uncertainties (see [21, p. 251]).Consequently, more reliable threshold selection approaches are needed, which will also be investigated in this paper.
We aim to develop an efficient and reliable verification technique for FDI algorithm using reachability analysis.Particularly, the reachability analysis tool for continuoustime systems using zonotopes to represent uncertainties in [19] and [20] is employed due to its computational efficiency and reliability in comparison with optimization-based interval analysis [15].This technique has recently received much attention in control system verification [20].But to the best of our knowledge, we are the first to introduce this technique to fault diagnosis system verification.The basic philosophy of the proposed solution is to calculate the effect of all types of uncertainties on FIS under normal and selected faulty conditions, and conduct system verification and threshold selection by inspecting the differences between the calculated reachable sets.Specifically, interval models are first employed to describe the physical plant, which can accommodate the effect of uncertainties and noises.Suppose a candidate observer-based FDI system has been predesigned for the uncertain model along with the generated FIS.The effect of system uncertainties on FIS under normal case and selected faulty cases are evaluated using reachability analysis, resulting in normal FIS reachable set and faulty FIS reachable sets, respectively.It should be highlighted that different from the work in [12], [16], and [17] where interval/set theory is used to develop new FDI algorithms under uncertainties, the reachable set theory in this paper is used for a different purpose, i.e., develop a new verification technique for the existing observer-based FDI algorithms.Furthermore, in interval/setmembership-based diagnosis approaches, they do not actually evaluate the FIS of an existing FDI observer since the system state intervals rather than FDI observer state intervals are calculated.Based on the aforementioned reachable sets, several objectives can be achieved (see Section III for details).
1) It can be qualitatively verified whether a candidate FDI algorithm works or not in the presence of described uncertainties.2) An appropriate dedicated threshold can be quantitatively selected based on the normal reachable set SO that it is robust to described uncertainties.
3) The minimum detectable fault by the candidate FDI can be determined by checking the intersection of normal and faulty reachable sets with different fault amplitudes.The proposed verification approach is demonstrated by evaluating an FDI system of a motor in the presence of parameter uncertainties, unknown load, and sensor noises, where a fault estimation-based approach is adopted to diagnose amplifier, velocity, and current sensor faults.The results show that the proposed solution cannot only verify a candidate FDI algorithm, but also help optimally choose a right threshold and determine the minimum detectable fault under the prescribed level of uncertainties.

II. PROBLEM FORMULATION AND REACHABILITY ANALYSIS
In this section, the problem of FDI will be first formulated including model-based fault diagnosis, fault estimation-based approach and the challenges therein.Then reachability analysis for uncertain systems will also be introduced, which plays a key role in FDI verification and threshold selection.
A. FDI Problem 1) Model-Based FDI System: Consider an uncertain continuous-time linear system in the presence of unknown disturbance, actuator/sensor faults, and sensor noises, given by where x, u, f a , y d , and f d are the system state, control input, actuator fault, desired output, and unknown disturbance, respectively, and A, B, D, and R are their distribution matrices.The term Dy d is introduced to facilitate integral control design in state space models to remove steady-state control error, and so system (1) does not represent the physical plant but is a modified system for controller and FDI observer design.y m , f s , and n are the measurements, sensor faults, and sensor noises, respectively, and C m , S, and N are their distribution matrices.A is the system matrix perturbation to account for the effect of parameter uncertainties and each element is either zero or a bounded interval.It is supposed that feedback control law u = −K c y m + K d y d is predesigned since this paper mainly focuses on FDI rather than robust control.
The objective of model-based FDI is to detect the presence of fault f a and f s (i.e., fault detection) and determine what kind of fault has occurred when faults are detected (i.e., fault isolation) based on system model (1).There are a vast number of model-based FDI methods, some simply designed directly based on the normal model by ignoring uncertainties and noises while others considering uncertainties in their algorithm development.The focus of this paper is about the evaluation of a model-based FDI method, not about the development of a new FDI algorithm.Moreover, the system verification and threshold selection approach to be developed is applicable to both residual and fault estimation-based approaches.Without loss of any generality the fault estimation-based algorithm [6] is selected as a candidate FDI algorithm, which will be introduced in the following section.
2) Fault Estimation-Based Diagnosis System: Fault estimation-based approaches perform fault diagnosis by directly estimating faults based on observer theory, such as Kalman filter type observers [23], disturbance observers [3], [6] among others.The diagnosis logic is that when faults are estimated, one can directly tell whether a fault has occurred or not and which fault has occurred by checking whether the fault estimate deviates from a predefined interval centered at zero.
To obtain fault estimates for system (1), the concept of state augmentation is applied, which obtains fault estimates by augmenting faults as additional states [24].Since in real applications the external disturbances (e.g., unknown load in a motor) f d has an enormous effect on state estimation and consequently FIS, we reduce its effect through estimating it by augmenting it as an additional state as well.More advanced algorithms such as unknown input decoupling observer [23] for state or fault estimation can be considered in the future.
Letting x = (x, f a , f s , f d ) be the augmented state, system (1) can be equivalently represented by where , and C = C 0 S 0 are the augmented system matrices.A state observer can be designed for uncertain system (2) to obtain the estimate of augmented state x including the FIS (i.e., the fault estimates) Remark 1: There are some assumptions on Ā in (2) so that the observer action can be preserved.The uncertainties should not affect the observer stability.Many results are available to handle robust observer design such as zero sensitivity observer, robust observer with D-stability and H ∞ observer [25].In this paper, pole placement is adopted based on nominal model of (2).The performance of FDI observer (3) under prescribed level of uncertainties including stability can be automatically verified by checking the reachable set of state estimation errors using the proposed reachability analysis-based verification framework in Section III [see Fig. 2(bottom)].
Remark 2: The gain matrix Ko will usually have influence on state estimation and fault diagnosis performance such as fault detection time, robustness against uncertainties, sensitivity to fault and minimum detectable fault.Interested readers may refer to [9] for further information, where the effect of observer gain matrix on interval observer-based fault detection has been rigorously investigated.
Combing ( 2) and ( 3), the error dynamic e = x − x can thus be obtained When extended state x is obtained, one can obtain the actuator fault f a and sensor fault f s , which serve as the FIS.After FISs are available, one should further choose an appropriate threshold evaluating FIS such that a Boolean decision can be made-normal or faulty.The overall diagram of fault estimation-based diagnosis approach is shown in Fig. 1, 3) Challenging Issues: Given a candidate FDI observer (3), one can see from (4) that state estimation error and consequently fault estimates are inevitably corrupted by uncertainties Āx and noises.It should be noticed that this phenomenon is inevitable in all kinds of model-based fault diagnosis algorithms such as residual-based [14] and parameter estimation-based [4].This is also the reason why robustness is usually seen as the most important and challenging issue of model-based FDI algorithms [4].
This phenomenon results in several challenges to the application of FDI algorithm in real engineering: first, how to verify whether an FDI algorithm with a given threshold is still valid or not in the presence of all kinds of uncertainties, i.e., qualitative algorithm verification; second, how to choose an appropriate threshold such that the FDI is robust against uncertainties and noises, i.e., quantitative robust threshold selection; third, how to quantitatively determine the minimum detectable fault by a given FDI algorithm.This paper proposes solutions to the aforementioned challenges.The following reachability analysis is introduced, as it plays a key role in FDI algorithm verification and robust threshold selection in this paper.IEEE TRANSACTIONS ON SYSTEMS, MAN, AND CYBERNETICS: SYSTEMS

B. Reachability Analysis
The tool for reachability analysis of continuous-time uncertain systems is based on known technique in [19] and [20], which has been applied to control algorithm verification [20].Reachability analysis is typically performed iteratively for short time intervals τ k := [t k , t k+1 ] with t k := ki, where k ∈ N and i ∈ R + denote time step and step-size.The time step-size i is a tuning parameter making a tradeoff between computation efficiency and computation accuracy for reachability analysis.The iterative computation of reachable sets requires set-based addition and multiplication defined as follows.
Definition 1 (Set-Based Addition/Multiplication): The rule for set-based addition and multiplication are defined as Besides, there are multiple ways to represent a set including polytopes [26], zonotopes [20], ellipsoids [27], and support function [28].The zonotopes are preferred to represent a reachable set in this paper since they can efficiently represent reachable sets in high-dimensional spaces while operations required for reachability analysis can efficiently be applied.The definition of a zonotope is given as follows.
Definition 2 (Zonotope): Given a center c ∈ R n and so-called generators g (i) ∈ R n , a zonotope is defined as which can be put into a compact form as Z = (c, g (1) , . . ., g (p) ).The order of a zonotope is defined as ρ := (p/n), where p is the number of generators.The zonotope can also be seen as the set-based addition of line segments [−1, 1]g (i) .
Besides, other functions (such as the convex hull, Cartesian product) of two zonotopes are referred to [20].
1) Computation Consider an uncertain continuoustime linear system described by a differential inclusion where the uncertain system matrix A ⊂ I n×n , initial state x(0) ∈ X 0 ⊂ R n , u c ∈ R n is the known input, and u(t) ∈ V ⊂ R n is the uncertain input.The reachability can over-approximately (but in a tight way) obtain the reachable set R([0, r]) is a solution of (5), t = r, x(0) ∈ X 0 }, which can be seen as the exact reachable set.It should be noticed that the computation of exact reachable set is an open problem and consequently over-approximation is usually preferred.
The detailed algorithm and its implementation are referred to [20].However, for the completeness of the paper, its basic algorithm structure for a time interval τ k is given as follows.Suppose that the reachable set of the affine dynamics , the reachable set of the particular solution due to the uncertain input u(t) is R d p (u(t), t), and the partial reachable set correcting the initial assumption that trajectories are straight lines between t k and t k+1 is R d .According to [19] and [20], the reachable set for a time interval τ k is computed in the following steps.5) is given by A = A 0 + A, where the normal matrix and uncertain matrix are given by The initial state is chosen as x 0 ∈ 1 3×1 + [−0.1, 0.1].Each element of the known control input u c is a step input with amplitude 1 at 3 s (termed control 1) where its initial value is chosen as 0 (termed control 0) and the amplitude is increased to 2 at 6 s (termed control 2).Each element of the uncertain input u(t) lies in the interval [−0.05, 0.05].The time stepsize i is chosen as 0.005 s and the order of zonotope ρ is 600.Under this setting, the computation time for reachability analysis is about 2.5 s using a desktop with Intel Core i5-3570 CPU@3.4GHz and RAM 8 GB.The state reachable sets of x 1 and x 3 (gray area) and 500 stochastic simulation trajectories with computation time 35 s by Monte Carlo stochastic sampling (black lines) are shown in Fig. 2.
One can see from Fig. 2 that: 1) all the stochastically simulated trajectories lie in the calculated reachable set, which verifies the effectiveness of reachability analysis tool and 2) when the input amplitude is increased, the corresponding width of reachable sets also gets larger.This is due to the fact that when input amplitude is increased, the state amplitude and consequently the effect of uncertainties Ax on system states will get large as well.

III. ALGORITHM VERIFICATION
AND ROBUST THRESHOLD In this section, the main results of the paper are provided including qualitative verification of a candidate FDI algorithm with a given threshold and quantitative robust threshold selection for a given FDI system.As a by-product, the minimum detectable fault by a given FDI algorithm can also be determined.These results are discussed in detail in the following sections.

A. Fault Diagnosis Algorithm Verification
Algorithm verification answers the question that whether a given algorithm is still valid or not in a realistic environment [10] (i.e., in the presence of all kinds of system uncertainties such as parameter uncertainties, external disturbances, and sensor noises).It can be seen as a bridge between academia and industry since any algorithm should go through verification and validation to see whether certain performance specification are satisfied before being applied in industry.Although in the past few decades there are a lot of modelbased FDI algorithms proposed in academia (see [3], [4]), little attention has been paid to FDI algorithm verification and, as a result, model-based FDI algorithms have not been extensively applied in industry.The commonly used algorithm verification tool in engineering is stochastic simulation, i.e., simulating the real system by choosing a large number of system models with different parameter combinations.However, as pointed out in [20], the number of required simulations grows exponentially with the number of state, input, and parametric variables and more importantly, it is still possible to miss the most critical parameter combinations.Consequently, efforts should be made to bridge this gap, i.e., proposing a systematically formal approach for FDI algorithm verification.The research in this paper is conducted by following this observation.
The basic steps of the proposed qualitative verification approach for a candidate FDI algorithm is summarized as follows (see also Fig. 3).
1) Different intervals are used to capture system uncertainties (e.g., parameters, unknown inputs, and sensor noises), and an uncertain interval model [e.g., (1)] is obtained to represent the physical system; a candidate FDI system is also selected and developed [e.g., (3)].
2) The FIS reachable set is calculated in the absence of the faults with system uncertainties described in step 1.
3) The FIS reachable set is calculated under selected faulty cases, i.e., with system uncertainties and selected faults of interest.4) Qualitatively determine whether a candidate FDI is still valid or not under described level of uncertainties by comparing the normal and faulty FIS reachable sets.The rules are summarized in Proposition 1. Proposition 1: If the faulty FIS reachable set is well separated from the normal FIS reachable set under described level of uncertainties, then one can qualitatively conclude that the FDI is valid under the prescribed uncertainties.
The calculation of FIS reachable sets under normal and faulty cases are discussed in Section III-C, while the quantitative robust threshold selection and minimum detectable fault determination are discussed in Section III-B.

B. Robust Threshold and Minimum Detectable Fault
The threshold is selected to evaluate FIS and consequently a Boolean decision-normal or faulty is made.Without taking into account uncertainties (i.e., ideal condition), the threshold is usually set to be a small value close to zero since state estimation error and consequently FIS approaches to zero in the steady-state under normal case.However, we can see from (4) that state estimate error and consequently fault estimates are subject to the effect of uncertainties Āx and sensor noise Nn, which means the fault estimates are not zero even under normal conditions.So a threshold should be carefully selected such that it is robust against these uncertainties.
The commonly used robust threshold selection approach is approximating all the uncertainties using unknown input with limited bounds or a bounding function [21].There are two problems within this approach.First, it is not an easy task to find an appropriate unknown input to cover the uncertainties due to the time-varying and unknown nature of state x.Second, it is pointed out in [21, p. 251] that this approach could lead to a conservative or even useless threshold, since valuable information about the structure of model uncertainties has not been taken into account.To this end, we choose the robust threshold by quantitatively analyzing FIS reachable sets under normal and selected faulty cases as discussed in Section III-A.The results on robust threshold selection is summarized in Proposition 2.
Proposition 2: A robust threshold can be chosen based on the upper and lower bounds of the calculated normal FIS reachable set.If the fault estimate lies in the normal reachable set, it is concluded that no fault appears, and if the fault estimate deviates from the normal reachable set, it can be concluded that a fault has occurred.
In addition, by comparing the normal reachable set and faulty reachable sets under different fault amplitudes, the minimum detectable fault by a candidate FDI system can also be determined.If the intersection of the FIS reachable sets under normal and faulty cases is empty, it can be concluded that the FDI algorithm can effectively detect the presence of fault, while if the intersection is not empty (e.g., the amplitude of fault is too small), then the FDI algorithm fails to detect this particular type of fault.The results on minimum detectable are summarized in Proposition 3. Proposition 3: Denote f m the minimum fault amplitude that can be detected by a candidate FDI system, it can be determined in the following way: the intersection of the FIS reachable sets under normal and faulty cases is not empty for the scenario f ≤ f m , and the intersection of these two sets is empty for the case f ≥ f m .
The overall diagram of the proposed FDI verification procedure is shown in Fig. 3.

C. FIS Reachable Set Computation
In this section, FIS reachable set computation under normal and selected faulty cases are discussed.The aforementioned FIS reachable set computation problem is first transformed into the problem of state reachable set computation so that the existing tool on reachability analysis introduced in Section II-B can be applied.
To facilitate analyzing the effect of uncertainties Āx and sensor noises Nn on the augmented state x and consequently fault estimates, we derive the dynamic of x and put it into the format of ( 5) such that the technique of state reachable set computation can be adopted.On the one hand, from e = x − x, the state estimation error dynamic (4) can be put into the following form: On the other hand, we can see from (3) that the time varying input u = −K c y m + K d y d and y m are involved in derivation of x, which makes the reachable set computation complicated.To simplify the problem, we substitute the control input u and y m into (3), such that the x dynamics are given by ẋ Combing ( 6) and ( 7), we can obtain the following composite dynamics including extended state estimation error e and extended state estimate x: Using the under-brace notations in ( 8), the composite dynamics ( 8) can be put into a compact form, given by where χ = [e T , xT ] T , A is the uncertain system matrix, u c is the known input, and u(t) is the uncertain input.One can see that ( 9) falls into the same format as that of the uncertain system ( 5) in Section II-B, which means the existing reachable set computation tool in Section II-B can be applied.Then we can obtain the reachable set of χ and consequently FIS reachable sets under normal case and selected faulty cases, since the FIS are chosen as the fault estimates in fault estimation-based approach which are the elements of χ .Remark 3: We can see from ( 8) that the distribution matrices of known input y d and uncertain input n are different.To make the reachable set computation problem of system ( 8) solvable based on the existing tool [i.e., putting (8) into ( 9)], the original uncertain input n has been transformed into u(t).This process will result in conservativeness in reachable set computation due to the fact that the dimension of uncertain input has been increased and consequently the single-useexpressions [19] of interval computation cannot be achieved.However, future work can be done to handle this problem and consequently remove the conservativeness.

IV. CASE STUDY: MOTOR FAULT DIAGNOSIS
In this section, a motor system is chosen as case study since it is simple but enough to illustrate the basic idea of the algorithm proposed in this paper.In fact, it is quite challenging to verify the FDI algorithm for this system (especially in the presence of multiple faults simultaneously) as shown later, the FIS is not only affected by system parameter variation and sensor noises but also by unknown external load.Both actuator fault (i.e., amplifier) and sensor fault (i.e., velocity and current sensor faults) diagnosis are considered.The overall diagram of motor FDI system is shown in Fig. 4, where the closed-loop control system using state space feedback with integral action is introduced in Section IV-A, the actuator fault estimator, sensor fault estimator, and their corresponding robust threshold selection are given in Sections IV-B and IV-C, respectively.

A. Motor System
A linear model for the motor system with actuator and sensor faults can be represented by (1) ω, I, (θ − θ d )dt] T are the system states representing position, velocity, current, and position error integral [note: x 4 = (θ − θ d )dt is introduced so that an integral control can be designed to remove the steady-state tracking error], u is the control voltage, which is designed as where θ m , ω m , and i m are the measurement value of position, velocity, and current, respectively.The parameters The meaning of the aforementioned parameters and their normal values are referred to the Appendix.We suppose there are four key parameters with uncertainties, i.e., k d , k T , k e , and resistance R. The uncertain levels are given as follows: where k dn , k Tn , k en , and R n denote the normal values.These parameter uncertainties can be grouped into interval matrix A. Without loss of generality, the sensor noises n p , n v , and n i lie in a bounded interval ∈ [−0.001, 0.001] with uniform distribution.

B. Amplifier Fault Diagnosis
In this section, fault estimation-based approach is applied to the problem of amplifier fault diagnosis.In this scenario, we suppose no sensor fault occurs to satisfy the observer observability and only actuator fault and unknown load are treated  as the additional states in generalized state observer design, i.e., x1 = (x, f a , f d ).The FDI observer is designed as where the system matrices are given by , and C1 = C 0 0 .
The observer gain Ko1 can be found in the Appendix.The reference position is chosen as a step signal at 0 s with amplitude of 1 rad, while the unknown external load and amplifier fault (note: additive step-type fault is considered since other types of fault can be equivalently transformed into additive one) with different amplitudes are plotted in Fig. 5.
The initial state of both control system and observer system are supposed to be zero vector.The reachable set computation results for both normal and amplifier faulty cases are shown in Figs. 6 and 7 where the gray areas are the calculated reachable set and the black lines are 500 stochastic simulations of the amplifier fault estimate using Monte Carlo stochastic sampling for each uncertain variable satisfying a uniform distribution in the predefined bounded interval.
We can see from Figs. 6 and 7 that reachable set computation can effectively capture the fault estimates, since all simulated fault estimates lie in the computed reachable set for both normal and faulty cases.When the unknown load is changed from 1 to 0.5 Nm at 3 s, the amplitude of reachable set reduces accordingly.This is because when load amplitude is reduced the state amplitude and so the effect of uncertainties Ax on fault estimate will be reduced.First, based on the two reachable sets in Fig. 6, we can qualitatively verify that the fault estimation-based diagnosis algorithm is valid in the presence of described system uncertainties and sensor noises, since the FIS reachable set under the normal case is close to zero and the FIS reachable set under the faulty case substantially deviates from the normal FIS reachable set.Besides, based on the principle of the proposed threshold selection approach, the threshold can be chosen as the upper and lower bounds of the normal reachable set, i.e., [0.1, 0.1].
Second, we can see from Fig. 6 that in the case of fault amplitude f a = 1 V, the FIS reachable set (after 6 s) does not intersect with the normal FIS reachable set (after 6 s) and consequently the fault diagnosis algorithm can effectively detect the presence of fault.However, in the case of fault with a smaller amplitude f a = 0.1 V, as shown in Fig. 7 the fault reachable set intersects with the normal reachable set and so the fault diagnosis algorithm fails to detect the presence of the actuator fault with this particular amplitude.This phenomenon is reasonable since when the fault amplitude is too small, its effect on the fault estimate and consequently FIS is limited and cannot be distinguished from the effect of system uncertainties.
To determine the minimal fault amplitude that can be detected by the given fault diagnosis algorithm under described uncertainties, simulations are further done by iteratively reducing the fault amplitude such that the lower bound of FIS reachable set under faulty case (in steady state) coincides with the upper bound of FIS reachable set under normal case (fault with positive value is assumed, similar process can be applied to the case of fault with negative value).This process has been done and the final results are shown in Fig. 8, where similar to previous cases the amplifier fault is supposed to occur at 6 s.
From Fig. 8, we can first determine the lower and upper limits of the FIS under normal case, i.e., [0.1, 0.1], which serves as the threshold interval.The minimal fault amplitude that can be detected by the given fault diagnosis algorithm under described uncertainties is 0.16 V, which means the amplifier faults with amplitude lower than 0.16 V are not detectable (or false alarm will appear) due to the effect of all kinds uncertainties and unknown external load.
Remark 4: The reachable set computation is carried out in MATLAB 2012 using the algorithm in Section II-B in conjunction with the existing toolboxes including "INTerval LABoratory" and Multi-Parametric Toolbox.The time stepsize of reachable set computation is chosen as i = 0.0015 and the order of zonotope ρ is 600.Under this parameter setting, the computation time is about 60 s using the same software as in Section II-B; the increase in computation time is mainly due to the increased vector dimension from three to twelve.To further improve the computation precision, one can reduce the step-size or increase the order of zonotope, however, it will result in higher computation load.

C. Sensor Fault Diagnosis
In this section, we consider the problem of multiple sensor fault diagnosis, where velocity and current sensors are considered simultaneously.The fault diagnosis algorithm verification and threshold selection for multiple sensors are complex and challenging, as will be shown by the simulation results that the FISs are affected by many factors such that reference signal, uncertainties and external unknown load.Moreover, the velocity sensor fault has effect on current sensor fault estimate and vice versa, which makes the problem even more challenging.To satisfy the observability, it is supposed that no actuator fault occurs and only velocity sensor fault, current sensor fault and unknown load are treated as additional states, i.e., x2 = (x, f s , f d ).The fault diagnosis observer is given by ẋ2 where the system matrices are given by , and C2 = C S 0 and the observer gain matrix Ko2 is referred to the Appendix.The unknown external load, velocity and current sensor fault profiles are shown in Fig. 9.Other simulation scenario is the same as that of amplifier fault diagnosis in Section IV-B.The reachable set computation results for velocity and current sensor fault estimates are plotted in Figs. 10 and 11, where the gray areas are the calculated reachable set and black lines are the 500 stochastic simulations of the sensor fault estimates using Monte Carlo stochastic sampling for uncertain variables, which are supposed to satisfy a uniform distribution in predefined bounded intervals.The computation time increased to about 100 s, since the vector dimension has increased to fourteen.
Similar to the case of amplifier fault diagnosis, we can see from Figs. 10 and 11 that all stochastic simulations lie in the calculated FIS reachable sets, which verifies the effectiveness of reachability analysis.The initial reachable sets before 1 s is not centered at zero, this is due to the transient effect of step reference position.At 3 s, the load amplitude is reduced from 1 to 0.5 Nm and so the effect of uncertainties on FIS is reduced, i.e., the width of reachable set becomes smaller.From Fig. 10, we can obtain that in the presence of velocity sensor fault at 5 s, the FIS substantially deviates from the FIS under normal case and so velocity sensor fault can be alarmed.At 7 s, there is a transient jump in the reachable set of velocity fault estimate due to the presence of current sensor fault.Because the fault estimates for velocity and current sensor faults are coupled.We can see from Fig. 11 that at 5 s the velocity sensor fault also has effect on current sensor fault reachable set.In the presence of current sensor fault at 7 s, there is a substantial jump in FIS reachable set, which again verifies the effect of fault diagnosis system for this type of fault under described uncertainties.
The simulation results in Figs. 10 and 11 quantitatively show the effect of reference position, external unknown load, system uncertainties, and sensor noises on FISs.As a result, from the calculated FIS reachable sets, we can choose an appropriate threshold interval (for example, [−0.08, 0.06] for velocity sensor fault and [−0.18, 0.16] for current sensor fault) such that good robustness can be achieved.

V. CONCLUSION
This paper considered the problem of fault diagnosis system verification and robust threshold selection for a typical modelbased fault diagnosis algorithm.Due to the presence of uncertainties, disturbances and noises, the FISs deviate from zero even under normal condition, which brings challenges to algorithm verification and threshold selection.Reachability analysis is drawn to calculate the reachable sets of FIS under normal and selected faulty cases.Based on the calculated reachable sets, a candidate fault diagnosis system can be qualitatively verified whether or not it still works under described uncertainties.A robust threshold can be quantitatively selected which is robust to uncertainties while minimizing the false alarm rate or missed detection rate.So a much better tradeoff under the described uncertainties can be achieved.In addition, by comparing these two reachable sets under different fault amplitudes, the minimum detectable fault by a given fault diagnosis algorithm can also be determined.
A case study of amplifier and sensor fault diagnosis of a motor is given to illustrate its principle.It has shown that the combination of reference set-points change, external load variations, parameter uncertainties and sensor noises makes the verification of a simple residual-based FDI algorithm quite challenging.It can be seen that the proposed method is not only able to verify a candidate FDI but also help to optimally choose a right threshold under the prescribed level of uncertainties.As a result, a better tradeoff between false alarm and missed detection can be achieved.
This paper mainly focused on proposing the idea of fault diagnosis algorithm verification and robust threshold selection using reachability analysis for uncertain systems.With the advent of more advanced reachable set computation tools and increasing computation power, it is expected that the proposed solution would find a wide range of applications.For example, it can be applied to optimize the observer gain matrix for a given model-based fault diagnosis algorithm.Furthermore, it can also be applied to performance verification for other model-based fault diagnosis approaches (e.g., robust residualbased and parameter estimation-based), such that the most suitable algorithm for a given problem can be determined.

APPENDIX
The meanings of motor parameters and corresponding normal values are summarized in Table I

Fig. 1 .
Fig. 1.Diagram of fault estimation-based diagnosis approach including three elements: closed-loop control system, FIS generator, and FIS evaluation.

Fig. 2 .
Fig. 2. State reachable sets (gray areas) of x 1 (left) and x 3 (right) under three different control amplitudes and 500 Monte Carlo stochastic simulated trajectories (black lines).

Fig. 3 .
Fig. 3. Proposed procedure for qualitative FDI system verification, robust threshold selection, and minimum detectable fault determination.
are the control parameters to be designed (see the Appendix) with K d := k θ .y d θ d is the desired output position, d is the load, and y m is the measurement output.f a and f s are the actuator and sensor fault, respectively, where f s = f sv , f sc T with f sv and f sc being the velocity and current sensor faults.System matrix A, control input matrix B, desired output matrix D, load matrix R, and measurement matrix C m are given as follows:

Fig. 5 .
Fig. 5. Profile of external load (real line) and amplifier fault at 6 s with different amplitudes dotted line f a = 0.1 V and dashed dotted line f a = 1 V.

Fig. 6 .
Fig. 6.Reachable set of amplifier fault estimate under normal case and fault amplitude f a = 1 (gray area); and 500 Monte Carlo stochastic simulations of fault estimate (black lines).

Fig. 7 .
Fig. 7. Reachable set of amplifier fault estimate under normal case and fault amplitude f a = 0.1 (gray area); 500 Monte Carlo stochastic simulations of fault estimate (black lines).

Fig. 8 .
Fig. 8. Reachable set of amplifier fault estimate with fault amplitude f a = 0.16 V and f a = −0.16V (gray areas); the upper and lower bounds of normal FIS reachable set (two dashed lines); and 500 Monte Carlo stochastic simulations of fault estimate (black lines).

Fig. 10 .
Fig. 10.Reachable set of velocity sensor fault estimate (gray area) and 500 stochastic simulations of fault estimate.

Fig. 11 .
Fig. 11.Reachable set of current sensor fault estimate (gray area) and 500 stochastic simulations of fault estimate.
1) Starting from R d (t k ), compute the reachable set R d a (t k+1 ).2) Obtain the convex hull of R d (t k ) and R d (t k+1 ) to approximate the reachable set for the time interval τ k .

TABLE I NORMAL
VALUES OF MOTOR PARAMETERS