Syrine Derbel
email: syrine.derbel@ensea.fr

Florentina Nicolau

Nabih Feki

Jean Pierre Barbot

Mohamed Slim Abbes

Mohamed Haddar

Faults diagnosis via a dynamical sparse recovery method and application to a gearbox system

Keywords: Diagnostic, fault detection, sparse recovery diagnosis, gearbox system, gear power transmission, mechanical faults, sensor faults

With the ever increase of the complexity and the importance of industrial systems, diagnosis techniques allowing to detect, locate and identify any abnormalities in the system as early as possible, have attracted a lot of attention over the past years. In this paper, we present a diagnosis method for nonlinear dynamical systems, called sparse recovery diagnosis (SRD), based on a dynamical algorithm that estimates a sparse fault vector from few system measurements. The term sparse means that many faults can be considered but only few of them can occur simultaneously. In order to illustrate the performances of this diagnosis method, we apply it to a gear power transmission. This dynamical system is among the most important mechanical components in industrial systems. The gear power transmission model considered in this paper is composed by a two-stage gear for which we take into account the torsional eect of the gears. Dierent sensor and mechanical faults perturbing its operating mode will be modeled and detected by the SRD method.

Introduction

The diagnostic of industrial systems is one of the most important tasks that ensures the continuity of the production, the safety and the economic aspects in the industries. For these reasons, many diagnosis methods enabling early fault detection have been developed in the literature (see for example, the following surveys [START_REF] Hwang | A survey of fault detection, isolation, and reconguration methods[END_REF], [START_REF] Dong | A survey on distributed ltering and fault detection for sensor networks[END_REF]). The presence of faults is very often the cause of system performances degradation, and even the damage and the collapse of the whole system. An unexpected fault can aect either the sensors leading to measurement errors (in that case, the fault is called a sensor fault) or the components of the system leading to a change of the system properties (in that case, the fault is called an operating fault).

The faults diagnosis consists in three important steps: fault detection, fault location and fault identication. The rst step enables to detect the presence of a malfunction in the system. Fault isolation consists on determining the location of faulty components, while fault identication allows to identify the type, shape and size of the appeared fault. In general, diagnosis methods are based on the concept of redundancy (see , for example, [START_REF] Dong | A survey on distributed ltering and fault detection for sensor networks[END_REF], [START_REF] Gao | A survey of fault diagnosis and fault-tolerant techniques-part i: Fault diagnosis with model-based and signal-based approaches[END_REF]) which can be of two types: hardware or analytical redundancy. The key idea of the hardware redundancy is to compare duplicated measurements using multiple hardware (for instance, sensors) for the same components and input signals. This comparison allows to take the best diagnostic decision by using dierent methods such as majority voting, limit checking, etc. The analytical redundancy uses specic information given by a mathematical model or historic system data, or signals information. It includes the model-based methods, the signals-based methods and the data-driven methods. Model-based methods such as parity space techniques, observers' methods and parameters estimation have been the object of many papers (see, e.g., the surveys, [START_REF] Frank | Fault diagnosis in dynamic systems using analytical and knowledgebased redundancy: A survey and some new results[END_REF], [START_REF] García | Deterministic nonlinear observer-based approaches to fault diagnosis: A survey[END_REF], [START_REF] Isermann | Model-based fault-detection and diagnosis status and applications[END_REF]). Signals-based methods are discussed and applied in many works (see, e.g., [START_REF] Kimmich | Fault detection for modern diesel engines using signal-and process model-based methods[END_REF], [START_REF] Barszcz | Application of spectral kurtosis for detection of a tooth crack in the planetary gear of a wind turbine[END_REF], [START_REF] Lei | A review on empirical mode decomposition in fault diagnosis of rotating machinery[END_REF] and [START_REF] Dong | A survey on distributed ltering and fault detection for sensor networks[END_REF]). They contain time domain, frequency domain and time-frequency domain methods. Finally, data-driven techniques, including neural networks and fuzzy logic methods, have been the aim of several research works (see, e.g., [START_REF] Dai | From model, signal to knowledge: a data-driven perspective of fault detection and diagnosis[END_REF], [START_REF] Dong | A survey on distributed ltering and fault detection for sensor networks[END_REF]).

In this paper, we present a diagnosis method (that we will call sparse recovery diagnosis, or simply SRD method) for nonlinear dynamical systems. The SRD method is inspired from sparse recovery methods ( [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF], [START_REF] Becker | NESTA: A fast and accurate rst-order method for sparse recovery In[END_REF]) which are well know in the eld of signal processing (where they are often applied for medical imaging and image processing, see, e.g., [START_REF] Chen | Artifact suppressed dictionary learning for low-dose ct image processing[END_REF], [START_REF] Dong | Hyperspectral image superresolution via non-negative structured sparse representation[END_REF]). From a fault detection point of view, other methods based on similar theories (like compressing sensing theory, sparse decomposition and impulse signal recovery) have been already developed in the literature (e.g., [START_REF] Chen | Compressed sensing based on dictionary learning for extracting impulse components[END_REF], [START_REF] Tang | Sparse classication of rotating machinery faults based on compressive sensing strategy[END_REF], [START_REF] Du | Compressed-Sensing-Based Periodic Impulsive Feature Detection for Wind Turbine Systems In[END_REF], [START_REF] Frat | Compressive sensing for detecting ships with second-order cyclostationary signatures[END_REF], [START_REF] Lin | Impulse detection using a shift-invariant dictionary and multiple compressions[END_REF]). They are based on the extraction of the fault features through a sparse representation in a specic base (like Fourier, wavelet or cepstral transforms, etc, for a spectral representation). In many works on fault diagnosis via sparse recovery techniques, fault characteristics are directly extracted from a small number of random projections without reconstructing the vibration signal completely and with prior fault characteristic frequency knowledge required. On the contrary, the SRD method presented in this paper is based on an analytical system modeling and allows to detect, locate and identify dierent types of faults aecting the dynamical system. The analytic model contains the healthy and the faulty operating of the system and is given by a nonlinear state-space representation. Whereas the previous works on fault diagnosis via sparse recovery techniques consider the static case only, the SRD method allows to take into account operating faults that depend on time and/or state. This is one of its most important characteristics.

The sparsity concept is introduced by the fact that the number of faults considered in the system modeling is much greater than the number of monitoring sensors as well as by the simultaneous appearance of the faults. This means that only a restricted number among the modeled faults can appear simultaneously. This number depends essentially on the number of monitoring sensors. Indeed, the SRD method takes into account a very large number of possible faults that may act on the system, but it can be applied under the assumption that few of them can occur simultaneously in the system. This can be seen as a limitation of the method but it can be overcome by increasing the number of monitoring sensors (thus by adding new sensors) allowing to obtain more system information and, therefore, to diagnose more defects.

Another important point is that the SRD is an on line diagnosis method (where by on line we mean that all algorithms are solved in a continuous process and with the same computation step).

Since the SRD method is based on a state-space representation, it belongs to the class of model-based diagnosis methods. In the literature, many works apply or develop modelbased techniques in order to diagnose dierent faults (see , for instance, [START_REF] Isermann | Trends in the application of model-based fault detection and diagnosis of technical processes[END_REF], [START_REF] Combastel | Model-based and wavelet approaches to induction motor on-line fault detection In[END_REF], [START_REF] Zheng | A rolling bearing fault diagnosis method based on multi-scale fuzzy entropy and variable predictive model-based class discrimination[END_REF]). Very often they are based on residuals generation for the state variables or/and the output variables. The created residuals represent, in general, the dierence between the estimated healthy variables and the real ones. These methods allow to diagnose faults by checking the consistency between this dierence and a threshold that is xed to distinguish the healthy from the faulty operating of the system or by applying some temporal and/or frequency analysis. However, in some complex systems, model-based methods require high amount of measurements and computational eorts due to the existence of a large number of variables. Moreover, each fault should have its own diagnosis algorithm to be analyzed, detected and isolated, in particular, when similarities between dierent residuals occur. This may render the system even more complex or may lead to dealing with big-data problems in the case of many various faults. From this point of view, a main advantage of the SRD method is that it allows to consider many dierent faults in the system modeling and, moreover, their diagnosis is carried out by an on line diagnosis algorithm with nite-time convergence. Furthermore, the SRD method allows to detect, locate and isolate faults using a single algorithm. This is not always possible for others diagnosis methods (because of a dicult localization step that requires some specic knowledge of the system characteristics, due to the lack of the algebraic relations between dierent system variables).

To sum up, with the help of the SRD method, a complete diagnostic (detection, location and identication) can be obtained. To illustrate its performances, we apply it to diagnose faults in a gearbox system.

The gear power transmission has a crucial role in mechanical transmission systems and its diagnosis depends on how faults aect the system. Some faults (like misalignment, eccentricity) may produce amplitude modulations. On the other hand, frequency modulations can be caused by other types of faults (such as gears cracking or aking). Therefore, sometimes it may be dicult to identify the nature of the defects.

