Ahmed Jhinaoui
email: ahmed.jhinaoui@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Joseph Morlier
email: joseph.morlier@isae.fr

Subspace Instability Monitoring for Linear Periodically Time-Varying Systems

Keywords: Helicopter dynamics, Time-varying systems, Periodic motion, Resonance, Subspace methods, Statistical inference, Fault detection, Stability limits

Most subspace-based methods enabling instability monitoring are restricted to the linear time-invariant (LTI) systems. In this paper, a new subspace method of instability monitoring is proposed for the linear periodically time-varying (LPTV) case. For some LPTV systems, the system transition matrices may depend on some parameter and are also periodic in time. A certain range of values for the parameter leads to an unstable transition matrix. Early warning should be given when the system gets close to that region, taking into account the time variation of the system. Using the theory of Floquet, some symptom parameter of stability-or residual-is defined. Then, the parameter variation is tracked by performing a set of parallel cumulative sum (CUSUM) tests. Finally, the method is tested on a simulated model of a helicopter with hinged blades, for monitoring the ground resonance phenomenon.

INTRODUCTION

Instability monitoring is currently the subject of extensive research activities motivated by an increasing requirement on the design of highly reliable control systems, for numerous applications such as civil engineering and aeronautics. Under some assumptions, such structures can be modeled by a linear time-invariant (LTI) model. The goal is to detect any anomalous behavior on these structures, and to alarm the user before any dysfunction or destruction occurs. There is a large amount of literature on the subject of instability monitoring, which has been handled with different approaches (see [START_REF] Angeli | On-line fault detection techniques for technical systems: A survey[END_REF] for details). One particular method, namely the subspacebased method, consists in comparing the characteristics of a system at a reference state with a subspace matrix given by new data corresponding to an unknown, possibly different state, as explained in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]. To check whether a change has occurred or not, a quasi-distance between the two states -called residual-is defined. Then, a statistical test decides if this residual is significantly different from zero. If some threshold is exceeded, the system has changed with respect to the reference. An overview of the different approaches of this general approach can be found in [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]. This method has been largely investigated and successfully tested on LTI systems. In some complex systems, the characterizing model may change due to some internal commands or parameters (Mach number for aircrafts, blade rotational velocity for helicopters...). In [START_REF] Basseville | In-flight monitoring of aeronautic structures: vibrationbased on-line automated identification versus detection[END_REF], the system is labeled as LTI, even if it depends on those internal commands (this class of systems is called linear parameter-varying (LPV) in control theory, but it is still considered time-invariant for each fixed parameter). This assumption holds since for each command's value, data has a time-invariant behavior, and since the command varies slowly as a step function, generally. The problem will be to detect a significant change of the structural characteristics leading to instability as those commands fluctuate and some system eigenvalue starts to cross the Nyquist circle. This problem has been solved for aircraft flutter monitoring in [START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF]. The subject of this paper is the extension of this approach to periodically time-varying (varying, even for fixed internal command's values) systems such as wind turbines and rotor systems, replacing the residual in [START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF] by a robust one considered in [START_REF] Döhler | Robust subspace based fault detection[END_REF]. To the best of the knowledge of the authors, the bibliography on subspace methods for LPTV systems is not abundant. For instance, in [START_REF] Verhaegen | A class of subspace model identification algorithms to identify periodically and arbitrarily time-varying systems[END_REF], some attempts have been carried out in subspace-based identification. Some other works have dealt with recursive tracking of periodic subspaces, as in [START_REF] Buzzi | Iterative cyclic subspace tracking for blind adaptive multiuser detection in multirate cdma system[END_REF]. Otherwise, Bayesian filtering such as Kalman and particle filters, or robust and adaptive observers are also used. Two key features motivate the study reported herein:

• A simple, necessary and sufficient criterion of stability can be extracted from periodic data, using the theory of Floquet [START_REF] Dacunha | A unified floquet theory for discrete, continuous, and hybrid periodic linear systems[END_REF]). Therefore, a simple residual based on this criterion can be built • A set of time-invariant data subsequences that have the same time-varying behavior can be found for such periodic systems, as outlined in [START_REF] Meyer | A unified analysis of multirate and periodically time-varying digital filters[END_REF], which makes it possible to apply a detection algorithm similar to the classical time-invariant one [START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF]), for each of these sequences

Once the criterion of stability is found, a set of residuals, corresponding to each time-invariant subsequence, is defined. Then, parallel statistical tests are performed to decide if an eventual change in the system characteristics has occurred. The statistical test used herein is the cumulative sum (CUSUM) because it is adequate for online implementation and for the detection of the instant of change (see [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]). The paper is organized as follows: in Section 2, the essential elements of the subspace methods for identification and instability monitoring are recalled. Then, Section 3 is devoted to the design of the new LPTV-extended method of instability monitoring. Finally, in Section 4 an illustrative simulation case is studied in order to test the efficiency of the suggested method.

SUBSPACE-BASED METHODS FOR LTI SYSTEMS

Subspace Identification

Consider the discrete time model, describing a system in the state space form:

z k+1 = F z k + w k y k = Hz k + v k (1)
with the state vector z ∈ R n , the output vector y ∈ R r , the state transition matrix F ∈ R n×n and the observation matrix H ∈ R r×n . The vectors w and v are unmeasured noises, assumed to be white Gaussian. The eigenstructure of (1) is given by the roots (λ, φ λ ) of the equations, below:

det(F -λI) = 0, F φ λ = λφ λ
(2) Let p and q be chosen parameters such that min{pr, qr} ≥ n. From the output data (y k ), a Hankel matrix H p,q ∈ R (p+1)r×qr is built:

H p,q =     R 1 R 2 • • • R q R 2 R 3 • • • R q+1 . . . . . . . . . . . . R p+1 R p+2 • • • R p+q     (3)
where the covariances of the output data write

R i = E(y k y T k-i
) and E is the expectation operator. When the number of data N goes to infinity, the approxima-