The gearbox system considered in this paper is composed by a two-stage gear and is affected by dierent mechanical and sensor faults. With the help of the SRD method, detection, identication and localization of the system faults are carried out. The simulation results show the accurate reconstruction of the mechanical faults and their diagnostic is carried out by verifying the nature of the obtained defect signal and by using some frequential/temporal characteristics of the gear elements.

The paper is organized as follows: the next section is devoted to the general presentation of the SRD method. Then, in the third section, we present the healthy and faulty models of the gear power transmission system and explain how the SRD method can be applied for this particular dynamical systems. The fourth section shows simulations results implemented on Matlab/Simulink and discusses the good performances of the SRD method. Finally, some conclusions are presented. Nonlinear dynamics of the system g j (x) Nonlinear dynamics of the faults Consider the following nonlinear state-space representation describing the healthy operating mode of the system:

w j jth operating fault, 1 ≤ j ≤ p y i Measurement output, 1 ≤ i ≤ m h i (x) Function of the measurement output, 1 ≤ i ≤ m s i Fault of the ith sensor, 1 ≤ i ≤ m
ẋ = f (x) y i = h i (x), 1 ≤ i ≤ m. (1)
where x ∈ R n is the state vector, f (x) is the nonlinear dynamics, y ∈ R m is the m-dimensional measurement output (corresponding to the sensor outputs). The point denotes the derivative with respect to an independent variable which in general represents the time. The functions h i , for 1 ≤ i ≤ m, are supposed independent. The appearance of faults in physical systems is often an inevitable phenomenon that, in general, modies their dynamical behavior. They are modeled by the following extension of representation (1):

ẋ = f (x) + p j=1 g j (x)w j (t), y i = h i (x) + s i , 1 ≤ i ≤ m. (2)
where (w

1 (t), • • • , w p (t)) T and (s 1 , • • • , s m )
T are vectors of, respectively, the system operating faults and the sensor faults. All g j (x), 1 ≤ j ≤ p, are known vector elds related to the dynamics faults. We suppose that f , g j and h i are suciently smooth. All sensor faults s i are assumed piecewise-constant functions (explaining the lack of the t-argument in (2)) of the form:

s i = 0, if 0 ≤ t < t i , a i , if t ≥ t i , (3)
for 1 ≤ i ≤ m, with a i the fault magnitude and t i the appearance time of the fault. The total number of faults q = p + m is supposed much greater than the number of monitoring sensors (see Assumption 3.2). It is clear that if some faults are equal to zero (i.e., w j = 0, for some 1 ≤ j ≤ p, or s i = 0, for some 1 ≤ i ≤ m), then those faults do not aect the system and thus they do not disrupt the operating of (2). The principal notations used in this section are recalled in Table 1.

Example 3.1. In order to illustrate the dierent notations required for the introduction of the SRD method, we consider the following system:

       ẋ1 = x 2 + x 3 + 2x 3 w 1 + x 2 3 w 2 , ẋ2 = 2x 2 3 + x 3 x 1 + x 3 w 2 + x 3 w 3 + x 3 w 4 , ẋ3 = x 1 + x 4 + x 3 w 2 + 2w 4 , ẋ4 = x 2 + x 2 3 + x 3 w 2 , (4)
with the measurements:

   y 1 = x 1 + x 2 + s 1 , y 2 = x 2 + s 2 , y 3 = x 4 + s 3 , (5)
where x = (x 1 x 2 x 3 x 4 ) T ∈ R 4 is the state vector, (w 1 w 2 w 3 w 4 ) T ∈ R 4 is the vector of operating faults, (s 1 s 2 s 3 ) T ∈ R 3 is the vector of sensor faults and

g 1 =     2x 3 0 0 0     , g 2 =     x 2 3 x 3 x 3 x 3     , g 3 =     0 x 3 0 0     , g 4 =     0 x 3 2 0    
are the vector elds showing how the operating faults w j act on the system. Notice that they are state dependent.

In order to apply the SRD method, it is important to work in coordinates in which each output corresponds to one (and only one) state variable of the system. Hence, the following change of coordinates may be required: introduce x = ϕ(x) ∈ R n , where ϕ(x) is a local dieomorphism whose m rst components are given by ϕ i (x) = h i (x) and the last n -m are any original coordinates completing the functions h i , 1 ≤ i ≤ m, to a coordinate system. Notice that the local dieomorphism ϕ(x) always exists. In the new coordinates, system (2) can be written as:

ẋ = f (x) + p j=1 gj (x)w j (t), y i = xi + s i , 1 ≤ i ≤ m, (6)
with

f (x) = ∂ϕ ∂x (x) • f (x) = ∂ϕ ∂x (ϕ -1 (x)) • f (ϕ -1 (x)), and
gj (x) = ∂ϕ ∂x (x) • g j (x) = ∂ϕ ∂x (ϕ -1 ((x)) • g j (ϕ -1 (x)),
where ϕ -1 is the local inverse of ϕ(x). Observe that in the new x-coordinates, the measurements are simply the m rst states xi , 1 ≤ i ≤ m. Notice that for mechanical systems, this is often the case since usually we measure the states of the system (displacements, rotational speed, etc). In order to simplify the notation, from now on, we will drop the tildes.

We denote the total number of faults by q (that is, q = m + p) and dene the following global fault vector including the defects aecting both the dynamics and the sensors:

d(t) =          s 1 . . . s m w 1 (t) . . . w p (t)          ∈ R q . ( 7
)
Example 3.1 (continuation). For system (4) together with measurements (5), the dieomorphism ϕ(x) can be taken as:

       x1 = x 1 + x 2 , x2 = x 2 , x3 = x 4 , x4 = x 3 , (8)
with inverse

       x 1 = x1 -x2 , x 2 = x2 , x 3 = x4 , x 4 = x3 . (9)
In the new x-coordinates we obtain:

           ẋ1 = x2 + x4 + 2x 2 4 + (x 1 -x2 )x 4 + 2x 4 w 1 + (x 2 4 +x 4 )w 2 + x4 w 3 + x4 w 4 , ẋ2 = 2x 2 4 + x4 (x 1 -x2 ) + x4 w 2 + x4 w 3 + x4 w 4 , ẋ3 = x1 -x2 + x3 + x4 w 2 + 2w 4 , ẋ4 = x2 + x2 4 + x4 w 2 , (10)
with the measurements (notice that each output h i corresponds to a new state variable xi )

   y 1 = x1 + s 1 , y 2 = x2 + s 2 , y 3 = x3 + s 3 . (11)
The global fault vector of ( 10)-( 11) is dened by

d(t) = [s 1 s 2 s 3 w 1 w 2 w 3 w 4 ] T ∈ R 7 .
Denition 3.1 ((Relative degree with respect to the operating faults)). The relative degree of the output component h i (x) = x i , 1 ≤ i ≤ m, is the positive integer ρ i dened as follows:

L g j L f h i = 0, for all 0 ≤ < ρ i -2, 1 ≤ j ≤ p, L g j L ρ i -1 f h i = 0, for some 1 ≤ j ≤ p, (12)
where L f h i denotes the th-order Lie derivative of the function h i along the vector eld f (see Denition 3.2 below).

Denition 3.2 ((Lie derivative)). The Lie derivative of the function h i , 1 ≤ i ≤ m, along a vector eld f can be computed by:

L f h i (x) = ∂h i (x) ∂x • f (x).
Recursively, we can dene, for ≥ 1, the th-order Lie derivative by:

L f h i (x) = L f (L -1 f h i (x)) = ∂L -1 f h i (x) ∂x • f (x), where L 0 f h i (x) = h i (x).
Hence, the term L g j L f h i of ( 12) means

L g j L f h i (x) = ∂L f h i (x) ∂x • g j (x).
Based on Denition 3.1, the relative degree ρ i with respect to the operating faults can be calculated for each h i . The integer ρ i indicates the number of times that we have to dierentiate before the eect of the operating faults can be seen on the output components h i . In other words, ρ i is the order of the rst derivative of h i that explicitly involves some of the operating faults w j and provides information of w j that will be used to detect them. Indeed, recall that y i = h i (x) + s i = x i + s i , for 1 ≤ i ≤ m (see [START_REF] Besançon | Nonlinear observers and applications[END_REF]), where all s i are piecewise-constants thus, by Denition 3.1, for 1 ≤ i ≤ m, we have:

y (ρ i ) i = L ρ i f h i (x) + p j=1 (L g j L ρ i -1 f h i (x))w j (t). (13)
In order to extract the sensor faults, we replace the states x i , for 1 ≤ i ≤ m, in (13) by their measured values y i = x i + s i (the real measurements aected by the faults)

y (ρ i ) i =L ρ i f h i (y, x) + p j=1 (L g j L ρ i -1 f h i (y, x))w j (t), (14)
where x denotes the unmeasured states x = (x m+1 , . . . , x n ).

Notice that expression (14) now involves explicitly s i (via y i ) and this dependence may be nonlinear. In this paper, we treat only the case when the right-hand side of ( 14) is ane with respect to s i , 1 ≤ i ≤ m, (i.e., of form [START_REF] Derbel | Application of homogeneous observers with variable exponent to a mechatronic system[END_REF]). This assumption may seen restrictive, but if ( 14) is nonlinear with respect to s i , then a Taylor approximation of the rst order can be used (see [START_REF] Torki | A novel FDI sparse recovery method: Application on PMSG wind turbine[END_REF]). With this assumption, ( 14) can be written as follows:

y (ρ i ) i =L ρ i f h i (x) + m j=1 σ ij (x)s j + p j=1 (L g j L ρ i -1 f h i (x))w j (t). (15)
Example 3.1 (continuation). For system [START_REF] Chen | Dynamic simulation of spur gear with tooth root crack propagating along tooth width and crack depth[END_REF], we compute the Lie derivative of each h i , for 1 ≤ i ≤ 3, along the vector elds g j , for 1 ≤ j ≤ 4, and obtain ρ i = 1, for 1 ≤ i ≤ 3.

In order to extract the sensor faults, we replace the states x1 , x2 and x3 in (10) by their measured values given by [START_REF] Chen | Compressed sensing based on dictionary learning for extracting impulse components[END_REF]. We obtain the following equations represented in form [START_REF] Derbel | Application of homogeneous observers with variable exponent to a mechatronic system[END_REF]:

       ẏ1 = L f h 1 (x) + x4 s 1 + (1 -x4 )s 2 + 2x 4 w 1 + (x 2 4 +x 4 )w 2 + x4 w 3 + x4 w 4 , ẏ2 = L f h 2 (x) + x4 s 1 -x4 s 2 + x4 w 2 + x4 w 3 + x4 w 4 , ẏ3 = L f h 3 (x) + s 1 -s 2 + s 3 + x4 w 2 + 2w 4 , (16)
where L f h i , for 1 ≤ i ≤ 3 , is the healthy operating dynamics of the system expressed in [START_REF] Chen | Dynamic simulation of spur gear with tooth root crack propagating along tooth width and crack depth[END_REF].

Denote by χ the measurement vector:

χ =    y (ρ 1 ) 1 -L ρ 1 f h 1 . . . y (ρm) m -L ρm f h m    ∈ R m . (17)
Physically, χ ∈ R m represents the dierence between the estimated time-derivatives of the measurements y i and, the known dynamics of the healthy system. These derivatives can always be estimated if the system is observable with respect to the outputs y i = h i (x), which is guaranteed by the observability assumption below (see Assumption 3.1).

Observability is a measure of how well internal states of the system can be inferred from knowledge of its external outputs [START_REF] Hermann | Nonlinear controllability and observability In[END_REF]. A system is said to be observable if one can determine the behavior of the entire system from its measurements and eventually from the known inputs. More precisely, the current state values can be determined in nite-time using only the outputs (we send the reader to [START_REF] Hermann | Nonlinear controllability and observability In[END_REF] for dierent notions of observability and formal denitions). If a system is not observable, this means that the current values of some of its state variables cannot be determined through output sensors. So observability can been seen as the ability to create enough smart sensors in order to recover all states of the system. The denition of observability as well as the criteria for checking it uses the state-space representation. To system (2), we associate the coditribution

O(x) = span{dL j f h i (x), 1 ≤ i ≤ m, j ≥ 0}.
System (2) is locally observable with respect to its outputs

y i = h i (x), 1 ≤ i ≤ m, if locally dim O(x) = n, ( 18
)
where n is the state-space dimension. This condition is called observability rank condition.

Assumption 3.1 ((Observability)). We assume that system (2) is observable with respect to the outputs y i = h i (x), for 1 ≤ i ≤ m, that is, it veries the observability rank condition [START_REF] Dong | Hyperspectral image superresolution via non-negative structured sparse representation[END_REF].

If the system is observable, dynamical equations can be designed to estimate the state of the system from the outputs measurements. Such module is called a state observer or simply an observer of the system. A state observer is typically computer-implemented and allows to fully reconstruct the system state from its output measurements. If derivatives of the output are required (as in our case, where we need y ρ i 1 . . . y ρm m ), then the term dierentiator is used. Since Assumption 3.1 is veried, nonlinear dierentiators can be applied in order to estimate in nite-time the derivatives (y

(ρ 1 ) 1 , • • • , y (ρm)
m ), for more details, see [START_REF] Besançon | Nonlinear observers and applications[END_REF], and thus to compute χ which plays an important role in the reconstruction (based on the measurements y i and their derivatives) of the fault vector d(t). Indeed, from relation [START_REF] Derbel | Application of homogeneous observers with variable exponent to a mechatronic system[END_REF], it follows that the vectors χ and d(t) are related by an equation of the form:

χ = Φ(x)d(t), (19)
where Φ(x) is the (m × q)-matrix given by:

Φ(x) =    σ 11 . . .σ 1m L g 1 L ρ 1 -1 f h 1 . . . L gp L ρ 1 -1 f h 1 . . . . . . . . . . . . σ m1 . . .σ mm L g 1 L ρm-1 f h m . . . L gp L ρm-1 f h m    (20)
Example 3.1 (continuation). Based on ( 16), we deduce the measurement vector χ ∈ R 3 :

χ =   ẏ1 -L f h 1 (x) ẏ2 -L f h 2 (x) ẏ3 -L f h 3 (x)   , (21)
and the matrix Φ(x) ∈ R 3×7 :

Φ(x) =   x4 (1 -x4 ) 0 2x 4 (x 2 4 + x4 ) x4 x4 x4 -x 4 0 0 x4 x4 x4 1 -1 1 0 x4 0 2   . (22)
Equation ( 19) will be used to recover the unknown faults from the measurements information. Indeed, the goal of the fault detection problem is to solve the algebraic system of m equations, given by ( 19), with q unknowns, i.e., the components of the fault vector d(t) ∈ R q . This problem is well known as a left invertibility problem (see, for instance, [START_REF] Takahara | Structure of left invertibility problem of linear systems[END_REF], [START_REF] Nijmeijer | Right-invertibility for a class of nonlinear control systems: A geometric approach[END_REF], [START_REF] Respondek | Right and left invertibility of nonlinear control systems[END_REF]).

The left invertibility problem has been widely studied and applied to many dierent problems (going from fault diagnosis [START_REF] Martinez-Guerra | Fault diagnosis viewed as a left invertibility problem[END_REF] to cryptographic applications [START_REF] Vo Tan | Left invertibility, atness and identiability of switched linear dynamical systems: a framework for cryptographic applications[END_REF], etc). It consists in determining the causal factors (i.e., the unknown inputs signals or, like in our context, the fault vector) from the knowledge of a set of observations (i.e., the measurement outputs and its derivatives). In the case when q ≤ m (i.e., the number of measurements is greater or equal than the number of the unknown inputs corresponding to the faults), the inverse problem can be solved either by inverting, if possible, directly the matrix Φ(x) (when m = q and Φ(x) invertible) or by solving the system of equations based on the dimension reduction (if m > q). In the over complete matrix (when m < q, which is usually the case in the context of fault detection), the left invertibility problem either has no solution or has an innity of solutions. However, if the vector d is sparse, then, in the particular case when relation [START_REF] Du | Compressed-Sensing-Based Periodic Impulsive Feature Detection for Wind Turbine Systems In[END_REF] is static, that is, the matrix Φ and the fault vector d are constant (do not depend on time or state) reconstruction strategies enabling to exactly reconstruct the unknown signal and can be deployed (see, for example, [START_REF] Candes | Sparse signal and image recovery from compressive samples[END_REF], [START_REF] Yang | Image super-resolution via sparse representation[END_REF]).

In the context of fault detection, techniques based on sparse recovery methods have already been developed (see for example, [START_REF] Tang | Sparse classication of rotating machinery faults based on compressive sensing strategy[END_REF], [START_REF] Du | Compressed-Sensing-Based Periodic Impulsive Feature Detection for Wind Turbine Systems In[END_REF]) with the particularity that the faults are not introduced mathematically in the system and their features are recovered from the frequency content in which the most coecients are sparse (for instance, using Fourier or wavelet transforms). Therefore, the originality of our work is rstly, to apply the SRD method in the dynamical case (that is, when Φ(x) and d(t) depend on time or state) and secondly, to consider a system modeling based on a state-space representation that describes the healthy and the faulty behavior.

Another particularity of the SRD method is that all algorithms (corresponding to the optimization problem given by [START_REF] Kimmich | Fault detection for modern diesel engines using signal-and process model-based methods[END_REF] and to the observers equations, see, for instance [START_REF] García | Deterministic nonlinear observer-based approaches to fault diagnosis: A survey[END_REF], [START_REF] Besançon | Nonlinear observers and applications[END_REF], [START_REF] Derbel | Application of homogeneous observers with variable exponent to a mechatronic system[END_REF]) are solved in a continuous process (that is, in the same iteration loop) and with the same computation step. In other words, with the help of a single simulation of all algorithms, the diagnostic of numerous dierent faults is obtained which is not always possible for other diagnostic methods (that often require several analyzes to be able to diagnose dierent faults). Remark 3.1. One of the diculties when adapting the sparse recovery method to nonlinear dynamical systems is the fact that the matrix Φ(x) is state dependent and moreover, its choice is not always unique, as explained below. Indeed, very often, when modeling the system, the vector d(t) includes as much information as possible related to the faults, that is, the fault signatures w j (i.e., the quantity that is always zero when the operating system is healthy, and is nonzero when the corresponding fault acts on the system) as well as for each w j the common terms to all components of the vector g j through which the fault w j acts on the system. But d(t) can also be dened with the help of the signatures w j only (like in expression [START_REF] Candes | An introduction to compressive sampling[END_REF]). So for nonlinear dynamical systems, the choice of the fault vector d(t), and, therefore, of the matrix Φ(x) is not unique. The choice of Φ(x) is crucial for the SRD method. Indeed, Φ(x) has to satisfy some theoretical conditions (namely, the RIP condition that will be dened below, see Denition 3.4). If the matrix Φ(x) constructed in the classical way, does not verify that condition, some constant or variable coecients can be transferred from the vector d(t) to the matrix Φ(x) without loss of information on the fault vector (since we always keep in d(t) the fault signatures as explained above). 7 (notice that its 6th component contains the signature w 3 of the fault as well as the term x4 which is common to all g3 components, see (10)). Its associated matrix Φ(x) is given by:

Example 3.1 (continuation). Consider now the global fault vector

d(t)=[s 1 s 2 s 3 w 1 w 2 x4 w 3 w 4 ] T ∈ R
Φ =   x4 (1 -x4 ) 0 2x 4 (x 2 4 + x4 ) 1 x4 x4 -x 4 0 0 x4 1 x4 1 -1 1 0 x4 0 2   . ( 23
)
The goal of the SRD method is to reconstruct a sparse fault vector d(t) based on the measurement outputs y i , 1 ≤ i ≤ m, and their derivatives. The term sparse means that the dimension q of d(t) may be very large (i.e., there are many possible faults that are taken into account in the model), but only a few number of components may be non zero (i.e., only few faults act simultaneously on the system). In other words, the objective of this method is to reconstruct the most parsimonious solution of the fault vector d(t) with the help of the measurements vector χ and the matrix Φ(x) (i.e., where by most parsimonious, we mean the solution with the lowest number of non-zero components of d(t)). This can be reformulated as the following optimization problem ( [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF], [START_REF] Candes | An introduction to compressive sampling[END_REF], [START_REF] Becker | NESTA: A fast and accurate rst-order method for sparse recovery In[END_REF]):

min d∈R q d 0 , under the constraint χ = Φ(x)d(t), (24)
where d 0 is a pseudo-norm and corresponds to the number of non-zero values of d(t).

An equivalent description of problem ( 24) is to minimize a cost function constructed by leveraging the observation error χ -Φ(x)d(t) by respecting the parsimony constraint on d(t) via a balancing parameter λ (see, for example, [START_REF] Candes | An introduction to compressive sampling[END_REF]):

min d∈R q { 1 2 χ -Φ(x)d(t) 2 2 + λ d(t) 0 }. ( 25
)
The optimization problem [START_REF] García | Deterministic nonlinear observer-based approaches to fault diagnosis: A survey[END_REF] Notice that d(t) is s-sparse means that among the q = m+p modeled faults only s of them may act on the system. Recall that the number q of the total faults is much greater than s which in turn is smaller than m. For the SRD method, we need the following assumption: Assumption 3.2 (s-sparsity of the fault vector). We assume that d(t) is s-sparse and, moreover that s veries the following condition:

2s + 1 ≤ m. ( 26
)
Denition 3.4 (s-order RIP condition). A matrix Φ(x) is said to satisfy the s-order RIP condition, if for any s-sparse signal d(t), the following condition is veried:

(1 -δ s ) d 2 2 ≤ Φ(x)d(t) 2 2 ≤ (1 + δ s ) d 2 2 , (27)
where δ s ∈ (0, 1) is a constant parameter. This property will be called the s-order RIP condition.

In order to check if a matrix Φ(x) satises the s-order RIP condition, let Γ denote any set of s indices among 1 to q (where q is the total number of faults), i.e., Γ is any s-combination of the set {1, 2, . . . , q}. Notice that Γ is not unique and there are q! s!(q-s)! of them. Denote by Φ Γ (x) the sub-matrix of Φ(x) formed by the columns of Φ(x) indexed by the elements of the set Γ (i.e., corresponding to the active nodes s) and by Φ T Γ (x) its transposed matrix. The s-order RIP condition for a matrix Φ(x) is equivalent to the fact that all eigenvalues of the matrix product Φ T Γ (x)Φ Γ (x), for all possible s-combinations Γ, must be strictly between 0 and 2, that is:

1 -δ s ≤ eig(Φ T Γ (x)Φ Γ (x)) ≤ 1 + δ s (28)
for δ s ∈ (0, 1).

Remark 3.2. Recall that the matrix Φ(x) associated to system (2) and to fault vector [START_REF] Candes | An introduction to compressive sampling[END_REF] is not unique and its expression depends on the denition of d(t), see Remark 3.1. If no matrix Φ(x) (constructed using the procedure of Remark 3.1) veries the s-order RIP condition, we choose one of them and try to transform it in such a way that this property is satised for the transformed matrix. We explain next how to do it. First, notice that if the RIP condition is satised then all sub-matrices Φ Γ (x) are of full rank. Thus we have to eliminate the columns of the original matrix Φ(x) that are identical or proportional. In this case, their corresponding faults act in the same way on the system and they can be emerged by adding them and eliminating one of the columns. For instance, if the rst column Φ 1 of Φ(x) (corresponding to d 1 (t)) and the second one Φ 2 (corresponding to d 2 (t)) are like above (i.e., Φ 1 = Φ 2 or Φ 1 = aΦ 2 where a is a constant or a functional parameter), we will keep only one of them (say the rst one) and its new associated fault becomes

d 1 (t) + d 2 (t) if Φ 1 = Φ 2 , respectively, d 1 (t) + ad 2 (t), if Φ 1 = aΦ 2 .
Remark that by regrouping the faults, the dimension of the global fault vector d(t) decreases from q to q -1. After this procedure the dimensions of Φ(x) change as well (the number of columns decreases, but the number of lines remains unchanged) and the new matrix Φ(x) no longer contains identical or proportional columns (this does not mean that Φ(x) is of full rank, in fact, Φ(x) may possess columns that are linear combinations of at least two other columns. If the transformed matrix Φ(x) still does not verify the RIP condition, we start to eliminate those collinear columns and apply the same principle to regroup the corresponding faults (for instance, if Φ 1 = aΦ 2 + bΦ 3 , where a and b are constant or functional parameters, we keep Φ 2 and Φ 3 , but eliminate Φ 1 , the new faults corresponding to Φ 2 and Φ 3 being d 2 (t) + ad 1 (t) and d 3 (t) + bd 1 (t)). After each such column removal, we check if the new matrix Φ(x) verify the s-order RIP condition.

In order to check the s-order RIP condition, it may also be helpful to normalize each column of the matrix Φ(x) with respect to the norm • 2 (after normalization, we should have

Φ k (x) 2 = m i=1 |Φ ik (x)| 2 = 1, 1 ≤ k ≤ q,
where Φ k denotes the kth column of Φ). As a consequence of this normalization, each component of the fault vector has to be multiplied by the norm of its corresponding column in Φ(x). In some cases, the normalization of the matrix Φ(x) is not necessary to verify the s-order RIP condition (there are matrices verifying this condition whose columns are not normalized).

In the sequel, we continue to denote the transformed matrix by Φ(x), the new fault vector by d(t) and its number components by q. Suppose that: Assumption 3.3. (s-order RIP) The transformed matrix Φ(x) is assumed to satisfy the s-order RIP condition. Under Assumption 3.3, the parsimony problem [START_REF] García | Deterministic nonlinear observer-based approaches to fault diagnosis: A survey[END_REF] for the transformed matrix Φ(x) and the new fault vector d(t) becomes:

min d∈R q { 1 2 χ -Φ(x)d(t) 2 2 + λ d(t) 1 }, (29)
and its solution exists and is unique [START_REF] Yu | Dynamical sparse recovery with nite-time convergence[END_REF].

Recall that m < q (i.e., the number of monitoring sensors is signicantly smaller than the number of faults) and in general, even after transforming the original matrix Φ(x), when necessary, we still have m < q, where q denotes now the dimension of the new regrouped fault vector, thus, Φ(x) cannot be invertible (since the matrix is not square) and, the left invertibility cannot be applied. Nevertheless, under Assumptions 3.2 and 3.3, the fault vector d(t) can be detected and reconstructed as a solution of (29). To that end, the following dynamical algorithm based on sliding mode techniques [START_REF] Edwards | Sliding mode control: theory and applications[END_REF] is proposed reminding very much that of [START_REF] Yu | Dynamical sparse recovery with nite-time convergence[END_REF] with the dierence that the matrix Φ and the vector d are no longer constant (see also [START_REF] Balavoine | Convergence speed of a dynamical system for sparse recovery[END_REF] for the linear case):

τ u(t) = -u(t) + (Φ T (x)Φ(x) -Id) d(t) -Φ T (x)χ α , d(t) = ϕ λ (u(t)), (30)
where . α = |.| α sign(.), u(t) ∈ R q is the internal state vector, ϕ λ (u(t)) = max(|u|λ, 0)sign(u) is the continuous soft thresholding function, λ ∈ R q is a vector of constant positive parameters that has to be suitably chosen in function of the noise and the minimum possible absolute values of the faults, the vector d(t) represents the estimation of the sparse signal d(t), τ is a (q × q)-diagonal matrix with constant parameters τ i , 1 ≤ i ≤ q, determined by the physical properties of the implementing system, Id is the identity (q × q)-matrix and the exponential coecient α is such that α ∈ [0, 1]. In [START_REF] Yu | Dynamical sparse recovery with nite-time convergence[END_REF], it has been proven, in the case when the matrix Φ(x) and the vector d(t) are constant, that under Assumptions 3.2 and 3.3 and for α ∈ [0, 1), the state u converges, in nite-time, to its equilibrium point and thus the estimated fault vector d(t) converges in nite-time to d(t). For α = 1, only asymptotic convergence is ensured. If Φ(x) and d(t) depend explicitly on time or/and states, when adapting the proof of [START_REF] Yu | Dynamical sparse recovery with nite-time convergence[END_REF] to this dynamical case, the derivatives of Φ(x) and d(t) appear and have to be crushed with a suitable choice of λ. In this case, a practical stability can be guaranteed (see [START_REF] Chaillet | Uniform semiglobal practical asymptotic stability for non-autonomous cascaded systems and applications[END_REF] for more details).

Remark 3.3. In the case when dierent faults are regrouped, the estimated fault vector and d 2 (t) is activated, in general, the faults diagnosis is possible either by verifying the order of magnitude of the defect or by the type of the obtained fault signal. However, if the two added faults are activated simultaneously, the estimated fault vector given by the SRD method contains the sum of these two faults and we may not be able to identify or to localize them. These situations will be discussed in our case study in the next section.

d(t)
Remark 3.4. In some cases, the dynamics f and/or vector elds g j can depend explicitly on time (i.e., we have f = f (t, x(t)) and/or g j = g j (t, x(t)), for some 1 ≤ j ≤ p). Thus the system is a time-variant one. A classic way to transform it into a time-invariant system (that is, of form ( 2)) is to introduce an additional state (which is simply the time-variable) described by the equation:

ẋn+1 (t) = 1. (31)
The associated extended system is time-invariant and contains n+1 states variables x x n+1 ∈ R n+1 . The new dynamics is f (x, x n+1 ) 1 and the vector elds associated to the faults are g j (x, x n+1 ) 0 , for 1 ≤ j ≤ p. Moreover t is, in general, also an output (or, at least, a known state) of the system.

Case Study

F 1ec Eccentricity fault vector [N × m] F c Crack fault vector [N × m]
In this section, we will apply the SRD method to diagnose mechanical faults in a gear power transmission system. The considered system is a reduced model of two-stage gear that takes into account only four degrees of freedom (see Figure 1). The principal notations used in the system modeling are recalled in Table 2.

Healthy model

The healthy system (see Figure 1) is composed of two pinions and two wheels supported by three shafts, one for the input shaft (the motor), the second for the intermediate shaft (which connects both stages of the gearbox system) and the third for the output shaft (the load). The model is obtained by developing the Lagrange equations of the kinetic and potential energies corresponding to the two-stage gear (see also [START_REF] Fakhfakh | Eect of manufacturing and assembly defects on two-stage gear systems vibration[END_REF] and [START_REF] Derbel | Electromechanical system control based on observers[END_REF] for a single stage modelisation).

The Lagrange formalism leads to the set of dierential equations governing the system motion:

Mq + Cq + K(t)q = F, ( 32
)
where q is the vector of degrees of freedom q = [θ 1 θ 2 θ 3 θ 4 ] t , M is the mass matrix expressed by:

M =     I 1 0 0 0 0 I 2 + I s 0 0 0 0 I 3 + I s 0 0 0 0 I 4     , (33)
with I i the inertia moment of the ith gear, 1 ≤ i ≤ 4, and I s the inertia moment of the intermediate shaft, C the damping, F the applied forces vector F = [C m 0 0 C r ] T , where C r and C m are, respectively, the load and the motor torques of the system and K(t) presents the time-varying stiness matrix containing the gear-mesh stiness matrix K e (t) and the constant shaft torsional stiness matrix K s :

K(t) = K e (t) + K s . (34)
The gear-mesh stiness matrix K e (t) is dened by :

K e (t) =    k 1 (t)c 2 R 2 1 k 1 (t)c 2 R 1 R 2 0 0 k 1 (t)c 2 R 1 R 2 k 1 (t)c 2 R 2 2 0 0 0 0 k 2 (t)c 2 R 2 3 k 2 (t)c 2 R 3 R 4 0 0 k 2 (t)c 2 R 4 R 3 k 2 (t)c 2 R 2 4    , (35)
where R i , 1 ≤ i ≤ 4, represents the base radius of the ith gear and k j (t), j = 1, 2, is the time-varying stiness for the jth stage dened as a square wave depending on time. The constant c simply denotes c = cos(β), with β the helix angle of the gears.

The shaft torsional stiness matrix K s is expressed as:

K s =     0 0 0 0 0 k s -k s 0 0 -k s k s 0 0 0 0 0     , (36)
where k s is the constant torsional stiness of the shaft. Developing [START_REF] Lin | Impulse detection using a shift-invariant dictionary and multiple compressions[END_REF], we obtain the following state-space representation of the two-stage gearbox model:

                                 θ1 = Ω 1 , θ2 = Ω 2 , θ3 = Ω 3 , θ4 = Ω 4 , Ω1 = 1 I 1 [-K 11 θ 1 -K 12 θ 2 -C 11 Ω 1 -C 12 Ω 2 + C m ], Ω2 = 1 I 2 +Is [-K 21 θ 1 -K 22 θ 2 -K 23 θ 3 -C 21 Ω 1 -C 22 Ω 2 -C 23 Ω 3 ], Ω3 = 1 I 3 +Is [-K 32 θ 2 -K 33 θ 3 -K 34 θ 4 -C 32 Ω 2 -C 33 Ω 3 -C 34 Ω 4 ], Ω4 = 1 I 4 [-K 43 θ 3 -K 44 θ 4 -C 43 Ω 3 -C 44 Ω 4 + C r ], (37)
where K ij and C ij , for 1 ≤ i, j ≤ 4, are respectively, the coecients of the stiness matrix K(t) and damping matrix C, [θ Ω] T ∈ R 8 is the state vector whose components are the degrees of freedom of the system. For (37), we have the following measurements:

y 1 = θ 1 , y 2 = θ 2 , y 3 = θ 3 , y 4 = θ 4 . (38)

Faulty model

In this paper, we suppose that mechanical faults and sensor faults can aect the gearbox system. We model two types of mechanical faults: eccentricity and crack defects. The faulty model is given by the following equation of motion The eccentricity fault is dened as the gap between the geometrical and the rotational axes of the gear. It is composed by two components: the rst one F 1ec is related to a modication in the kinetic energy and the second one F 2ec is due to the change in potential energy:

Mq + C q + K(t)q = F + F ec + F c , (39)
F ec = F 1ec + F 2ec . (40)
We have:

-

F 1ec = [0 I 2 ω2 0 I 4 ω4 ] T
, where ω i is the angular velocity of the ith gear given by

ω i = -R i-1 R i ω i-1 + ėi R i cos(β)
, with e i = cos(β)[sin(ω i-1 t -φ i-1 ) i-1 -sin(ω i t -φ i ) i ], the deviation due to the eccentricity defect of the gear, φ i the initial eccentricity phases and i the eccentricity magnitude, i ∈ {2, 4},

-F 2ec = [k 1 (t)δeR 1 cos(β) k 1 (t)δeR 2 cos(β) k 2 (t)δe R 3 cos(β) k 2 (t)δeR 4 cos(β)] T , with
δe the initial deviation due to the eccentricity.

In our case, without loss of generality, the eccentricity eect due to the potential energy F 2ec is neglected, since the eccentricity fault aects the gear system on the inertia of the system (so, in our model, its impact is already taken into account in the expressions of ωi and we keep the term due to the kinetic energy F 1ec only).

The crack defect F c is given by

F c = [R 1 cos(β)∆k 1p (t) R 2 cos(β)∆k 1w (t) R 3 cos(β) ∆k 2p (t) R 4 cos(β)∆k 2w (t)] T .
The term ∆k jp (t) (respectively, ∆k jw (t)) characterizes the variation of the time-varying gear-mesh stiness k j (t) associated to the pinion (respectively to the wheel) of the jth stage of the gearbox system (in our application, j ∈ {1, 2}). These quantities are indeed the product between the time-varying gear-mesh variation and the deviation due to the crack. The fault vector F c will be indicated through the drop of the time-varying stiness of the gears k j (t), j ∈ {1, 2}. In the healthy functioning case of the gear (without faults), the global time-varying stiness is a square wave function that normally evolves, for a spur gear, between a maximum value that corresponds to the maximum number of segments in contact (two pairs of teeth) and a minimum value corresponding to the minimum number of segments in contact (one pair of teeth in contact), see Figure 2a, where t m1 is the meshing gear period for the rst stage of the gearbox system. The descent in the curves (see Figure 2a) denes the progressive exit of a tooth pair from the meshing zone and the rise corresponds to the gradual entry of a new pair of teeth in contact. In the presence of defects, the disturbances appear on the temporal evolution of the global timevarying stiness k 1 (t) (see Figure 2b, for the rst stage) at the moments of the passage of the defective tooth (here, Z 1 = 18) in the zone of contact (see, e.g., [START_REF] Chen | Dynamic simulation of spur gear with tooth root crack propagating along tooth width and crack depth[END_REF], [START_REF] Yu | The eects of spur gear tooth spatial crack propagation on gear mesh stiness[END_REF]).

In the faults modeling, we also take into account constant sensor faults ∆θ 1 , ∆θ 2 , ∆θ 3 and ∆θ 4 (corresponding, for instance, to oset faults), related to the two gear displacements. The faulty model considering the crack, eccentricity and sensor faults is given by the following system (compare it to [START_REF] Takahara | Structure of left invertibility problem of linear systems[END_REF], the healthy model):

                                         θ1 = Ω 1 , θ2 = Ω 2 , θ3 = Ω 3 , θ4 = Ω 4 , Ω1 = 1 I 1 [-K 11 θ 1 -K 12 θ 2 -C 11 Ω 1 -C 12 Ω 2 + C m +R 1 cos(β)∆k 1p (t)], Ω2 = 1 I 2 +Is [-K 21 θ 1 -K 22 θ 2 -K 23 θ 3 -C 21 Ω 1 -C 22 Ω 2 -C 23 Ω 3 + R 2 cos(β)∆k 1w (t) + I 2 ω2 ], Ω3 = 1 I 3 +Is [-K 32 θ 2 -K 33 θ 3 -K 34 θ 4 -C 32 Ω 2 -C 33 Ω 3 -C 34 Ω 4 + R 3 cos(β)∆k 2p (t)], Ω4 = 1 I 4 [-K 43 θ 3 -K 44 θ 4 -C 43 Ω 3 -C 44 Ω 4 + C r +R 4 cos(β)∆k 2w (t) + I 4 ω4 ], (41)
with the measurements:

y 1 = θ 1 + ∆θ 1 , y 2 = θ 2 + ∆θ 2 , y 3 = θ 3 + ∆θ 3 , y 4 = θ 4 + ∆θ 4 . ( 42
)
The appearance of these faults in the gearbox system modify directly the dynamical behavior of the system. It must be noted that the faulty system ( 41) is a non-autonomous system that depends explicitly on time (i.e., the drift and the vector elds associated to the faults involve explicitly the t-variable through the terms K ij , 1 ≤ i, j ≤ 4, ∆k jp (t) and ∆k jw (t), for j ∈ {1, 2}). In order to transform it into form [START_REF] Candes | Sparse signal and image recovery from compressive samples[END_REF], the state extension explained in Remark 3.4 is applied, that is, we add to system (41) the following equation ṫ = 1. From now on, when we say system (41), we will actually refer to the extended system associated to [START_REF] Yang | Image super-resolution via sparse representation[END_REF]. Moreover, notice that in system [START_REF] Takahara | Structure of left invertibility problem of linear systems[END_REF], the measurements correspond indeed to states of the system y i = θ i , for 1 ≤ i ≤ 4. Thus, we do not need to make the change of coordinates explained in the second section since system (41) is already in form [START_REF] Candes | Sparse signal and image recovery from compressive samples[END_REF]. We will denote the faulty dynamics of (41), i.e., the right-hand side of [START_REF] Yang | Image super-resolution via sparse representation[END_REF], by F (θ, Ω).

The total number of possible faults is much larger than 4 (the number of the monitoring sensors, see [START_REF] Tang | Sparse classication of rotating machinery faults based on compressive sensing strategy[END_REF]). Therefore, the SRD method can be applied. In the next part, we will apply it in order to detect, identify and locate possible faults.

Application of the sparse recovery method

We dene the global fault vector d(t) that contains, not only the faults signatures but as mush information as possible (see Remark 3.1):

d(t) =                   ∆θ 1 ∆θ 2 ∆θ 3 ∆θ 4 R 1 cos(β)∆k 1p (t) R 2 cos(β)∆k 1w (t) R 3 cos(β)∆k 2p (t) R 4 cos(β)∆k 2w (t) I 2 ω2 I 4 ω4                   ∈ R 10 . ( 43
)
In order to get relation [START_REF] Du | Compressed-Sensing-Based Periodic Impulsive Feature Detection for Wind Turbine Systems In[END_REF] between the measurement vector χ and the fault vector d(t), we rst compute the relative degrees of the outputs y 1 = θ 1 + ∆θ 1 , y 2 = θ 2 + ∆θ 2 , y 3 = θ 3 + ∆θ 3 and y 4 = θ 4 + ∆θ 4 as presented in (13) (recall that the sensor faults are supposed piecewise-constant). For instance, consider y 1 (the same calculation has to be done for y 2 , y 3 and y 4 ):

   y 1 = θ 1 + ∆θ 1 , ẏ1 = θ1 = Ω 1 , ÿ1 = Ω1 = F 5 (θ 1 , θ 2 , Ω 1 , Ω 2 ), (44)
where

F 5 = 1 I 1 [-K 11 θ 1 -K 12 θ 2 -C 11 Ω 1 -C 12 Ω 2 + C m + R 1 cos(β)∆k 1p (t)]
describes the faulty dynamics of Ω 1 , see [START_REF] Yang | Image super-resolution via sparse representation[END_REF].

In order to extract the sensor faults, we replace the states θ 1 and θ 2 in

F 5 (θ 1 , θ 2 , Ω 1 , Ω 2 )
by their measured values θ 1 + ∆θ 1 and θ 2 + ∆θ 2 (i.e., with the presence of the sensor faults):

ÿ1 = F 5 (θ 1 + ∆θ 1 , θ 2 + ∆θ 2 , Ω 1 , Ω 2 ), (45)
= 1 I 1 [-K 11 (θ 1 + ∆θ 1 ) -K 12 (θ 2 + ∆θ 2 ) -C 11 Ω 1 -C 12 Ω 2 + C m + R 1 cos(β)∆k 1p (t)]. (46)
Applying the same principle for the remaining outputs, we obtain:

ÿ2 = F 6 (θ 1 + ∆θ 1 , θ 2 + ∆θ 2 , θ 3 + ∆θ 3 , Ω 1 , Ω 2 , Ω 3 ), ( 47
) ÿ3 = F 7 (θ 2 + ∆θ 2 , θ 3 + ∆θ 3 , θ 4 + ∆θ 4 , Ω 2 , Ω 3 , Ω 4 ), (48)
ÿ4 = F 8 (θ 3 + ∆θ 3 , θ 4 + ∆θ 4 , Ω 3 , Ω 4 ), (49)
where F 6 , F 7 and F 8 describe, respectively, the faulty dynamics of Ω 2 , Ω 3 and Ω 4 , see [START_REF] Yang | Image super-resolution via sparse representation[END_REF] for their expressions. When developing, ( 46)-(49) become:

ÿ1 = L 2 f h 1 - K 11 I 1 ∆θ 1 - K 12 I 1 ∆θ 2 + R 1 cos(β)∆k 1p (t) I 1 , (50) ÿ2
= L 2 f h 2 - K 21 I 2 + I s ∆θ 1 - K 22 I 2 + I s ∆θ 2 - K 23 I 2 + I s ∆θ 3 + R 2 cos(β)∆k 1w (t) I 2 + I s + I 2 ω2 I 2 + I s , (51)
ÿ3 = L 2 f h 3 - K 32 I 3 + I s ∆θ 2 - K 33 I 3 + I s ∆θ 3 - K 34 I 3 + I s ∆θ 4 + R 3 cos(β)∆k 2p (t) I 3 + I s , (52) ÿ4
= L 2 f h 4 - K 43 I 4 ∆θ 3 - K 44 I 4 ∆θ 4 + R 4 cos(β)∆k 2w (t) I 4 + ω4 . (53)
Notice that the above relations are ane with respect to the sensor faults ∆θ i and they are of form [START_REF] Derbel | Application of homogeneous observers with variable exponent to a mechatronic system[END_REF].

Following [START_REF] Dong | A survey on distributed ltering and fault detection for sensor networks[END_REF], we deduce the measurement vector χ containing the available information about the system (recall that χ denotes the dierence between the derivatives of the outputs and the healthy operating dynamics of the system):

χ =         ÿ1 -1 I 1 [-K 11 θ 1 -K 12 θ 2 -C 11 Ω 1 -C 12 Ω 2 + C m ] ÿ2 -1 I 2 +Is [-K 21 θ 1 -K 22 θ 2 -K 23 θ 3 -C 21 Ω 1 -C 22 Ω 2 -C 23 Ω 3 ] ÿ3 -1 I 3 +Is [-K 32 θ 2 -K 33 θ 3 -K 34 θ 4 -C 32 Ω 2 -C 33 Ω 3 -C 34 Ω 4 ] ÿ4 -1 I 4 [-K 43 θ 3 -K 44 θ 4 -C 43 Ω 3 -C 44 Ω 4 + C r ]         , (54)
where ÿ1 , ÿ2 , ÿ3 and ÿ4 are estimated using dierentiators (see [15] where a sliding mode dierentiator is applied to a single stage of gears). It must be noted that dierentiators can be implemented since system [START_REF] Yang | Image super-resolution via sparse representation[END_REF] veries the observability property explained in Assumption 3.1. Algorithms very similar to those of [START_REF] Derbel | Application of homogeneous observers with variable exponent to a mechatronic system[END_REF] are used for the two-stage gear (in particular, the dierentiators parameters that have to be suitably chosen in order to ensure rapid and accurate convergence between the real and the estimated system states, have the same order of magnitude as those of [START_REF] Derbel | Application of homogeneous observers with variable exponent to a mechatronic system[END_REF]) and we do not present them here.

From this and denition (43) of d(t), we determine the matrix Φ(x) linking χ and d(t) as in [START_REF] Du | Compressed-Sensing-Based Periodic Impulsive Feature Detection for Wind Turbine Systems In[END_REF]:

Φ =      -K 11 I 1 -K 12 I 1 0 0 1 I 1 0 0 0 0 0 -K 21 I 2 +Is -K 22 I 2 +Is -K 23 I 2 +Is 0 0 1 I 2 +Is 0 0 1 0 0 -K 32 I 3 +Is -K 33 I 3 +Is -K 34 I 3 +Is 0 0 1 I 3 +Is 0 0 0 0 0 -K 43 I 4 -K 44 I 4 0 0 0 1 I 4 0 1 I 4      . ( 55
)
An important observation is that in our application, the matrix Φ(x) depends explicitly on time through the K-coecients.

Recall that the SRD method guarantees that at most s faults can be accurately detected and according to Assumption 3.2, the integer s should verify 2s + 1 ≤ m, where m is the number of the monitoring sensors (here, we have m = 4). Thus, the SRD method guarantees that at most one fault can be accurately detected for the gearbox model (however, we will see in the next section that, in general, we are able to accurately detect more than one fault). This point can be seen as a limitation of the method but it is justied by the fact that in the industries (and in particular for mechanical systems), the presence of several simultaneous faults is not commun because they generate large disturbances in the system and leads to the extreme case when the machine stops.

One of the most important conditions of the SRD method is that Φ(x) has to satisfy the s-order RIP property (here, with s = 1), see Denition 3.4. From (55), it can be noticed that some columns are proportional (i.e., the faults corresponding to these columns, here, the sixth and the ninth columns of Φ(x), respectively, the eighth and the tenth columns of Φ(x), act in the same way on the system). Therefore, we deduce immediately that Φ(x), given by (55), does not verify the s-order RIP condition. Moreover, even if we dene the global fault vector d(t) with the help of the faults signatures only (that is, in [START_REF] Yu | The eects of spur gear tooth spatial crack propagation on gear mesh stiness[END_REF], the terms R i cos(β), for 1 ≤ i ≤ 4, I 2 and I 4 are absent), its associate matrix Φ(x) would not satisfy the s-order RIP property either. Finally, observe that not all columns of the matrix Φ(x) are normalized with respect to the 2-norm. Thus, we have to modify the fault vector d(t), as explained in Remark 3.2, by multiplying the faults associated to the non normalized components of Φ(x) by their corresponding norm and regrouping the defects which generate the collinearity of the columns:

d(t) =                n 1 ∆θ 1 n 2 ∆θ 2 n 3 ∆θ 3 n 4 ∆θ 4 R 1 cos(β)∆k 1p (t) I 1 R 2 cos(β)∆k 1w (t) I 2 +Is + I 2 ω2 I 2 +Is R 3 cos(β)∆k 2p (t) I 3 +Is R 4 cos(β)∆k 2w (t) I 4 + ω4                ∈ R 8 , (56)
with

n 1 = Φ 1 (x) 2 , n 2 = Φ 2 (x) 2 , n 3 = Φ 3 (x) 2 and n 4 = Φ 4 (x) 2 , where Φ i (x)
denotes the ith column of (55) and n i is its 2-norm.

In this case, the normalization of the matrix Φ(x) is required since it contains dierent orders of magnitude. By normalizing Φ(x), we do not lose the fault vector information and the new matrix (after regrouping the faults and normalization) becomes:

Φ =        -K 11 I 1 -K 12 I 1 n 2 0 0 1 0 0 0 -K 21 (I 2 +Is)n 1 -K 22 (I 2 +Is)n 2 -K 23 (I 2 +Is)n 3 0 0 1 0 0 0 -K 32 (I 3 +Is)n 2 -K 33 (I 3 +Is)n 3 -K 34 (I 3 +Is)n 4 0 0 1 0 0 0 -K 43 I 4 n 3 -K 44 I 4 n 4 0 0 0 1        . (57)
It must be noted that with the new representation, we remain in the over complete matrix case (i.e., the number of monitoring sensors is smaller than the number of faults, m = 4, q = 8). The next section presents the simulations results obtained for dierent activated faults as well as for dierent cases of s. We will rst see that, in accordance with the theoretical results, when only one fault acts on the system, that fault is always correctly diagnosed by the SRD method. Second, we will see that the SRD method works also in many cases when two faults are activated.

Simulations and results

In order to evaluate the performances of the SRD method, simulations are carried on Matlab/Simulink. The system parameters are given in Table 3. The designed system was implemented with an ode1 Euler (we choose this solver because system (30) is not C 1 -smooth), and with the exponential coecient α of (30) being α = 0.5 (assuring a convergence in nitetime). Notice that the considered faults have very dierent order of magnitude and we have to take this into account when choosing the algorithm parameters λ and τ : λ i = 10 -2 , τ i = 10 4 for 1 ≤ i ≤ 4 (associated to sensor faults) and λ i = 10 7 , τ i = 10 -10 for 5 ≤ i ≤ 8 (associated to operating faults). This choice of parameters is due to in the system modeling, that is, we consider dierent types of faults with dierent orders of magnitude. Notice that if all faults have the same order of magnitude, the same parameters λ and τ can be used for all of them, see [START_REF] Torki | A novel FDI sparse recovery method: Application on PMSG wind turbine[END_REF]). Simulations results show the good performances of the SRD method to diagnose several mechanical faults in the gearbox system developed in the previous section. Scenario 1 presents 662 rpm the case when only one sensor fault occurs, i.e., s = 1, (recall that, according to Assumption 3, we can accurately detect one fault among the eight possible). Scenario 2 explains the case when crack faults occur. Scenario 3 considers two merged faults (crack and eccentricity) that are activated (notice that we are still in the case s = 1, since for the modied fault vector, the two original faults are merged into a single one). Scenario 4 illustrates several examples when two dierent faults occur simultaneously (i.e., s = 2). Scenario 5 shows the application of the SRD method when the measurements of the system are aected by a random noise. The non activated faults are set to zero and for all scenarios the corresponding estimates converge to zero. The last scenario (scenario 6) shows the application of SRD method with non-stationary operating conditions. For each scenario, the gures show all real defects d i (t), 1 ≤ i ≤ 8, and the estimated ones di (t), 1 ≤ i ≤ 8. Presenting all curves on the same gure allows us to highlight the sparsity of the fault vector as well as the fact that (after a short transient time) the estimates of the non-activated faults (those that are identically zero and do not act on the system) converge indeed to zero. Zoomed regions of the estimated faults are carried out in all scenarios in order to show the transient parts. We notice that at the beginning of the algorithm, all estimated faults di (t), 1 ≤ i ≤ 8, varies in time, then the non-activated ones converge to zero. The bold dotted signals represents the activated faults while the bold-solid line signal correspond to their estimation given by the SRD method. This scenario illustrates the case when one fault among the eight possible faults given by ( 56) is activated. The SRD method is applied four times to extract each of the following faults (for each case of Figure 4a-3d only one of them is activated and the others are set to zero):

s 1 = ∆θ 1 = 0, 0 ≤ t < 0.01s, 10 -2 , t ≥ 0.01s, (58)
s 2 = ∆θ 2 = 0, 0 ≤ t < 0.01s, 210 -2 , t ≥ 0.01s, (59)
that aect the system at t = 0.01,

s 3 = ∆θ 3 = 0, 0 ≤ t < 0.02s, 510 -3 , t ≥ 0.02s, (60)
that aects the system at t = 0.02s and

s 4 = ∆θ 4 = 0, 0 ≤ t < 0.015s, 10 -3 , t ≥ 0.015s, (61)
that aects the system at t = 0.015s. Figure 4a-3d show that the activated fault is always accurately recovered in nite-time by the SRD method, while the estimated values of the non activated faults converge to zero. Indeed, from each sub-gure, it is clear that the fault signal di related to the sensor fault of θ i , 1 ≤ i ≤ 4, converges to its non zero value, while the remaining di converge to zero. This conrms the good performances of the SRD method to diagnose sensor faults. This scenario presents the case when either the fth or the seventh fault of d(t) is activated. Both of them are operating faults (corresponding to the presence of the crack in the rst or in the second pinion of the gearbox system). It is clear from Figure ?? that the estimated fault d5 recovered by the SRD method converges in nite-time to the acti- vated fault and the non activated ones converge to zero. The signal d5 is characterized by a time-varying stiness variation that occurs at each contact of the defective gear and its identication is guaranteed by the signal type (square wave that corresponds to a time-varying stiness form). Its localization is obtained by the distance between two successive peaks of the time-varying stiness variation that corresponds to the teeth number of the defective gear. Here, this distance equals to Z 1 = 18 in Figure ?? and Z 3 = 20 in Figure 4b. The teeth numbers of the defective pinion i, for i ∈ {1, 3} are characteristics of the system which are given in Table 3.

Scenario 3 (s = 1, an activated merged fault corresponding to eccentricity and crack defects)

This scenario corresponds to the case when either the defect d 6 or d 8 , regrouping the original crack and the eccentricity faults of the gear, aect the system. We present two cases (for the rst one d 6 is activated and d 8 = 0, while for the second one d 6 = 0 and d 8 is activated). From Figure 5, it is clear that we detect the presence of these faults by the convergence of the estimated defects values to their real values. The distinction between crack and eccentricity 3) while the appearance of the eccentricity fault is modeled by a sinusoidal form signal (see Figure 5b and5d). The localization of the eccentricity fault is determined from the rotational speed of the defective gear. For instance, in Figure 5b, the rotational speed of the obtained signal is w = 2πf = 2π T = 2π 0.1165-0.0807 = 175.5 rad/s which corresponds to the rotational speed of the rst pinion (see Table 3).