tions 1 below hold: Ri = 1 N N k=i+1 y k y T k-i and Ĥp,q = 1 N N -p k=q Y k + Y k -T , where Y k + = (y T k • • • y T k+p ) T and Y k -= (y T k-1 • • • y T k-q )
T . The Hankel matrix can be factorized as below (see [START_REF] Benveniste | Single sample modal identification of a nonstationary stochastic process[END_REF]):

H p,q = O p C q (4)
where

O p =     H HF . . . HF p     , C q = F G, F 2 G, • • • , F q G and G = E(z k y T k )
. Therefore, the observability matrix O p is ob-

1
In the following, we drop theˆnotation for sake of simplicity tained from a thin singular value decomposition (SVD) of the matrix H p,q and its truncation at the desired order n [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]):

H p,q = U ∆V T = [ U 1 U 2 ] ∆ 1 0 0 ∆ 2 V 1 T V 2 T (5) O p = U 1 ∆ 1/2 1 (6)
with U 1 the left subspace of H p,q found in the n first columns of U , and ∆ 1 the n upper singular values of ∆.

The observation matrix H is found in the first block-row of O p . To get the transition matrix F , the shift invariance property of O p is used in resolving the least square solution of [START_REF] Basseville | In-flight monitoring of aeronautic structures: vibrationbased on-line automated identification versus detection[END_REF]):

O ↑ p F = O ↓ p (7)
where:

O ↑ p =     H HF . . . HF p-1     , O ↓ p =     HF HF 2 . . . HF p     (8)
Once F is computed, the eigenvalues λ and the eigenvectors φ λ can easily be derived from the resolution of (2). Then, the mode shapes (the observed eigenvectors), defined by Φ λ = Hφ λ , are deduced. Let the eigenstructure be the vector θ = Λ vec(Φ) considered as a canonical parametrization of ( 1), where Λ is the vector whose elements are the eigenvalues and Φ is the matrix of the observed eigenvectors Φ λ = Hφ λ .

Subspace Instability Monitoring

In this section, the statistical subspace-based instability monitoring from [START_REF] Basseville | In-flight monitoring of aeronautic structures: vibrationbased on-line automated identification versus detection[END_REF] and [START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF] is recalled, where a residual function for change detection is associated with the subspace identification described in the previous section. This function, or residual, compares the system parameter θ 0 of a reference with a block Hankel matrix computed on new data (y k ) k=1•••N corresponding to an unknown, possibly different parameter θ. To decide whether a significant change has occurred or not, a local approach is used to decide between two hypotheses:

H 0 : θ = θ 0 and H 1 : θ = θ 0 + δθ/ √ N ( 9
)
where δθ is unknown but fixed. Let K be a left kernel of the Hankel matrix at the reference H 0 p,q , such that:

K T (θ 0 )H 0 p,q = 0, K T (θ 0 )K(θ 0 ) = I (10)
To check whether new data agree with the reference state corresponding to θ 0 , the residual function below is introduced in [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]:

ζ N (θ 0 ) = √ N vec(K T (θ 0 )H p,q )
, where H p,q is the Hankel matrix built with new data from the current parameter θ, possibly different from θ 0 . From now on, we will use another residual that is robust to noise changes, recently suggested in [START_REF] Döhler | Robust subspace based fault detection[END_REF]:

ζ N (θ 0 ) = √ N vec(K T (θ 0 )U 1 ) (11) where K is the left kernel of the left subspace U 0 1 of H 0 p,q
at the reference, and U 1 the left subspace of Hankel matrix H p,q built on the current data. Let J (θ 0 ) be the matrix of partial derivatives of the vectorized residual with respect to the vector of parameters [START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF]) and Σ(θ 0 ) is the residual covariance matrix, which can be estimated as an empirical sum [START_REF] Basseville | Subspace-based fault detection algorithms for vibration monitoring[END_REF]). The matrices Σ(θ 0 ) and J (θ 0 ) associated to the residual in (11) are independent of the possibly time varying noise properties [START_REF] Döhler | Robust subspace based fault detection[END_REF]). If Σ(θ 0 ) is positive definite and assuming θ as in (9), the central limit theorem (CLT) insures that the residual ζ N (θ 0 ) in ( 11) is asymptotically Gaussian distributed [START_REF] Döhler | Robust subspace based fault detection[END_REF]):

ζ N (θ 0 ) -→ N →∞ N (J (θ 0 )δθ, Σ(θ 0 )) (12)
which implies that the mean of ζ N (θ 0 ) converges to 0 only under H 0 (i.e. δθ = 0). In other words, detecting if the eigenstructure has changed from the reference is equivalent to detecting a change in the mean of ζ N (θ 0 ). In order to detect the direction of the change, a normal-

ized residual ζ N (θ 0 ) is introduced, such that: ζ N (θ 0 ) = J T (θ 0 )Σ -1 (θ 0 )ζ N (θ 0 )
, which leads to the normal distribution:

ζ N (θ 0 ) -→ N →∞ N (Σ(θ 0 )δθ, Σ(θ 0 )) (13)
where: Σ(θ 0 ) = J T (θ 0 )Σ -1 (θ 0 )J (θ 0 ). This covariance matrix is positive definite. Therefore, a positive (resp. negative) change in the mean of the normalized residual ζ N (θ 0 ) reflects a positive (resp. negative) change in δθ. A recursive decision function, or residual, is introduced in order to design an on-line monitoring algorithm:

Z k = J (θ 0 ) T Σ(θ 0 ) -1 vec K T (θ 0 )U 1,k (14)
where U 1,k is the current left subspace of the Hankel matrix at the k-th sample data y k , obtained by an update of the SVD. Based on the recursive residual in (14) and considering that the Z k 's are independent Gaussian which describe themselves a change in θ by a change in their mean, a commonly used CUSUM test [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]) is applied. In general, only the part Λ of θ changes with respect to time. The eigenvectors are assumed to be constant, considering the small displacement hypothesis (SDH). Now, denote ρ as the norm of one of the eigenvalues λ's of F and notice that when one of the ρ's increases to 1, the system goes toward instability [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]). The CUSUM test described below allows to detect an increase in any component of the parameter by using as Jacobian matrix the corresponding submatrix of J(θ 0 ) : the sensitivity J(θ 0 ) is replaced in ( 14) by the column J(ρ 0 ) corresponding to the sensitivity with respect to ρ at ρ equal to ρ 0 :