In this scenario, a complete localization has been carried out: rstly, we are able to detect the presence of the merged crack and eccentricity faults, secondly, to identify and distinguish them from the obtained signals given by the SRD method, and, nally, to localize them with the help of some characteristics of the mechanical system presented in the obtained signals (e.g., the teeth number, the rotational speed of wheel/pinion). Recall the eciency of the SRD method is guaranteed for 2s + 1 ≤ m (in the case of the gearbox application, when only a single fault is activated). A natural question is what happens when at least two faults occur simultaneously. While the previous scenarios present the SRD method results for one active fault, Scenario 4 treats the case when two dierent operating faults appear simultaneously (see Figure 6)). Figure 6a illustrates the presence, at the same time, of the crack defect in the rst pinion and in the rst wheel. The SRD method gives a good convergence of the estimated faults to their real values and a complete diagnostic can be carried out. The fault distinction is always possible by the tooth number of the defective gear.

Figure 6b displays the case when two sensor faults related to the displacements θ 1 and θ 2 are activated. We can notice that these activated faults are also accurately recovered in nitetime via the SRD method while the estimated quantities of the non activated faults converge to zero. Thus, we are able to detect the presence of two dierent faults simultaneously and conclude that the condition 2s + 1 ≤ m given in Assumption 3.2 is sucient but not necessary for an accurate reconstruction of the faults using the SRD method.

Figure 6c presents the case when a crack in the rst pinion and a sensor fault are activated. It is clear that these faults are not recovered by the SRD method and their diagnostic is not possible. This is due to the fact that the faults have dierent orders of magnitude. It should be however noted to that our diagnosis method always allows the detection of an anomaly in the system by the non-zero signals given by the SRD method.

Obviously it is possible to diagnose more defects by increasing the number of sensors and therefore, providing more information about the system, and allowing to take into account more possible faults. However, the simultaneous presence of several defects in industrial systems leads to a fast decrease of the equipment life and to a large perturbation or even the stop of the operating system. Usually, the collected sensors signals are polluted by random noises due to the presence of many source of vibration. Scenario 5 considers the case when the measurements are aected by a white random noises. Figure 7a shows the presence of a sensor fault, which is ∆θ 1 , while in Figure 7b a crack fault in the rst pinion occurs. For both cases, the estimated faults given by the SRD method converge in nite-time to their real values. The diagnostic of these faults is always possible by analyzing the nature of the obtained signals (as explained in the previous scenarios). Thus, this diagnosis method is still able to identify mechanical and sensor faults by including random measurement noises in the dynamical state-space representation.

Scenario 6 (non-stationary operating conditions)

The last scenario illustrates the case of non-stationary operating conditions of the gearbox system which is subject to a variable resistant torque. We have considered a load uctuating in a saw-tooth shape (see Figure 41)).