S k (ρ 0 ) = Σ(ρ 0 ) -1/2 k j=q Z j (ρ 0 ), T k (ρ 0 ) = min q≤j≤k S j (ρ 0 ) and g k (ρ 0 ) = S k (ρ 0 ) -T k (ρ 0 ).
The two hypotheses H 0,k and H 1,k to decide between, at each instant k, are now:

H 0,k : Eg k (ρ 0 ) ≈ 0 H 1,k : Eg k (ρ 0 ) > ( 15
)
where is some empirically tuned threshold [START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF]) that allows to detect only the significant changes. If this threshold is exceeded, the test reacts and the current data are describing a new value for the parameter ρ > ρ 0 . The test described above is applied to all the norms of the eigenvalues (λ i ) i=1•••n , and when a change is detected in one of the ρ i , the system is considered to be unstable.

SUBSPACE INSTABILITY MONITORING IN LPTV CASE

Linear periodically time-varying systems are encountered in many different fields. Some of the most relevant applications that involve periodic behaviors are rotor systems such as helicopters and wind turbines, satellite control of attitude or also in communications. This class of systems is considered as an intermediate class bridging the timeinvariant case to the time-varying one. For the latter class, the stability analysis is generally based on complex algebraic-analytic tools (see [START_REF] Bourles | Linear Time-Varying Systems: Algebraic-Analytic Approach[END_REF] for details). Finding a criterion of stability is a crucial step in subspace instability monitoring, because a residual based on this criterion has to be built and then tracked by the statistical test presented in the last section. It seems difficult to get a simple criterion for time-varying systems in general. Hopefully, there exist specific methods for periodic systems that give relatively simple criteria about stability. One of these methods is the modal analysis using the transformation of Floquet, described hereafter.

Transformation of Floquet

Let consider the continuous LPTV system below, with a period T :

ẋ(t) = A(t)x(t) + w(t), A(t + T ) = A(t), ∀t ∈ R y(t) = Cx(t) + v(t) (16)
with x ∈ R n the state vector, A ∈ R n×n the state matrix, y ∈ R r the output vector and C ∈ R r×n the observation matrix. The vectors w and v are two additive noises assumed to be white Gaussian. The idea of the Floquet transformation is to replace this system by an equivalent autonomous (i.e. the matrix of transition is no more function of time) one. Theorem 1. If A is continuous in time-or at least piecewise continuous-and an initial condition x(t 0 ) = x 0 is fixed, then a solution of ẋ(t) = A(t)x(t) is guaranteed to exist. Let Φ(t) be the matrix whose n columns are n linearly independent solutions, Φ(t) is known as the Fundamental Transition Matrix (FTM) and Φ(t) = A(t)Φ(t), Φ(t + T ) = Φ(t)Φ(T ).

Proof. see [START_REF] Dacunha | A unified floquet theory for discrete, continuous, and hybrid periodic linear systems[END_REF] Corollary 2. (Floquet Transformation) The value of the fundamental matrix at t = T is called the Monodromy Matrix.

Q = Φ(T ) (17)
Let R = 1 T log(Q) and x(t) = Φ(t)e -Rt z(t). Then, the first-order ordinary differential equation above can be transformed into the equivalent equation:

ż(t) = Rz(t) (18)
The equation of observation for the new variable z is y(t) = C (t)z(t) where C (t) = CΦ(t)e -Rt and C (t + T ) = CΦ(t + T )e -R(t+T ) = CΦ(t)Φ(T )e -RT e -Rt = CΦ(t)e -Rt = C (t). Therefore, any system that writes as in ( 16) can be transformed into an equivalent autonomous system with a periodic matrix of observation:

ż(t) = Rz(t) + Φ(t)e -Rt -1 w(t) y(t) = C (t)z(t) + v(t), C (t + T ) = C (t), ∀t (19)
According to Floquet theory, the continuous system ( 19) is stable if and only if the real parts of eigenvalues of R are negative or, similarly, if the norms of the eigenvalues of Q in ( 17) are inferior to one. Let consider the discretized form of ( 19), at a sampling rate 1 τ (where T is assumed to be a multiple of τ ) [START_REF] Ma | Quantifying robustness of biomedical network models[END_REF]):

z k+1 = F z k + Γ k w k y k = H k z k + v k (20)
The discretized monodromy matrix F = e Rτ , the discretized matrix of observation H k = C (kτ ) and Γ k = (k+1)τ kτ e Rγ dγ are periodic. The discrete period is T d = T τ .

Residuals Generation

As outlined in [START_REF] Meyer | A unified analysis of multirate and periodically time-varying digital filters[END_REF], the data (y k ) have different varying dynamics in the periodic case. However, the subsequences (y k0+iT d ) i∈N and (z k0+iT d ) i∈N are shown to have time-invariant behaviors for any k 0 ∈ N. A total of T d different time-invariant subsequences exists. Let the eigenstructure below be the parameter function of one of these subsequences (denoted the j-th subsequence, with j ∈ [[1, T d ]]):

θ (j) =       Λ vec(Φ j ) vec(Φ j+1 ) . . . vec(Φ j+p )       ( 21
)
where Λ is the vector of the eigenvalues λ i of F and the mode shapes Φ are the observed eigenvectors of F ,

Φ j = [H j φ 1 , • • • , H j φ n ] with F φ i = λ i φ i .
The condition in the discrete form is that the eigenvalues of F are in the Nyquist circle. The Hankel matrices are built using the time-invariant subsequences as described hereafter. Lemma 3. For any j ≥ q and if (N T d + p) data samples are available, the j-th Hankel matrix H (j) p,q is defined as:

H (j) p,q = 1 N N -1 i=0 Y + j+iT d Y -T j+iT d (22)
This matrix can be factorized as in (4) to the product:

H (j) p,q = O (j) p C (j) q (23)
where

O (j)
p is the observability matrix for the j-th subsequence, such that:

O (j) p = H T j (H j+1 F ) T • • • (H j+p F p ) T T (24)
Notice that unlike the LTI case, the observation matrix H j+i is varying with i = 1, (where is the floor operator), in order to get the sum for approximately the same subsequence.

The observability matrix can be obtained, similarly to Section 2, from a thin singular value decomposition (SVD) of the Hankel matrix H (j)

p,q :

H (j) p,q = U (j) 1 U (j) 2 ∆ (j) 1 0 0 ∆ (j) 2 (V (j) 1 ) T (V (j) 2 ) T (25) O (j) p = U (j) 1 (∆ (j) 1 ) 1/2 (26) Writing O (j)
p in the modal basis gives:

O (j) p = Φ T j (Φ j+1 D) T • • • (Φ j+p D p ) T T . T = O (j) p . T (27)
where diagonal matrix D is defined as D = diag(Λ), and T is an invertible transformation. When p is large enough (if the rank of the observability matrix is n), the j-th subsequence of the periodic system ( 20) is described by its observability matrix O Then, a left kernel K (j) of the left subspace U (j) 1 of H (j) p,q is deduced:

(K (j) ) T U (j) 1 = 0, ∀j ∈ [[1, T d ]] (28)
The T d equations in (28) are satisfied if current output data describe well the reference parameters θ (j) 0 . These functions are the candidates to be the residuals that compare current data to reference parameters. In the periodic case, there are T d residuals, one for each of the T d time-invariant subsequences. At an instant k, these residuals write in recursive forms ∀j ∈ [[1, T d ]], as in ( 14):

Z (j) k = (J (j) (θ (j) 0 )) T (Σ (j) (θ (j) 0 )) -1 vec(K (j) T U (j) 1,k ) (29) where K (j) is the left kernel at the reference θ (j) 0 , U (j)
1,k the left subspace update of the Hankel matrix H (j) p,q , when the k-th data sample of the (j)-th subsequence is available, and J (j) (θ (j) 0 ) and Σ (j) (θ (j) 0 ) are the residual's sensitivity and covariance matrices as in [START_REF] Döhler | Robust subspace based fault detection[END_REF].

Parallel CUSUM Tests

Let ρ 0 be the norm of one of the eigenvalues contained in Λ. The goal herein is to decide whether the value of ρ 0 has increased. For each time-invariant subsequence, the hypotheses to decide between at the time sample k are:

H (j) 0,k : Eg (j) k (ρ 0 ) ≈ 0 H (j) 1,k : Eg (j) k (ρ 0 ) > (30)
where is some predefined threshold that allows to detect only the significant changes. The T d CUSUM tests to perform, corresponding to each time-invariant subsequence, are:

S (j) k (ρ 0 ) = Σ (j) (ρ 0 ) -1/2 k i=q Z (j) i (ρ 0 ) T (j) k (ρ 0 ) = min q≤i≤k S (j) i (ρ 0 ) g (j) k (ρ 0 ) = S (j) k (ρ 0 ) -T (j) k (ρ 0 ) (31) Z (j)
i (ρ 0 ) is the recursive residual with respect to ρ 0 and replacing, in (29), the sensitivity J (j) (θ (j) 0 ) by the sensitivity J (j) (ρ 0 ). Notice that when the monitoring starts, no information is available about the relation in time between the identified reference and the current data. In other words, if the monitoring starts at t = t 0 (with t 0 ≥ q) and (N T d +p) data are available, it is not known to which kernel K (j) the Hankel matrix H (t0)

p,q = 1 N N -1 i=0 Y + t0+iT d Y -T
t0+iT d corresponds, because the phase of identification and the phase of monitoring may have not the same time origin. The current data is synchronized with the reference parameter kernel (denoted K (j * ) ) as follows:

K (j * ) := argmin j=1•••T d ( (K (j) ) T H (t0) p,q ) (32)
To summarize the algorithm:

• using the identification algorithm described in [START_REF] Jhinaoui | Subspace identification for linear periodically time-varying systems[END_REF], a set of T d kernels K (j) , corresponding to the T d time-invariant subsequences is computed at a fixed reference θ 0 = (θ (j) 0 ) j=1,••• ,T d • the Jacobians and the covariances matrices are computed (see [START_REF] Döhler | Robust subspace based fault detection[END_REF]) • start the monitoring. If it starts at t = t 0 , the current data are synchronized with the identified kernels using (32) • for t ≥ t 0 , the CUSUM tests are performed for each subsequence, in parallel, as in (31) • the instant time of change is defined as the instant when one (or more) of the parallel CUSUM tests exceeds some fixed threshold [START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF])

Notice also, that for some complex systems, the mechanical model is function of some internal parameters, such as rotational velocity Ω for helicopters' blades. In this case, the system is periodic of period T = 2π Ω . And then, the period of the discrete system is T d = 2π

Ωτ . When Ω = Ω 1 , the system is LPTV and there exist T d1 = 2π

Ω1τ timeinvariant subsequences, and when Ω = Ω 2 > Ω 1 , the system is LPTV and there exist T d2 = 2π

Ω2τ time-invariant subsequences, with T d2 < T d1 . As in the case reported in [START_REF] Basseville | In-flight monitoring of aeronautic structures: vibrationbased on-line automated identification versus detection[END_REF] and assuming that the velocity is slowly varying as step function, for two different velocities Ω 1 and Ω 2 , the T d2 subsequences (y 1+iT d 2 ), (y 2+iT d 2 ), • • • , (y T d 2 +iT d 2 ) from Ω = Ω 2 and (y 1+iT d 1 ), (y 2+iT d 1 ), • • • , (y T d 2 +iT d 1 ) from Ω = Ω 1 are considered as coming from the same subsystem, whereas, the T d1 -T d2 other subsequences vanish when the velocity changes to Ω = Ω 2 .

ILLUSTRATIVE EXAMPLE

Helicopter Model