Figure 9 shows the time evolution of the acceleration of the rst stage gear by applying a variable load torque. We notice the presence of peaks that correspond to the peaks of the applied variable load. A sensor fault on the rst pinion is activated at t = 0, 01s. Figure 10 illustrates the good performances of the SRD method in non-stationary conditions. [START_REF] Candes | Sparse signal and image recovery from compressive samples[END_REF]

Conclusions

In this paper, a diagnosis method, called sparse recovery diagnosis (SRD), for nonlinear dynamical systems is presented. Then it is applied to the mechanical system of a twostage gear model. A healthy and faulty modeling of the considered system are presented and implemented via Matlab/Simulink. The SRD is applied in order to diagnose several mechanical faults. Eccentricity, crack and sensor faults are modeled and recovered in nitetime by this diagnosis method. The diagnosis method introduced in this paper is an eective diagnosis method but it is not always able to detect and distinguish all types of defects. For this reason, combining dierent diagnosis methods is suggested to obtain an optimized and a perfect diagnostic for mechatronic systems. Our future work will be focused on the comparison between other diagnosis techniques and the SRD algorithm as well as on the theoretical generalization of the SRD method for a class of mechanical systems with the presence of dierent possible faults.

nL

  Dimension of the system state m Number of monitoring sensors p Number of system faults q Total number of faults (q = p + m) ρ iRelative degree of the output y i f h i th-order Lie derivative of the function h i along the vector eld f