The proposed method is applied to simulation data of a helicopter on the ground. The goal is to detect the ground resonance before it occurs. Similar to that proposed in [START_REF] Byers | Embedded absorbers for helicopter rotor lag damping[END_REF], the present mechanical model is developed to characterize the dynamic behavior of a helicopter with a hinged rotor (see Fig. 1). The fuselage is considered to be a rigid body with mass M .

M(t) Ẍ(t) + C(t) Ẋ(t) + K(t)X(t) = 0 (33)
The mass, the damping and the stiffness matrices M, C and K are periodic of period T = 2π Ω and X =

[ x y φ 1 φ 2 • • • φ N b ]
T . Let X = X T ẊT T , then (33) can be written in a state space form:

Ẋ (t) = A(t)X (t) (34)
Where the periodic state transition matrix writes:

A(t) = 0 I -M -1 (t)K(t) -M -1 (t)C(t) (35)
The equation of observation is y(t) = CX (t) where C = [ 0 I ]. The numerical values used for simulation are those reported in [START_REF] Byers | Embedded absorbers for helicopter rotor lag damping[END_REF] and N b = 4.

Simulation

Using Matlab, time series data are simulated from the mechanical model above with a sampling rate 1 τ = 50Hz. The scenario consists in simulating an angular velocity's acceleration from Ω 0 = 3.1rad/s to Ω f inal = 4.5rad/s (close to resonance) with a step of 0.1rad/s. For each value of Ω, 1000 samples are simulated and the eigenvalues of F in (20) are computed. The imaginary parts of these eigenvalues are negligible. Therefore, the real parts ρ 1 • • • ρ n , which are positive, are close approximations of the norms of these eigenvalues. The evolution of the ρ i 's with respect to airspeeds Ω is plotted in Fig. 2. The resonance onset can be observed when one of the ρ i 's exceeds 1. Notice that the ρ i 's vary very close to instability.

Results

At the chosen reference corresponding to Ω 0 = 3.1rad/s, the system is identified on a large data set from the subspace identification algorithm described in [START_REF] Jhinaoui | Subspace identification for linear periodically time-varying systems[END_REF]. The kernels, Jacobians and covariance matrices of the T d0 = 102 time-invariant data subsequences are computed at this reference [START_REF] Döhler | Robust subspace based fault detection[END_REF]). At the beginning of the monitoring phase, the Hankel matrix computed on current data is synchronized with respect to its proper reference kernel (as given in ( 32)). This synchronization is performed using (32) whose minimum is plotted in Fig. 3. The x-axis represents the index of the reference kernels, and the y-axis the value of the cost in equation ( 32). It corresponds to the correct phase (as defined by the simulation) at an index equal to 50. Once The problem of instability monitoring for linear periodically time-varying systems is addressed. The most challenging point is to define a simple criterion of stability for these LPTV systems. This problem is solved using the Floquet theory. Then, new instability monitoring tests are developed, similarly to the time-invariant case. The suggested method is successfully tested on simulation data of a helicopter model. Current works are carried out to apply the same method to real industrial applications.

  ,••• ,T d , the observability matrices should be given for all the T d subsequences (y 1+iT d ) i , (y 2+iT d ) i , (y 3+iT d ) i , • • • , (y T d +iT d ) i . For each subsequence (j) j=1•••T d , the Hankel matrix is computed as in (22).

  Fig. 1. Helicopter's mechanical model system writes:

Fig. 2 .

 2 Fig. 2. Real parts ρ 1 to ρ 6 vs. angular velocity (N.B: ρ 1 ρ 2 )