Figure 1 :

 1 Figure 1: Torsional model of a two-stage gear

Figure 2 :

 2 Figure 2: Time-varying stiness evolution divided by the meshing gear period t m1

Figure 3 :

 3 Figure 3: Scenario 1: Sensor faults

5. 2

 2 Scenario 2 (s = 1, an activated operating fault) Crack fault for the second pinion: second stage

Figure 4 :

 4 Figure 4: Scenario 2: Crack faults

Figure 5 :

 5 Figure 5: Scenario 3: Merged faults

5. 4

 4 Scenario 4 (s = 2, two faults occur simultaneously) Two activated faults: eccentricity and sensor faults

Figure 6 : 4 :

 64 Figure 6: Scenario 4: Two activated faults s = 2

5. 5

 5 Scenario 5 (s = 1, an activated fault in the presence of noises) Crack fault of rst pinion

Figure 7 : 5 :

 75 Figure 7: Scenario 5: Presence of noise

2 )Figure 9 : 3 Figure 10 :

 29310 Figure 8: Time evolution of load

Table 1 :

 1 Principal notations

	2 Sparse recovery diagnosis method
	3 Problem statement
	Name Description
	x	State vector
	f (x)

Table 2 :

 2 Nomenclature of the two-stage gear system

	Name Description	Unit
	θ i	Angular position of the ith gear	[rad]
	R i	Bases radius of the ith gear	[m]
	I i	Inertia moment of the ith gear	[kg × m 2 ]
	I s	Inertia moment of the shaft	[kg × m 2 ]
	k e (t)	Time dependent gear-mesh stiness [N/m]
	k s	Torsional stiness of the shaft	[N × m/rad]
	q	Vector of degrees of freedom	[rad]
	F	External torques vector	[N × m]

Table 3 :

 3 Gear parameters

	The rst stage parameters	Values
	Module	2.5 mm
	Tooth number of the pinion Z 1	18
	Tooth number of the wheel Z 2	26
	Face width	20mm
	Pressure angle	20 deg
	Helix angle β (spur gear)	0 deg
	Base radius of the pinion R 1	21 mm
	Base radius of the wheel R 2	30 mm
	Rotational speed N 1	1674 rpm
	Rotational speed N 2	1159 rpm
	The second stage parameters	Values
	Module	2 mm
	Tooth number of the pinion Z 3	20
	Tooth number of the wheel Z 4	35
	Face width	38mm
	Pressure angle	20 deg
	Helix angle β (spur gear)	0 deg
	Base radius of the pinion R 3	18 mm
	Base radius of the wheel R 4	32 mm
	Rotational speed N 3	1159 rpm
	Rotational speed N 4