Fig. 3 .

 3 Fig. 3. Cost function vs. indexes of reference kernels the synchronization is done, the monitoring phase begins. For each of the T d time-invariant subsequences, a CUSUM test is started. The number of tests T d is decreasing with respect to the angular velocity Ω which means that some time-invariant subsequences have vanished. At Ω f inal , only 70 subsequences (then, 70 parallel tests) will remain. A choice of = 1000 was taken based on the fluctuation of the tests in the reference zone for Ω[START_REF] Mevel | Fast in-flight detection of flutter onset -a statistical approach[END_REF]). If one of the tests reacts (if some ρ i 's have significantly increased. Here it is ρ 5 that increases), monitoring is stopped. Two of these test evolutions are reported in Fig.4. All tests have the same behavior and their reactions to the increase in parameter value are clearly observed. They all react and the earliest response corresponds to Ω = 4.1rad/s, which is close to the zone of resonance.

  Proof. see[START_REF] Jhinaoui | Subspace identification for linear periodically time-varying systems[END_REF] Fact 4. For notation simplicity, T is assumed to be a multiple of τ . In general, this is not the case and the timeinvariant subsequences are the ensembles (y jτ +iT τ ) i and the data Y + j+iT d and Y - j+iT d are replaced by Y +

		jτ +iT	and
	Y -jτ +iT	τ
	τ

• • • , p.

This work was supported by the European project FP7-NMP CP-IP 213968-2 IRIS. The authors acknowledge Leonardo Sanches for insightful discussions and comments.


Ahmed Jhinaoui
email: ahmed.jhinaoui@inria.fr

Laurent Mevel
email: laurent.mevel@inria.fr

Joseph Morlier
email: joseph.morlier@isae.fr

Extension of Subspace Identification to LPTV Systems: Application to Helicopters

In this paper, we focus on extending the subspace identification to the class of linear periodically time-varying (LPTV) systems. The Lyapunov-Floquet transformation is first applied to the system's state-space model in order to get the monodromy matrix (MM) and, thus, a necessary and sufficient condition for system stability. Then, given two successive covariance-driven Hankel matrices, the MM matrix is extracted by some calculus of a simultaneous singular value decomposition (SVD) and a least square optimization. The method is illustrated by a simulation application to the model of a hinged-blades helicopter.

INTRODUCTION

Over the last forty decades, subspace identification methods have enjoyed of some popularity and numerous applications of these methods have emerged in civil engineering, aeronautics and many other fields. The stochastic subspace identification (SSI) consists in obtaining the modal parameters (natural frequencies, modal damping ratios and mode shapes) of a system subject to vibrations, by some geometric manipulations and projections of data given by sensors' measurements or input measurements. There exist, then, two types of identification algorithms: output-only and input-output algorithms. A comprehensive overview of different approaches of SSI can be found in [START_REF] Ljung | System Identification: Theory for the User[END_REF] Unfortunately, the most part of research interest has been given to linear time-invariant (LTI) systems. In contrast, the literature on time-varying (TV) case is not abundant. The attempts to extend these methods to TV can be categorized into two main classes, as outlined in [START_REF] Liu | Identification of linear time-varying systems[END_REF]. The first class consists in identifying the considered system recursively using adaptive algorithms which is appropriate only for slowlyvarying dynamics and when a priori information about the variation behavior is available . Whereas, the second class suggests to find a set of output sequences that have the same time-varying behavior [START_REF] Verhaegen | A class of subspace model identification algorithms to identify periodically and arbitrarily time-varying systems[END_REF], which makes it possible to apply a classical time-invariant identification algorithm to these sequences. We suggest a new algorithm of identification in this paper, with three main features: first the a covariance-based method is proposed, which makes the computation non cumbersome. Second, we combine the method proposed to the theory of Floquet that gives a simple criterion for the stability of a periodic system. With the Lyapunov-Floquet transformation, the state transition matrix is replaced by the monodromy matrix. The third feature is to find out this monodromy matrix without any approximation (not up to a matrix of similarity). The paper is organized as follows: in Section 2, a typical SSI algorithm is presented. Section 3 gives the essential elements of the Floquet theory. Then, Section 4 is devoted to the design of the new LPTV-extended method. Finally, in Section 5 the efficiency of the method is tested on a numeric simulation of a model of a helicopter with a hinged-blades rotor.

Covariance-driven Subspace Identification

In this section, we present a typical SSI algorithm [START_REF] Ljung | System Identification: Theory for the User[END_REF]. Let consider the discrete-time state space model of a given system:

x k+1 = F x k + v k+1 y k = H x k + w k (1)
where x ∈ R n is the state vector, y ∈ R r the output vector or the observation, F ∈ R n×n the state transition matrix and H ∈ R r×n the observation matrix. The vectors v and w are noises that are assumed to be white Gaussian and zero-mean. We choose the parameters p and q, and we build the covariance-driven Hankel matrix:

H p+1,q =      R 0 R 1 • • • R q-1 R 1 R 2 • • • R q . . . . . . . . . . . . R p R p+1 • • • R p+q-1     
The covariances of the output data write R i = E(y k y T k-i ), where E is the expectation operator. Remark 1. : R i can be approximated by Ri = 1 N ∑ N k=i+1 y k y T k-i if N is the number of output measurements we have (N 1). Then, the estimated Hankel matrix:

Ĥp+1,q = Hank( Ri ) = 1 N ∑ N k=q Y k + Y k -T
Where:

Y k + = (y k • • • y k+p ) T , Y k -= (y k • • • y k-q+1 ) T
Let G = E(x k y k ) be the correlation between the state and the observation.

O p+1 =      H HF . . . HF p     
and C q = G FG • • • FG q-1 are the (p + 1)th order observability matrix and the qth order controllability matrix. The computation of the R i 's leads to the decomposition:

H p+1,q = O p+1 C q
Therefore, the observability matrix O p+1 can be obtained with a thin SVD of the Hankel matrix H p+1,q and its truncation at the desired model order n:

H p+1,q = U 1 U 2 Σ 1 0 0 Σ 2 V 1 T V 2
T and:

O p+1 = U 1 Σ 1 1 2
Where Σ 1 contains the first n singular values.

The observation matrix H is extracted from the first r rows of the observability matrix O p+1 . The state transition matrix F is obtained from a least squares solution of:

O ↑ p+1 F = O ↓ p+1
Where:

O ↑ p+1 =      H HF . . . HF p-1      , O ↓ p+1 =      HF HF 2 . . . HF p     
We get the eigenstructure of the system (1) from the resolution of det(F -Iλ ) = 0 and Fφ λ = λ φ λ , (λ , φ λ ) denote the eigenvalues and the eigenvectors of the system.

Floquet Theory

The theory of Floquet is a mathematical theory of ordinary differential equations (ODEs). Introduced in 1883, it is the first complete theory for the class of periodically time-varying systems. In this section, we briefly review some of its essential elements that are related to our study. More details can be found in [START_REF] Dacunha | A unified floquet theory for discrete, continuous, and hybrid periodic linear systems[END_REF].

Let consider the periodic differential system:

ẋ(t) = A(t)x(t) (2)
Where x ∈ R n is the state vector. The state transition matrix A ∈ R n×n is continuous in time (or at least, piecewise continuous) and periodic, of period T > 0. If an initial condition x(t 0 ) = x 0 is fixed, a solution of ( 2) is guaranteed to exist. Let Φ(t) be the matrix whose n columns are n linearly independent solutions of (2), Φ(t) is known as the Fondumental Transition Matrix (FTM). It has the properties below:

Φ(t) = A(t)Φ(t), Φ(t + T ) = Φ(t)Φ(T ), ∀t

Stability Analysis

The value of the fundemantal matrix at

t = T , Q = Φ(T ) is called the Monodromy Matrix. Let R = 1 T log(Q).
According to Floquet's theory, the dynamical system ( 2) is stable if and only if the eigenvalues of R are negative or, similarly, if the norm of the eigenvalues of Q is less than one.

Floquet Transformation

The Floquet transformation, also called the Lyapunov-Floquet transformation, gives an underlying autonomous system (a system with a constant state transition matrix) that is equivalent to the initial periodic system i.e: the transformation is invertible. If we make the change of variable x(t) = Φ(t)e -Rt z(t), the theory insures that:

ż(t) = Rz(t) (3)
Now if the equation of observation of the considered system is: y(t) = Cx(t), the equivalent equation for the new variable z is:

ẏ(t) = Cx(t) = CΦ(t)e -Rt z(t) (4)
Remark 2. : It is easy to demonstrate that Φ(t + T )e -R(t+T ) = Φ(t)e -Rt using the fact that R = 1 T log(Φ(T )). Therefore, any periodic system can be transformed into an equivalent autonomous system with an equivalent periodic matrix of observation.

We have presented above, a subspace method for time-invariant systems. In this section, we suggest an extension of this method to the linear periodically time-varying case. Let consider the periodic state-space system, of a period T :

ẋ(t) = A(t)x(t) y(t) = Cx(t) (5)
As outlined in Section 3, an equivalent representation can be given by the Floquet's transformation:

ż(t) = Rz(t) y(t) = CΦ(t)e -Rt z(t) (6)
Sampling at rate 1 τ (τ > 0) yields the discrete time model below:

z k+1 = F z k + v k+1 y k = H k z k + w k , z k = z(kτ) (7)
The discrete state transition and observation matrices are: F = e τR and H k = CΦ(kτ)e -Rkτ . The output and the state are corrupted by the noises v and w assumed to be white Gaussian and zero means. If µ is an eigenvalue of the continuous system R, the correspondent eigenvalue λ of F is such that: λ = e τ µ . According to the theory of Floquet, the continuous system is stable when the eigenvalues µ are negative. Therefore, the discrete system is stable when λ < 1 (or, if λ is complex, its real part is under 1). The discrete observation matrix H k is periodic of period T d = T τ + 1 ( denotes the floor operator). We build the instantaneous product

Y k + Y k -T : Y k + Y k -T =      H k H k+1 F . . . H k+p F p      [ z k y T k Fz k-1 y T k-1 • • • • • • Fz k-q+1 y T k-q+1 ]
Following the lines of [START_REF] Park | Canonical forms on discrete linear periodically time-varying systems and a control application[END_REF], the ensembles z k z k+T d z k+2T d z k+3T

G (k) = E(z k y T k ) = 1 N ∑ N-1 i=0 z k+iT d y T k+iT d .
A Hankel matrix is formed from these cross-products, at sample k:

H (k) p+1,q = ∑ N-1 i=0 Y k+iT d + Y k+iT d -T = ∑ N-1 i=0      H k+iT d H k+1+iT d F . . . H k+p+iT d F p      z k+iT d y T k+iT d • • • F q-1 z k-q+1+iT d y T k-q+1+iT d (8)
We know that the matrix of observation H is periodic. Then

H k+iT d = H k , ∀i > 0.
Using this property, the matrix of Hankel writes:

H (k) p+1,q =      H k H k+1 F . . . H k+p F p      ∑ N-1 i=0 z k+iT d y T k+iT d • • • F q-1 ∑ N-1 i=0 z k-q+1+iT d y T k-q+1+iT d =      H k H k+1 F . . . H k+p F p      G (k) • • • • • • F q-1 G (k-q+1)
As in the time-invariant case, the matrix of Hankel H

p+1,q can be decomposed in a product of two matrices:

H (k) p+1,q = O (k) p+1 C (k) q (9)
Where:

O (k) p+1 =      H k H k+1 F . . . H k+p F p      and C (k) q = G (k) • • • F q-1 G (k-q+1)
The observability matrix O k p+1 can be obtained as in Section 2 via an SVD of H (k) p+1,q and its truncation at the desired model order n. The observation matrix H k , at k, is obtained from the first r rows of

O (k) p+1 .
The extraction of the transition matrix from one observability matrix is no longer possible. In order to get F, we have to compute two successive matrices of Hankel. At (k + 1), H (k+1) p+1,q writes:

H (k+1) p+1,q =      H k+1 H k+2 F . . . H k+p+1 F p      G (k+1) • • • • • • F q-1 G (k-q+2)
Now if we have the two successive observability matrices O p+1 , the transition matrix F is -as in Section 2-the Moore-Penrose (least square) solution of:

O ↑ (k+1) p+1 F = O ↓ (k) p+1 (10)
Where:

O (k) p+1 ↓=      H k+1 F H k+2 F 2 . . . H k+p+1 F p      , O (k+1) p+1 ↑=      H k+1 H k+2 F . . . H k+p+1 F p-1      Remark 3. : Notice that the range space of O (k) p+1 and O (k+1)
p+1 must be in the same basis. In fact, the left part of the singular value decomposition gives the observability matrix up to some similarity matrix i.e for example: in the time-invariant case,

O p+1 = TU 1 Σ 1 1 2
where T is a unitary matrix (we have omitted this discussion until now, for simplicity). This matrix will be simplified by its inverse, with the Moore-Penrose resolution. And the obtained F corresponds, then, to the real transition matrix with no similarity.

In the case of time-variant systems, if the two successive SVD are done separately, we will have two different similarity matrices T k and T k+1 . When it comes to the least square resolution, no simplification will be possible:

F = (O ↑ (k+1) p+1 ) † T -1 k+1 T k O ↓ (k) p+1 (11)
The symbol † denotes the Moore-Penrose pseudo-inverse. None of the similarity matrices T k or T k+1 is known. A solution for this is to force an SVD in the same basis for the two successive Matrices of Hankel:

H (k) p+1,q H (k+1) p+1,q = U 1,k U 2,k U 1,k+1 U 2,k+1 Σ 1,k 0 0 Σ 2,k V T k,k+1 (12)
The size of covariance-driven matrices of Hankel is (p + 1)r × qr. Therefore, the complexity of the SVD does not depend on N, and the computation is not cumbersome.

The range spaces of the observability matrices, up to the same similarity matrix, are

O (k) p+1 = U 1,k Σ 1,k 1 2 , O (k+1)
p+1 = U 1,k+1 Σ 1,k 1 2 (13)
Once F and H k are computed, the modal structure of the system (7) are obtained from the resolution of det(F -Iλ ) = 0, Fφ λ = λ φ λ and ϕ λ = H k φ λ .

Algorithm: To sum up, here are the steps of the suggested LPTV identification method:

• NT d output-measurements are available. We set the data in ensembles of timeinvariant series y k y k+T d y k+2T d y k+3T d • • • for each sample k

• we compute, at k and (k + 1), the two successive matrices of Hankel H

p+1,q and H (k+1) p+1,q using the formula (8)

• we make a simultaneous singular value decomposition of those two matrices as in (12) • we compute the range spaces of the observability matrices, at k and (k + 1), as in (13) • we resolve the equation (10) using the Moore-Penrose pseudo-inverse

• H k is obtained as the first r rows of O (k) p+1
• given F and H k , we compute the eigenstructure of the system

Application

In this section, we give an illustrative example in order to test the new suggested identification algorithm. The system we study herein is a helicopter with a hingedblades rotor. The chosen modeling is the modeling for the analysis of the ground resonance phenomenon. We present the dynamical equations of motion that describe this phenomenon, and develop the linear periodically time-varying mechanical model in a similar way as described in [START_REF] Sanches | Instability zones identification in multiblade rotor dynamics[END_REF], but adding damping to the structure this time. The helicopter's fuselage is considered to be a rigid body with mass M, attached to a flexible LG (landing gear) which is modeled by two springs K x and K y , and two viscous dampers C x and C y as illustrated in Fig. 1. The rotor spinning with a velocity ω, is articulated and the offset between the MR (main rotor) and each articulation is noted a. The blades are modeled by a concentrated mass m at a distance b of the articulation point. Torque stiffness and a viscous damping K β and C β are present into each articulation. The moment of inertia around the articulation point is I z . The degrees of freedom are the lateral displacements of the fuselage x and y, and the out-of-phase angles β k=1•••N b , with N b the number of blades.

Let Z = x y β 1 • • • β N b
T . A linear model for the considered system under free vibrations is defined by:

M (t) Z (t) + C (t) Ż (t) + K (t)Z (t) = 0 (14)
Where:

M (t) =        (M + Nm) 0 -mbsin(ωt) 0 (M + Nm) mbcos(ωt) -mbsin(ωt) mbcos(ωt) (mb 2 + I z ) -mbsin(ωt + α) mbcos(ωt + α) 0 . . . . . . . . . -mbsin(ωt + α) • • • -mbsin(ωt + (N -1)α) mbcos(ωt + α) • • • mbcos(ωt + (N -1)α) (mb 2 + I z ) • • • 0 . . . . . . . . .      C (t) =        C x 0 -2mbωcos(ωt) 0 C y -2mbωsin(ωt) 0 0 C β 0 0 0 . . . . . . . . . -2mbωcos(ωt + α) • • • -2mbωcos(ωt + (N -1)α) -2mbωsin(ωt + α) • • • -2mbωsin(ωt + (N -1)α) 0 • • • 0 . . . . . . . . .      and
K (t) =        K x 0 mbω 2 sin(ωt) 0 K y -mbω 2 cos(ωt) 0 0 K β 0 0 0 0 . . . . . . . . . mbω 2 sin(ωt + α) • • • mbω 2 sin(ωt + (N -1)α) -mbω 2 cos(ωt + α) • • • -mbω 2 cos(ωt + (N -1)α) 0 • • • 0 . . . . . . . . .     
The system can be written as in [START_REF] Dacunha | A unified floquet theory for discrete, continuous, and hybrid periodic linear systems[END_REF], with x = Z T Ż T T and

A(t) = 0 I -M (t) -1 K (t) -M (t) -1 D(t) , A(t) = A(t + 2π ω ). It is periodic of pe- riod T = 2π ω .
The numerical values used for the application are reported in Fig. 2. The helicopter model is simulated for ω = 4rad.s -1 . To get precise results, the subspace identification should be processed on very large dataset (to overcome the problem of bias due to the noise i.e 1 N ∑ k v k ≈ 0 and 1 N ∑ k w k ≈ 0). A total of N = 80000 data points are generated, with a sampling frequency of 50Hz. The order of the system is known and is equal to n = 12. The identification algorithm is applied to the data as explained in Section 4. The summary of the identified eigenvalues, and the difference between them and the real values, are given in Fig. 3. The differences in the obtained real part (damping ratios) are less than 0.1 %. For the imaginary part (frequencies), the differences are larger, but still not significant considering the large uncertainty on the estimates. The identification is more accurate when N is larger.

Structural variable Value

m 31.9Kg M 2902.9Kg K β 200N/m K x 3200N/m K y 3200N/m C β 15Ns/m C x 300Ns/m C y 300Ns/m a 0.2m b 2.5m I z 259Kg/m 2
The most important variable the real part of the eigenvalues, because it is the criterion for deciding whether the system is stable or not. Therefore, the suggested method of identification is sufficiently accurate for the purpose of stability analysis.

Identified mode

Real mode damping frequency damping frequency 0.9974 0.0766 0.9998 0.0200 0.9974 -0.0766 0.9998 -0.0200 0.9989 0.0702 0.9998 0.0214 0.9989 -0.0702 0.9998 -0.0214 0.9995 0.0556 0.9994 0.0359 0.9995 -0.0556 0.9994 -0.0359 0.9985 0.0130 0.9995 0.0330 0.9985 -0.0130 0.9995 -0.0330 0.9971 0.0310 0.9994 0.0335 0.9971 -0.0310 0.9994 -0.0335 0.9999 0.0471 0.9994 0.0335 0.9999 -0.0471 0.9994 -0.0335 Fig. 3 Identified modes vs. real modes

Conclusion

The problem of identification for linear periodically time-varying is addressed. An extension of the covariance-driven SSI algorithm is proposed and tested with simulation data. Future works encompasses a generalization of the subspace detection methods to the same class of systems.

Acknowledgments

This work was supported by the European project FP7-NMP CPIP 213968-2 IRIS.

Fig. 1

 1 Fig. 1 Mechanical model of a helicopter with 3 blades

Fig. 2

 2 Fig. 2 Structural properties for hinged-blades helicopter with 4 blades

  • • • and y k y k+T d y k+2T d y k+3T d • • • are time-invariant series. If we have, for example, NT d output measurements (N 1), the output-input covariance of the k-th invariant ensemble can be estimated by