Qinghua Zhang
email: qinghua.zhang@inria.fr

Stochastic Hybrid System Actuator Fault Diagnosis by Adaptive Estimation

Keywords: Fault diagnosis, hybrid system, random uncertainty, actuator fault, joint state-parameter estimation

Based on the interacting multiple model (IMM) estimator for hybrid system state estimation and on the adaptive Kalman filter for time varying system joint state-parameter estimation, a new algorithm, the adaptive IMM estimator, is proposed in this paper for actuator fault diagnosis in stochastic hybrid systems. The working modes of the considered hybrid systems are described by stochastic state-space models, and the mode transitions are characterized by a Markov model. Actuator faults are modeled as parameter changes, and the related fault diagnosis problem is solved by the proposed adaptive IMM estimator through joint stateparameter estimation. Numerical examples are presented to illustrate the performance of the proposed method.

INTRODUCTION

In order to improve the safety, the reliability and the performance of complex industrial systems, the topic of fault diagnosis is attracting more and more researchers during the last decades [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]; [START_REF] Gertler | Fault detection and diagnosis in engineering systems[END_REF]; [START_REF] Chen | Robust Model-Based Fault Diagnosis for Dynamic Systems[END_REF]; [START_REF] Ding | Model-Based Fault Diagnosis Techniques -Design Schemes Algorithms and tools[END_REF]; [START_REF] Isermann | Fault-Diagnosis System[END_REF]; [START_REF] Korbicz | Fault Diagnosis: Models, Artificial Intelligence, Applications[END_REF]; [START_REF] Blanke | Diagnosis and Fault-Tolerant Control[END_REF]; [START_REF] Simani | Model-based Fault Diagnosis in Dynamic Systems Using Identification Techniques[END_REF]). Most of these studies are model-based, assuming that the considered systems are described by differential equations or by their discrete time counterpart. However, complex industrial systems may have behaviors that cannot be described by a single set of differential equations, because of working mode changes. Typically, each of the working modes is modeled by a different set of differential equations. Such systems are known as hybrid systems. It is thus necessary to develop methods for hybrid system fault diagnosis in order to cope with complex industrial systems.

In general a hybrid system involves both continuous dynamics and discrete events that may be combined in a more or less complex manner (Blom and Lygeros (2006); Lunze and Lamnabhi-Lagarrigue ( 2009)). In this paper it is assumed that each considered hybrid system has a finite number of working modes described by discrete time stochastic state-space models subject to parameter changes, and that at every time instant one of the working modes is active. If the active mode is known all the time, then the hybrid system can be treated as a (discontinuous) time varying system, and in this case some fault diagnosis methods designed for time varying systems can be applied [START_REF] Chung | A game theoretic fault detection filter[END_REF]; [START_REF] Chen | A generalized least-squares fault detection filter[END_REF]; [START_REF] Chen | Optimal stochastic fault detection filter[END_REF]; [START_REF] Zhong | Optimal fault detection for linear discrete time-varying systems[END_REF]; [START_REF] Zhang | Statistical detection and isolation of additive faults in linear time-varying systems[END_REF]). It is more challenging to study hybrid systems whose active mode is unknown. For hybrid systems with a This work has been supported by the ITEA2 MODRIO project. deterministic mode switching mechanism, fault diagnosis has been studied in [START_REF] Bemporad | Moving horizon estimation for hybrid systems and fault detection[END_REF]) through the mixed logic dynamic formalism (see also [START_REF] Ferrari-Trecate | Moving horizon estimation for hybrid systems[END_REF]), in [START_REF] Belkhiat | Design of robust fault detection observers for discretetime linear switched systems[END_REF]) with state observers, and in [START_REF] Wang | Model-based Health Monitoring of Hybrid Systems[END_REF]) with a bond graphbased approach. When the deterministic mode switching mechanism can be determined from observations within a finite time, a hybrid observer-based method has been proposed in [START_REF] Wang | Robust state estimation and fault diagnosis for uncertain hybrid nonlinear systems[END_REF]). Stochastically switching hybrid systems are considered in [START_REF] Cinquemani | Fault detection in a class of stochastic hybrid systems[END_REF]), with each monitored fault modeled as one of the modes of the hybrid system. More generally, particle filters can be applied to stochastic hybrid systems for fault diagnosis [START_REF] Guo | Fault diagnosis of hybrid systems using particle filter based hybrid estimation algorithm[END_REF]), but such solutions are numerically expensive.

For the hybrid systems considered in this paper, each mode is described by a state-space model subject to Gaussian noises, and the mode switching mechanism is characterized by a Markov model. Such systems exhibit stochastic behaviors both in each mode and during mode transitions. In this framework, actuator faults are modeled as parameter changes, typically the gain losses of actuators. The magnitude of each parameter change, when it occurs, is an unknown real value, possibly belonging to some bounded interval. It is thus impossible to model such faults as modes of a hybrid system, as each fault would correspond to a particular value of the parameter change magnitude and infinitely many modes would have to be considered.

In the fault-free case, the hybrid systems considered in this paper have been well studied for the problem of state estimation [START_REF] Blom | The interacting multiple model algorithm for systems with markovian switching coefficients[END_REF]; [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation[END_REF]; [START_REF] Hwang | State estimation for hybrid systems: applications to aircraft tracking[END_REF]), notably with the well-known Interacting Multiple Model (IMM) estimator [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation[END_REF]). Such results cannot be directly applied to the systems considered in this paper for state estimation, because of unknown parameter changes. On the other hand, for time varying state-space systems, efficient joint state-parameter estimation methods exist, for instance the adaptive Kalman filter, also known as adaptive observer for continuous time systems [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (MIMO) linear time varying systems[END_REF]; [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]). The main idea of the present paper is to combine the IMM estimator and the adaptive Kalman filter, in order to design an adaptive IMM estimator for joint state-parameter estimation of hybrid systems. Then the resulting algorithm can be directly applied to the actuator fault diagnosis problem considered in this paper.

PROBLEM STATEMENT

The stochastic hybrid system considered in this paper is modeled at two levels. At the top level, the system has a finite number of working modes. At each time instant, one of the modes is active, and random transitions between different modes are characterized by a Markov model. At the bottom level, each mode of the system is described by a stochastic linear state-space model subject to actuator faults formulated as parameter changes.

Markov transition model

Assume that a hybrid system has r working modes, labeled by M 1 , M 2 , . . . , M r . At the initial time instant k = 0, the prior probability that (the mode labeled by1 ) M j is active is

P {M j } = µ j (0) (1) with known probabilities µ 1 (0), µ 2 (0), . . . , µ r (0) satisfying r j=1 µ j (0) = 1.
The mode switching mechanism is characterized by a Markov process with the transition probabilities P {M j (k)|M i (k -1)} = p i,j .

(2) where k is the discrete time instant, and p i,j are known transition probabilities independent of k and satisfying r j=1 p i,j = 1 for i = 1, 2, . . . , r.

Stochastic state-space mode model

In each of the possible r modes, say M j , the considered hybrid system is described by the linear state-space model

x(k) = A(M j )x(k -1) + B(M j )u(k -1) + w(k) + Φ(k -1; M j )θ (3a) y(k) = C(M j )x(k) + v(k) (3b) where x(k) ∈ R n is the state, u(k) ∈ R l the in- put, y(k) ∈ R m the output, A(M j ), B(M j ), C(M j ) are mode-dependent matrices of appropriate sizes, w(k) ∈ R n , v(k) ∈ R m
are mutually independent white Gaussian noises of covariance matrices Q(M j ) ∈ R n×n and R(M j ) ∈ R m×m respectively, and the term Φ(k; M j )θ represents actuator faults with a known matrix sequence Φ(k; M j ) ∈ R n×p and a constant (or piecewise constant) vector θ ∈ R p . At the initial time instant k = 0, the initial state x(0), under the assumption of each possible mode M j , is assumed to be a Gaussian random vector

x(0) ∼ N (x j (0|0), P j (0|0)).

(4)

A typical example of actuator faults represented by Φ(k; M j )θ is in the case where θ corresponds to actuator gain loss coefficients. Assume that each of the l actuators, say the one corresponding to the q-th component of the vector u(k), is affected by a gain loss represented by a coefficient (1 -θ q ), where θ q is the q-th component of θ. Then the input term changes from its nominal form B(M j )u(k) to the faulty form

B(M j )(I l -diag(θ))u(k) = B(M j )u(k)-B(M j )diag(u(k))θ
where I l is the l × l identity matrix. In this case, p = l and Φ(k;

M j ) = -B(M j )diag(u(k)).
(5)

Actuator fault diagnosis

The problem of actuator fault diagnosis considered in this paper is to characterize actuator parameter changes (typically gain loss coefficients associated with Φ(k; M j ) as expressed in ( 5)) by estimating them from the mode prior probabilities µ j (0), the mode transition probabilities p i,j , the mode-dependent matrices

A(M j ), B(M j ), C(M j ), Q(M j ), R(M j ), Φ k (M j )
, and the input-output data sequences u(k), y(k).

In this considered framework, the actual active mode M j at each time instant k is unknown. This is the main cause of difficulty for stochastic hybrid system fault diagnosis compared to the case of single mode (non hybrid) systems.

As a matter of fact, if the active mode was known at each time instant k, then the considered hybrid system would be equivalent to a time varying state-space system, for which the actuator fault diagnosis problem, similar to the one formulated in this paper, has already been studied. See, for instance, [START_REF] Zhang | Statistical detection and isolation of additive faults in linear time-varying systems[END_REF]).

INTERACTING MULTIPLE MODEL ESTIMATOR FOR STATE ESTIMATION

In this section, one of the two basic elements for developing the proposed fault diagnosis method, the well-known Interacting Multiple Model (IMM) estimator for hybrid system state estimation, is shortly recalled.

The problem of state estimation, for the considered faultfree (θ = 0 in (3)) stochastic hybrid system, is to characterize the probability distribution of the state vector x(k), from the mode prior probabilities µ j (0), the mode transition probabilities p i,j , the mode-dependent matrices A(M j ), B(M j ), C(M j ), Q(M j ), R(M j ), and the inputoutput data sequences u(k), y(k).

Optimal multiple model estimator and simplifications

In general, mode transitions can happen at every time instant k. To ease the introduction of multiple model estimators, let us first consider a much simpler case: the considered system always remains in one of the r possible modes, but the actual mode is unknown. This is not really a hybrid system and will be referred to as a static unknown mode system. In this case, as the actual mode is unknown, we have to try each of them. Under the assumption of each possible mode, say M j , the system is characterized by the faut-free (θ = 0) stochastic state-space model (3), to which the Kalman filter can be applied for stat estimation. Therefore, at every time instant, r Kalman filters are run in parallel, each assuming a particular working mode. Each of these Kalman filters provides a state estimation. The overall state estimation can be made from a weighted average of the r state estimates with weighting coefficients equal to the posterior probabilities of the r possible modes given input-output observations. This algorithm is known as the static multiple model estimator (Bar-Shalom et al. ( 2001)). Now let us consider true hybrid systems with mode transitions that may happen at each time instant k. In this case, in principle it is no longer sufficient to run r Kalman filters.

As mode transitions can happen at every time instant, all the possible mode sequences up to the current instant k should be considered. Within the instants from 1 to k, there are r k different possible mode sequences. In principle, for the optimal dynamic multiple model estimator, r k Kalman filters should be run in parallel, each corresponding to one of the possible mode sequences. Following this approach, the number of Kalman filters increases exponentially with k. In practice it is not reasonable to implement such solutions, it is then necessary to make simplifications, leading to heuristic solutions, notably the IMM estimator.

The IMM estimator

At each time instant k = 1, 2, 3, . . . , the IMM estimator performs one iteration composed of the same computation steps based on the input-output data, on the modedependent system matrices, and for k = 1 on the initialization data or for k > 1 on the results of the last iteration.

The computation steps at iteration k are as follows.

• Compute the mixing probabilities for mixing the state estimates and covariances of the last iteration. • Compute the mixed state estimates and covariance matrices from the last iteration of the r parallel Kalman filters weighted by the mixing probabilities. • For each of the r assumed active modes at instant k, perform a Kalman filter iteration delivering a state estimate and a covariance matrix. • Update mode probabilities from the likelihoods of the r Kalman filters. • Deliver the algorithm output of the current iteration, by averaging the estimates of the r Kalman filters weighted with the updated mode probabilities.

See (Bar-Shalom et al. ( 2001)) for more details about the IMM estimator.

TIME VARYING SYSTEM ADAPTIVE KALMAN FILTER

In this section, the second basic element for developing the proposed fault diagnosis method, the adaptive Kalman filter, is shortly introduced.

Consider time varying state-space systems in the form of

x(k) = A(k)x(k -1) + B(k)u(k -1) + w(k) + Φ(k -1)θ (6a) y(k) = C(k)x(k) + v(k).
(6b) Here it is assumed that the time varying system matrices A(k), B(k), C(k), Φ(k -1) and noise covariance matrices Q(k), R(k) are known at every time instant k. The stochastic hybrid systems formulated in Section 2 are also statespace systems with system and covariance matrices evolving in time, but the evolution characterized by a Markov model is not known exactly at each time instant.

It is also assumed for system (6) that the initial state x(0) is a random vector following the Gaussian distribution N (x(0|0), P (0|0)), and an initial guess of the unknown parameter vector θ, namely θ(0), is given.

The adaptive Kalman filter is designed for joint estimation of the state vector x(k) and the parameter vector θ of system (6). If the parameter vector θ ∈ R p was known in (6), then the basic Kalman filter would be applicable to its state estimation. In order to perform joint state-parameter estimation, the adaptive Kalman filter presented below incorporates a parameter estimation mechanism, complimenting the part originating from the basic Kalman filter.In addition to the variables originating from the basic Kalman filter, the new algorithm involves new matrices

Υ(k) ∈ R n×p , Ω(k) ∈ R m×p , S(k) ∈ R p×p , Γ(k) ∈ R p×m
and a forgetting factor λ ∈ (0, 1). In particular, the recursively computed Υ(k) and S(k) need initializations Υ(0) = 0 (the n × p zero matrix) and S(0) = αI p with some α > 0.

At each time instant k, the adaptive Kalman filter performs the following computations.

P (k|k-1) = A(k)P (k-1|k-1)A T (k) + Q(k) (7a) Σ(k) = C(k)P (k|k -1)C T (k) + R(k) (7b) K(k) = P (k|k -1)C T (k) [Σ(k)] -1 (7c) P (k|k) = (I n -K(k)C(k))P (k|k -1) (7d) Υ(k) = (I n -K(k)C(k))A(k)Υ(k -1) + (I n -K(k)C(k))Φ(k -1) (7e) Ω(k) = C(k)A(k)Υ(k -1) + C(k)Φ(k -1) (7f) S(k) = 1 λ S(k -1) - 1 λ S(k -1)Ω T (k) λΣ(k) + Ω(k)S(k -1)Ω T (k) -1 Ω(k)S(k -1) (7g) Γ(k) = S(k)Ω T (k) λΣ(k) + Ω(k)S(k)Ω T (k) -1 (7h) ỹ(k) = y(k) -C(k) A(k)x(k -1|k -1) -B(k)u(k -1) -Φ(k -1) θ(k -1) (7i) θ(k) = θ(k -1) + Γ(k)ỹ(k) (7j) x(k|k) = A(k)x(k -1|k -1) + B(k)u(k -1) + Φ(k -1) θ(k -1) + K(k)ỹ(k) + Υ(k)( θ(k) -θ(k -1)). ( 7k
)
Similar algorithms corresponding to the continuous-time counterpart of this adaptive Kalman filter are more frequently studied [START_REF] Zhang | Adaptive observer for multiple-inputmultiple-output (MIMO) linear time varying systems[END_REF]; [START_REF] Li | An adaptive observer for joint estimation of states and parameters in both state and output equations[END_REF]). Because of the noises in the state and output equations, the state and parameter estimation errors do not converge to zero when k tends to infinity. If the noises are ignored, the estimation errors converge exponentially to zero [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF]). In the presence of the noises, by taking mathematical expectations of all the terms in the error equations, the mathematical expectations of the estimation errors converge exponentially to zero.

ADAPTIVE IMM ESTIMATOR FOR HYBRID SYSTEM ACTUATOR FAULT DIAGNOSIS

Now let us go back to the problem of hybrid system actuator fault diagnosis formulated in Section 2 through equations ( 1)-( 4).

If in (3) the mode-dependent matrices A(M j ), B(M j ), C(M j ), Q(M j ), R(M j ), Φ(k; M j ) were known at each time instant (in other words, if the active mode was known at every time instant), then the system would fit into the framework of time varying systems formulated in ( 6), and it would be possible to directly apply the adaptive Kalman filter. Of course, in the presently considered case of stochastic hybrid systems, the active mode is unknown. Like in the case of optimal state estimator discussed in Section 3.1 for fault-free hybrid systems, in principle it is possible to try all the possible mode sequences up to the current instant k and to combine somehow all the resulting state and parameter estimates. Again the number of possible mode sequences (r k ) increases exponentially with time.

Let us follow the same ideas as in the basic IMM estimator to avoid the exponentially increasing complexity of the optimal estimator. At each time instant, instead of considering all the possible past mode sequences, summarize the past estimates with weighted averages by using mixing probabilities µ i|j (k -1|k -1), and run r parallel adaptive Kalman filters, each assuming a different currently active mode M j . The state and parameter estimates of the hybrid system are then obtained by taking the weighted average of the r estimates delivered by the r adaptive Kalman filters, with weights equal to the posterior probabilities of the corresponding modes.

At each instant k, the adaptive IMM estimator for joint state-parameter estimation consists of the following steps.

Calculation of mixing probabilities

The mixing probabilities for mixing results of the last iteration (of instant k -1) are computed as

µ i|j (k -1|k -1) = 1 cj p i,j µ j (k -1) (8)
where µ j (0) (for instant k = 1) are mode prior probabilities, µ j (k-1) (for k > 1) are mode probabilities updated at instant k -1 (see the Mode probability update step below), and cj = r i=1 p i,j µ i (k -1).

Intermediate results mixing

During the last iteration (at instant k -1), r adaptive Kalman filters were run in parallel, each assuming a different active mode at instant k -1, yielding state estimates xi (k-1|k-1), parameter estimates θi (k-1), the matrices P i (k -1|k -1), S i (k -1), Υ i (k -1), all indexed by the corresponding assumed mode M i .

The mixed quantities are then computed as

xj (k -1|k -1) = r i=1 xi (k -1|k -1)µ i|j (k -1|k -1) θj (k -1) = r i=1 θi (k -1)µ i|j (k -1|k -1) Ŝj (k -1) = r i=1 S i (k -1)µ i|j (k -1|k -1) Υj (k -1) = r i=1 Υ i (k -1)µ i|j (k -1|k -1) P j (k -1|k -1) = r i=1 µ i|j (k -1|k -1){P i (k -1|k -1) + [x i (k -1|k -1) - xj (k -1|k -1)] • [x i (k -1|k -1) - xj (k -1|k -1)] T }.

Mode-matched filtering

For each of the possible modes M j at instant k, an adaptive Kalman filter is implemented as follows.

P j (k|k-1) = A(M j ) P j (k-1|k-1)A T (M j )+Q(M j ) Σ j (k) = C(M j )P j (k|k -1)C T (M j ) + R(M j ) K j (k) = P j (k|k -1)C T (M j ) Σ j (k) -1 P j (k|k) = (I n -K j (k)C(M j ))P j (k|k -1) Υ j (k) = (I n -K j (k)C(M j ))A(M j ) Υj (k -1) + (I n -K j (k)C(M j ))Φ(k -1; M j ) Ω j (k) = C(M j )A(M j ) Υj (k -1) + C(M j )Φ(k -1; M j ) S j (k) = 1 λ Ŝj (k -1) - 1 λ Ŝj (k -1)[Ω j (k)] T • λΣ j (k) + Ω j (k) Ŝj (k -1)[Ω j (k)] T -1 • Ω j (k) Ŝj (k -1) Γ j (k) = S j (k)[Ω j (k)] T λΣ j (k) + Ω j (k)S j (k)[Ω j (k)] T -1 ỹj (k) = y(k) -C(M j ) A(M j ) xj (k -1|k -1) -B(M j )u(k-1)-Φ(k-1; M j ) θj (k-1) θj (k) = θj (k -1) + Γ j (k)ỹ j (k) xj (k|k) = A(M j ) xj (k -1|k -1) + B(M j )u(k -1) + Φ(k -1; M j ) θj (k -1) + K j (k)ỹ j (k) + Υ j (k)( θj (k) - θj (k -1)
). The likelihood of the mode M j , given the input-output data up to instant k, is evaluated through the innovation ỹj (k)

Λ j (k) = 1 (2π) m det(Σ j (k)) • exp - 1 2 (ỹ j (k)) T (Σ j (k)) -1 ỹj (k) .

Mode probability update

The mode probabilities are updated as

µ j (k) = 1 c Λ j (k)c j with cj = r i=1 p i,j µ i (k -1) c = r j=1 Λ j (k)c j .

Algorithm outputs

The outputs of the adaptive IMM estimator at instant k are the state estimate Like the basic IMM estimator, this adaptive algorithm is also based on heuristic simplifications of the optimal estimator to avoid the exponentially increasing complexity.

NUMERICAL EXAMPLE

In order to illustrate the proposed solution for actuator fault diagnosis in hybrid systems, let us consider a system with 4 modes (r = 4 ). Each mode is described by a third order state-space model (n = 3) subject to white Gaussian noises, with one input (l = 1) and 2 outputs (m = 2). The mode-dependent system matrices A(M j ), B(M j ), C(M j ) are randomly generated such that each mode is stable (the eigenvalues of A(M j ) are inside the unit cycle), observable and controllable. The noise covariances matrices are chosen as Q(M j ) = 0.1I 3 and R(M j ) = 0.05I 2 for all the 4 modes. The mode transition probabilities p i,j are randomly generated. During each simulation trial, a gain loss of 50%, corresponding to a jump of θ (a scalar parameter, i.e., p = 1) from 0 to 0.5 at the time instant k = 500 is simulated, and the simulation runs till k = 1000. The adaptive IMM estimator proposed in this paper is then applied to the simulated system for joint state-parameter estimation. The result of parameter estimation for one of the trials is presented in Fig. 1, and that of state estimation in Fig. 2.

In these figures, the results obtained with the adaptive Kalman filter (see Section 4) are also presented as a reference for the purpose of comparison. In practice the adaptive Kalman filter is not applicable, as its computations require the true mode transition sequence, which is unknown in practice. Because the adaptive IMM estimator uses less information, it cannot perform as well as the adaptive Kalman filter, but yet the results are quite similar.

In order to statistically evaluate the performance of the proposed method, 1000 simulated trials are performed, each corresponding to a different random realization of mode transition probabilities, mode-dependent system matrices, state and output noises. At each time instant k, the histogram of the parameter estimation error based on the 1000 simulated trials is generated, and all the histograms are depicted as a 3D illustration in Fig. 3 for the adaptive IMM estimator, and in Fig, 4 for the adaptive Kalman filter, again for the purpose of comparison. The histograms are normalized so that they are similar to probability density functions.

CONCLUSION

Based on the existing basic IMM estimator for hybrid system state estimation and on the adaptive Kalman filter for time varying system joint state-parameter estimation, a new algorithm, the adaptive IMM estimator has been proposed in this paper for actuator fault diagnosis in stochastic hybrid systems. In the presented numerical illustrations, the results of the adaptive IMM estimator are quite close to those of the adaptive Kalman filter, which represents a performance upper bound that cannot be attained in practice. The computational burden of the adaptive IMM estimator is essentially equal to that of the adaptive Kalman filter multiplied by the number of modes.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The simulated "true" parameter (θ, black), parameter estimates by the adaptive IMM estimator (adIMM, blue) and by the adaptive Kalman filter (AdKF, red).

Fig. 3 .

 3 Fig. 3. Histogram per instant k of the parameter estimation error of the adaptive IMM estimator.

For shorter statements, "the mode labeled by M j " is often simply written as "M j " in this paper.


F Loete

Q Zhang

M Sorine

Experimental Evaluation of the Inverse Scattering Method for Electrical Cable Fault Diagnosis

Keywords: Inverse scattering, cable health monitoring, soft fault diagnosis, frequency domain reflectometry. 

Recently published theoretic and experimental results have shown the ability of inverse scatteringbased methods to detect and to locate soft faults in electric cables, in particular, faults implying smooth spatial variations of cable characteristic parameters. The purpose of the present paper is to further experimentally evaluate the inverse scattering method for retrieving spatially distributed characteristic impedance from reflectometry measurements. With high quality coaxial cables connected in parallel, composite cables of piecewise constant characteristic impedance profiles are built in order to evaluate the accuracy of the inverse scattering method and its robustness in the presence of impedance discontinuities.

INTRODUCTION

The number of wired connections in modern engineering systems has been rapidly increased during the last decades, to the point that the reliability of electrical connections becomes a crucial issue. For example, in a modern automotive vehicle, the total length of on-board cables can reach up to 4 km. These cables, composed of twisted pairs, coaxial cables and simple wires, are subject to aging or degradations caused by severe environmental conditions. In this area, reliability becomes a safety issue. In some other domains, cable defects may have catastrophic consequences [START_REF] Furse | A Critical Comparison of Reflectometry Methods for Location of Wiring Faults[END_REF]). It is thus a crucial challenge to design smart embedded diagnosis systems able to detect wired connection defects in real time. This fact has motivated research projects on methods for fault diagnosis in electric transmission lines.

Many diagnosis methods based on Time or Frequency Domain

Reflectometry (TDR or FDR) have been developed to detect, locate and characterize defects in electrical harnesses (Furse et al. [2005], [START_REF] Lelong | On line wiring diagnosis by modified spread spectrum time domain reflectometry[END_REF], [START_REF] Smail | Detection of Defects in Wiring Networks Using Time Domain Reflectometry[END_REF], [START_REF] Furse | Frequency domain reflectometry for on board testing of aging aircraft wiring[END_REF], [START_REF] Smith | Fault location on aircraft wiring using spread spectrum time domain reflectometry[END_REF], [START_REF] Furse | An inexpensive distance measuring system for navigation of robotic vehicles[END_REF], [START_REF] Chung | Application of phase detection frequency domain reflectometry for locating faults in an F-18 flight control harness[END_REF], Loete et al. [2012a]). These methods are based on the same principle as radar technology, but applied to guided waves in electric cables. A high frequency electrical signal is sent down a cable, where it is reflected at impedance discontinuities or at other inhomogeneities along the cable. The reflected wave is analyzed to detect, locate and characterize the defects. It has been reported that this technology is able to detect and to locate hard faults (open or short circuits) up to an accuracy of about 3 cm [START_REF] Smith | Fault location on aircraft wiring using spread spectrum time domain reflectometry[END_REF]). On the other hand, it is a much more difficult problem to deal with soft faults which imply only slight and localized characteristic impedance changes in electric cables. Some recently reported methods amplify fault signatures to improve the detection of the small echoes caused by soft faults [START_REF] Franchet | The use of the pseudo wigner-ville transform for detecting soft defects in electric cables[END_REF], [START_REF] Franchet | On a useful tool to localize jacks in wiring network[END_REF]], El Sahmarany et al. [2012a,b], [START_REF] References Abboud | A matched-pulse approach for soft-fault detection in complex wire networks[END_REF], [START_REF] Abboud | A non-iterative method for locating soft faults in complex wire networks[END_REF]). However, to our knowledge, few satisfactory experimental results on reflectometry-based methods for the diagnosis of such soft faults have been reported (Loete et al. [2012a]), (Loete et al. [2012b]).

The characteristic impedance

is the main parameter characterizing an electric cable. When a healthy and uniform cable is considered, the characteristic impedance is usually specified as a single value for the whole cable. In order to deal with localized defects, it is necessary to consider distributed characteristic impedance all along the cable, denoted as with indicating the position along the cable in some coordinate system.

In traditional reflectometry-based methods, the reflected waveforms are visualized and analyzed in the time domain (TDR) or in the frequency domain (FDR), and sometimes a mixture of them is used. As these methods do not directly analyze the characteristic impedance of the cable, it may be difficult to interpret their results for soft faults diagnosis, in particular when such faults cause only small variations of the characteristic impedance along the cable. It is thus useful to develop inverse methods retrieving the spatially distributed characteristic impedance profile from reflectometry measurements (in opposition to direct methods simulating reflectometry measurements from specified cable characteristics), but few methods of this nature have been reported [START_REF] Smail | Detection and location of defects in wiring networks using Time Domain reflectometry and neural networks[END_REF]). The inverse scattering theory [START_REF] Jaulent | The inverse scattering problem for LCRG transmission lines[END_REF]) provides a powerful theoretic framework for the estimation of distributed characteristic impedance from reflectometry measurements, in particular for smoothly varying characteristic impedance profiles along cables. Based on this theory, the ability of directly visualizing and analyzing the distributed characteristic impedance of a cable represents an important progress in the interpretation and processing of reflectometry measurements, notably for the purpose of soft fault diagnosis. The theoretic study and numerical simulations reported in [START_REF] Zhang | Inverse scattering for soft fault diagnosis in electric transmission lines[END_REF] and [START_REF] Tang | An inverse scattering approach to soft fault diagnosis in lossy electric transmission lines[END_REF] have shown that the inverse scattering method is a promising tool for soft fault diagnosis through the estimation of distributed characteristic impedance. The experimental results reported in Loete et al. [2012b] have confirmed the ability of this method for detecting and locating soft faults. In these experiments, the tested twisted pair cable is untwisted at a small portion in order to create a smooth local variation of the characteristic impedance. Such experiments provide a simple solution for hard wire simulation of soft faults, but it is not easy to accurately control the characteristic impedance profile of the untwisted portion in order to evaluate the accuracy of the inverse scattering method.

The purpose of the present paper is to evaluate the accuracy of the method through new experiments designed with parallel coaxial cables whose characteristic impedances are known with a good accuracy. As in theory the inverse scattering method assumes smooth characteristic impedance profiles, the piecewise constant characteristic impedance profiles used in the experiments realized with coaxial cables allow to evaluate the robustness of the method when the smoothness assumption is not satisfied. This paper is organized as follows. After briefly recalling in Section 2 the inverse scattering method for electric cable monitoring, experimental results with coaxial cables will be presented in Section 3 to evaluate the accuracy and the robustness of the inverse scattering method. Section 4 then concludes the paper.

THE INVERSE SCATTERING METHOD

This section shortly recalls the inverse scattering method for lossless cables. More details can be found in [START_REF] Zhang | Inverse scattering for soft fault diagnosis in electric transmission lines[END_REF]. For the electric cables of less than 10 meters studied in this paper, losses are neglected. A lossless cable driven by a harmonic voltage wave can be modeled by the frequency domain telegrapher's equations

, , 0 (1a)
, , 0 (1b)
where is the angular frequency (denoted so because it is strongly related to the wavenumber in the inverse scattering theory), indicates the position along the cable, , and , denote the voltage and the current of frequency at the position , and are distributed inductance and 1 is a monotonic (and thus invertible) function, corresponding to the inverse function . By abuse of capacitance along the cable; is the imaginary unit. The boundary conditions associated to (1) represent a generator at the left end and a load at the right end of the cable.

Define the electrical distance, the characteristic impedance, the reflected and direct power waves, as follows 1 , ∈ 0,

(2) ν , , , (3a)
ν , , , (3b)
Some direct computations then lead to the following Zakharov-Shabat equations [START_REF] Jaulent | The inverse scattering problem for LCRG transmission lines[END_REF]):

ν , ν , ν , (4a) ν , ν , ν , (4b) where (5)
At the left end (corresponding to 0), the reflection coefficient ν , 0 /ν , 0 is typically measured with a network analyzer. It is known [START_REF] Zhang | Inverse scattering for soft fault diagnosis in electric transmission lines[END_REF]) that the characteristic impedance can be computed from measured at the left end of the cable through the following steps:

1. Compute the Fourier transform of the reflection coefficient 1 2 exp 2. Solve the integral equations (known as Gel'fand-Levitan-Marchenko equations) for their unknown kernels , and , :

, , 0 , , 0
3. Compute the potential function through notation, V , will be simply written as V , , so will other similar notations depending on .

, 4. By solving equation ( 5) for , compute 0 2 where 0 is equal to the internal impedance of the signal generator connected to the left end of the cable (typically integrated in a network analyzer).

EXPERIMENTAL EVALUATION OF THE INVERSE SCATTERING METHOD

Since equation ( 4) contains the spatial derivative of the characteristic impedance, it is assumed in theory that is derivable. Consequently, the inverse scattering method presented in section II is, in theory, limited to smooth characteristic impedance profiles, thus excluding the case of impedance discontinuities. The ability of the inverse scattering method for detecting and locating characteristic impedance smooth variations has been experimentally validated Loete et al. [2012b]. However, impedance discontinuities are frequently present in practice, typically at the connection between the network analyzer and the tested cable if there is an impedance mismatch or at the junction between several branches. The aim of this section is to study the robustness of the inverse scattering method in presence of impedance discontinuities and to quantify its accuracy in such situations.

The experiments presented in this section are designed with 50 Ω coaxial RG-58 cables with low loss characteristics. They also have good noise immunity, and a well-controlled spatially homogeneous 50 Ω characteristic impedance. The tested configurations are very simple, but they are in no way favorable to the inverse scattering method, because of impedance discontinuities. The accuracy of the characteristic impedance of the coaxial cables allows us to quantitatively evaluate the error of the inverse scattering method. The measurements are carried out over the 4 MHz -2204 MHz bandwidth with 4 MHz frequency steps.

Coaxial cables will be used to build composite lines with piecewise constant characteristic impedance profiles. However, the applied inverse scattering algorithm is exactly the same one as in the previous section, which does not assume piecewise constant characteristic impedance profiles.

Remark that this study is for the purpose of evaluating accuracy and the robustness of the inverse scattering method at the presence of the impedance discontinuities. The fact that the tested cable has a piecewise constant characteristic impedance profile is ignored by the inverse scattering method, which has been designed for arbitrary smooth characteristic impedance profiles. In practice if it is known in advance that the cable has a piecewise constant characteristic impedance profile, some other methods designed for this particular case may be more appropriate. However, the piecewise constant profile assumption is not reasonable if the tested cable is possibly affected by soft faults.

Characteristic impedance drop on a single segment

Two identical 50 Ω coaxial cables in parallel (Fig. 1) are used as an equivalent cable of 25 Ω. Similarly, three and four identical cables in parallel form respectively equivalent 16.67 Ω and 12.5 Ω cables. Such a composite segment is then connected to simple 50 Ω coaxial cables as shown in Fig. 1 and illustrated in Fig. 2, in order to build a cable with piecewise constant characteristic impedance. The advantage of this configuration is the accuracy of the characteristic impedance (50 Ω, 25 Ω, etc.) in each segment, based on the high quality coaxial cables used in the experiment. Fig. 3 shows the inverse scattering results carried out on the cables illustrated in Fig. 2. The location of the impedance discontinuities, as well as the impedance level of the composite segment, are estimated with a good accuracy. Nonetheless, one can see some artifacts in the part immediately following the segment composed of parallel cables, especially for strong impedance discontinuities. It is possible to remedy this drawback by increasing the bandwidth of the reflectometry measurements, but in practice the network analyzer has a limited bandwidth. By adding zeros to complete high frequency measurements (zero padding), these artifacts are also reduced, because the inverse scattering algorithm then uses a smaller integration step size. For example, in one of the experiments, the same reflection coefficient is processed twice by the inverse scattering method, first with zero padding up to 15 GHz and then up to 30 GHz. The results are shown in Fig. 4 confirming that a higher zero padding improves the result. However, the practice of zero padding is limited by the computer memory size when the data are processed. The results shown in Fig. 3 have been obtained with zero padding up to 30 GHz. Fig. 3 (z) profiles obtained by inverse scattering on the line including an (equivalent) segment of (-) 50 Ω (--) 25 Ω (▪▪) 16.6 Ω and (▪ ─) 12.5 Ω.

Fig. 4 Comparison of two (z) profiles obtained by inverse scattering on the same line including an (equivalent) segment of 25 Ω with zero padding up to 15 GHz (-) and 30 GHz (--). These results show that a higher zero padding leads to a more accurate result

Double characteristic impedance drop

This example is intended to study the ability of the inverse scattering method to estimate more complex discontinuous characteristic impedance profiles.

In Fig. 5, the estimated characteristic impedance profile exhibits spikes at each junction between the different segments of the line (e.g at the injection point, at the connection of the 50 Ω coaxial cable and the segment of doubled cables etc…). These are caused by the 50 Ω connectors which are slightly improperly adapted to the line and furthermore locally add some extra lengths to the line (each straight or Tee connector adds about 2cm to the 50 Ω coaxial line). This explains the small two first impedance spikes observed at 0 and 0.6 . The spike at the third connection point (between the 25 Ω and the 16.66 Ω segments) is all the more observed since 3 BNC Tees are stacked in order to connect the doubled and tripled cables segments.

As shown in Fig. 5, in the middle of a 50 Ω coaxial cable, a segment of doubled cables and another one of tripled cables are inserted, with equivalent characteristic impedances of 25 Ω and 16.67 Ω. The result of the inverse scattering method is also shown in Fig. 5. The characteristic impedance profile is estimated with a good accuracy, despite the multiple discontinuities present in the characteristic impedance profile

Accuracy of the inverse scattering method in the case of piecewise constant characteristic impedance profiles.

Remind that, in theory, the case of characteristic impedance with discontinuities is not suitable for the inverse scattering method. Nevertheless, the reported experiments show that this method can tolerate discontinuities to some extent. The results of section 3.1 show that, the higher is the step jump of the characteristic impedance, the larger is the error of the inverse scattering method. To complete these results, numerical simulations have been made to further characterize the errors. The simulated cable configurations are similar to that of Fig. 2, but the characteristic impedance level jump in the middle segment can be freely chosen for numerical simulations. The reflection coefficients of the simulated cables with various impedance level jumps ΔZ are simulated in the 50Ω ; 350Ω range and then are given as input to the inverse scattering method. The relative accuracy of the inverse scattering method was calculated as:

% ∆ ∆ ∆
where ∆ is the impedance level jump of the middle segment estimated by the inverse scattering method. The results are plotted in Fig. 6. As expected, the error is increasing with ΔZ.

Therefore, in presence of strong impedance jumps, the error term can rapidly reach quite high values. Nevertheless, when ∆ ∈ 50Ω ; 100Ω , the error is less than 5%. So one must be careful when interpreting impedance values computed in presence of strong discontinuities. When both smooth variations and discontinuities are present in the characteristic impedance profile of a cable, the reflection coefficient measured at one end of the cable is the result of all these imperfections of the cable. As waves injected in a cable are essentially reflected at the discontinuities, too strong discontinuities may mask the effects of the reflections caused by smooth variations. It is thus important to avoid strong discontinuities when smooth variations should be monitored. In the experiments reported in Loete and al. [2012b], it was shown that an impedance jump of 60 Ω (impedance mismatch between the network analyzer and the 110 Ω cable), and another one of 110 Ω (impedance mismatch between the 110 Ω cable and the 220 Ω load) are well tolerated when smoothly varying characteristic impedance profiles are estimated, with a satisfactory accuracy for the purpose of soft fault detection, location and characterization.

CONCLUSION

The experimental results presented in this paper confirm the previously reported theoretic results and numerical simulations about the inverse scattering method. These results provide a new way for the interpretation of reflectometry measurements. To characterize the errors of the inverse scattering method, high quality coaxial cables were then used to build composite lines with piecewise constant characteristic impedance profiles. It is shown that the inverse scattering method can tolerate to some extent impedance discontinuities.

The results reported in this paper, together with those of Loete et al. [2012b], confirm the ability of the inverse scattering method for distributed characteristic impedance estimation, with a satisfactory accuracy for the purpose of soft fault diagnosis in electric cables. Moreover, the implemented inverse scattering algorithm is numerically efficient: for each of the experiments presented in this paper, the processing of the data by the inverse scattering algorithm takes less than a second on a typical notebook computer.

Fig. 1 .

 1 Fig. 1. Composite line including two RG-58 50 Ω coaxial cables in parallel.

Fig. 2 .

 2 Fig. 2. Topology of the composite lines used. Equivalent cables of characteristic impedance smaller than 50 Ω are realized by putting multiple 50Ω coaxial cables in parallel.

Fig. 5

 5 Fig. 5 Above: Two segments of doubled and tripled cables inserted in the middle of a 50 Ω coaxial. Below: Result of the inverse scattering method.

Fig. 6

 6 Fig. 6 Inverse scattering method accuracy in case of impedance discontinuities. (▲)Accuracy obtained from the inversion of experimentally measured reflection coefficient. (■) Accuracy obtained from the inversion of simulated reflection coefficient.


Keywords: Hybrid system, descriptor system, IMM estimator, descriptor system Kalman filter

Stochastic hybrid systems have been largely studied in the literature in the framework of Markov mode transitions and linear Gaussian state space mode models. In order to generalize such results to hybrid systems involving phenomena characterized by differentialalgebraic equations, this paper studies the problem of state estimation for stochastic hybrid systems with modes described by descriptor equations. The proposed state estimation algorithm follows the interacting multiple model (IMM for short) approach, but the classical state space system Kalman filters are replaced by descriptor system Kalman filters. Because of the difficulty for computing the innovation of each descriptor system Kalman filter, a new method for likelihood evaluation is proposed, as one important step of the new IMM algorithm. Numerical examples are presented to illustrate the performance of the proposed algorithm.

Hybrid Descriptor System State Estimation through an IMM Approach

Liangquan Zhang * , * * Qinghua Zhang * , * * *

INTRODUCTION

Differential equations have been widely used in the study of dynamic systems, in particular, finite dimensional state space equations of the form ẋ(t) = f (x(t), u(t))

(1) are often used for modeling engineering systems, with x(t) and u(t) denoting respectively the state vector and the input vector of the considered system, ẋ(t) = dx(t)/dt, and f is a function characterizing the dynamic behavior of the system. Despite the large success of such theories in engineering practice, some complex systems cannot be appropriately described in this framework. Such exceptions include differential-algebraic systems and hybrid systems.

Algebraic constraints in engineering systems typically result from singularities of differential equations. For example, a train following a turning track is constrained by the geometrical form of the tack if the mass of the earth is considered infinitely large compared to the mass of the train. Such a system can be described by differentialalgebraic equations (DAE) of the form g( ẋ(t), x(t), u(t)) which can represent a wider class of systems than the classical state space systems of the form (1). After linearization and discretization in time, such equations can be approximated by implicit discrete time state space equations of the form

E k+1 x(k + 1) = A k x(k) + B k u(k) + υ(k),
(2) where x(k), u(k) and υ(k) are respectively the discrete time state, input and the modeling errors indexed by k = 1, 2 . . . , and E k+1 , A k , B k are time varying matrices of appropriate sizes. With a possibly rank deficient matrix This work has been supported by the ITEA2 MODRIO project. E k+1 , refer to (2) which is known as a descriptor equation [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]].

On the other hand, some complex systems have different working modes, for example, the starting mode, the normal working mode, or some reduced regime mode in case of component failures. If in each of these modes the system is described by some differential equations, the over-all functioning of the system, including the mode switching mechanism and the behavior within each mode, is usually modeled as a hybrid system [ [START_REF] Bar-Shalom | Estimation with applications to tracking and navigation[END_REF]].

The purpose of the present paper is to study state estimation for hybrid systems with working modes described by descriptor equations of the form (2).

The study of descriptor systems, notably about state estimation, has a rich literature [START_REF] Yeu | Fault detection and isolation for descriptor systems using sliding mode observer[END_REF]; [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]; [START_REF] Marx | Robust faulttolerant control for descriptor systems[END_REF]; [START_REF] Koenig | Design of proportional-integral observer for unknown input descriptor systems[END_REF]; [START_REF] Gao | Descriptor observer approaches for multivariable systems with measurement noises and application in fault detection and diagnosis[END_REF]] (for more information see references therein). Besides this topic, state estimation for hybrid systems has been largely studied with various applications such as target tracking [START_REF] Bar-Shalom | Tracking and data association[END_REF]; [START_REF] Blom | The interacting multiple model algorithm for systems with Markovian switching coefficients[END_REF]], signal processing [START_REF] Doucet | Particle filters for state estimation of jump Markov linear systems[END_REF]], and fault diagnosis [START_REF] Hanlon | Multiple-model adaptive estimation using a residual correlation Kalman filter bank[END_REF]; [START_REF] Koutsoukos | Monitoring and diagnosis of hybrid systems using particle filtering method[END_REF]]. These reported results concern hybrid systems involving modes described by state space equations. To our knowledge, no study has been reported about hybrid systems with working modes described by descriptor equations.

In the case of stochastic hybrid systems with working modes described by classical linear Gaussian state space equations, it is known that the complexity of the optimal state estimator increases exponentially with time. In practice, heuristic algorithms of reduced complexity are used, without convergence proof. Among such algorithms, is the interacting multiple model (IMM for short) estimator [START_REF] Bar-Shalom | Estimation with applications to tracking and navigation[END_REF], C. E. [START_REF] Seah | State estimation for stochastic linear hybrid systems with continuous-state-dependent transitions: An IMM approach[END_REF]].

It is based on multiple Kalman filters, each assuming a particular working mode. At each iteration, in addition to the classical prediction and update steps of the Kalman filter, all the Kalman filters are modified based on the knowledge about the mode transition mechanism modeled as a Markov chain.

In the present paper, the algorithm proposed for hybrid descriptor system state estimation follows also the IMM approach. The main novelties in this algorithm reside in the replacement of the classical state space system Kalman filter by the descriptor system Kalman filter, and in a new method for the evaluation of the likelihood of each Kalman filter at each iteration. Of course, like the classical IMM estimator for state-space hybrid systems, no convergence proof of the algorithm is known. This paper is organized as follows: In Section 2 we present a formal formulation of the estimation problem for a stochastic linear hybrid descriptor system. In Section 3, we make a comparison between the Kalman filter for classical situation and the descriptor one. We state the IMM algorithm in the framework of descriptor system in 4. In Section 5, we provide numerical examples. Finally, conclusions are presented in Section 6.

PROBLEM FORMULATION

The stochastic hybrid systems considered in this paper are modeled at two levels. At the top level, each system has a finite number of working modes. At each time instant, one of the modes is active, and random transitions between different modes are characterized by a Markov model. At the bottom level, each mode of the system is described by a stochastic descriptor system model.

Top level Markov Transition Model

Let {M 1 , M 2 , . . . , M r } be the set of r possible modes of a hybrid system, and m (k) ∈ {1, 2, . . . , r} denote the index of the mode at time instant k = 0, 1, 2, . . .. The mode evolution, corresponding to the sequence m(0), m(1), • • • , is a Markov chain described by the mode transition matrix

Π = {p ij } i,j=1,••• ,r , where p i,j = P [ m(k + 1) = j| m(k) = i] , with known transition probabilities p i,j independent of k and satisfying r j=1 p i,j = 1, i = 1, 2, . . . , r.
For the initial time instant k = 0, the prior probability that M j is active is

P [m(0) = j] = µ j (0) with known probabilities µ 1 (0), µ 2 (0), . . . , µ r (0) satisfying r j=1 µ j (0) = 1.

Bottom level Stochastic Descriptor System Model

At the bottom level, the state x(k) ∈ R n , the input u(k) ∈ R q and the output y(k) ∈ R p of the hybrid system satisfy the stochastic descriptor equations

E Mj x(k + 1) = A Mj x(k) +B Mj u(k) + W Mj (k), (3)
y(k) = C Mj x(k) +D Mj u(k) + V Mj (k), (4) where the time index k = 0, 1, 2, • • • ., the descriptor system matrices E Mj ∈ R l×n , A Mj ∈ R l×n , B Mj ∈ R l×q , C Mj ∈ R p×n , D Mj ∈ R p×q ,
and the state and output white noises

W Mj (k) ∈ R l , V Mj (k) ∈ R p , with W Mj (k) ∼ N (0, Q Mj ), V Mj (k) ∼ N (0, S Mj ) which are mutually independent to each other.
When the mode of the hybrid system stays unchanged, say at M j , the mode index m(k) = j, then the system described by equation ( 3) behaves like a stochastic descriptor system characterized by the matrices A j , B j , C j , D j , E j , Q j , S j . Within each mode, it is assumed that descriptor system (3)-( 4) is observable and controllable by the state noise. See [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]] for the definitions of descriptor system observability and controllability.

Other other hand, the mode M j may evolve, so that the mode index sequence m(0), m(1), • • • forms a Markov chain, as described at the top level. In particular, during every mode change, say from M m(k) = M j to M m(k+1) = M i , the evolution of the state x(k+1) and the output y(k+ 1) are also described by equation ( 4), given a realization of the mode sequence Markov chain.

Let Z k = [u(0), y (0) , u(1), y(1), . . . , u(k), y(k)] denote the measurements up to time k.

The following assumptions are made in this paper.

(A1): Assume that x(0) ∼ N (0, P 0 ) is independent of W Mj (k) and V Mj (k);

(A2): For j ∈ {1, 2, . . . , r} , H Mj = E Mj C Mj is full column rank; (A3): For j ∈ {1, 2, . . . , r} , the matrices Q Mj and S Mj are symmetric positive definite.

KALMAN FILTER FOR DESCRIPTOR SYSTEMS

In this section, let us consider the single-mode descriptor system

E k+1 x(k + 1) = A k x(k) + B k u(k) + W (k), (5)
y(k + 1) = C k+1 x(k + 1) + D k+1 u(k + 1)
+V (k + 1). ( 6) First recall the classical Kalman filter for state space systems corresponding to the case l = n and

E k+1 = I n×n .                  x(0|0) = E(x 0 ), P (0, 0) = V ar(x 0 ), P (k + 1|k) = A(k)P (k|k)A T (k) + Q(k), G(k + 1) = P (k + 1|k)C T (k + 1) •(C(k + 1)P (k + 1|k)C T (k + 1) + S(k + 1)) -1 , x(k + 1|k) = A(k)x(k|k) + B(k)u(k), x(k + 1|k + 1) = x(k + 1|k) + G(k + 1)(y(k + 1) -D(k + 1)u(k + 1) -C (k + 1) x(k + 1|k))
, where the innovation process, also known as the prediction error, v(k

+ 1) = y(k + 1) -ŷ(k + 1|k) with ŷ(k + 1|k) = C(k + 1)x(k + 1|k) + D(k + 1)u(k + 1),
is a white Gaussian sequence. The associated likelihood function can be computed by using the innovation v(k+1).

Now let us consider the more general descriptor system (5)-( 6). Its Kalman filter [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]] writes (d in subscript short for descriptor systems),

                                                         xd (0|0) = E(x 0 ), Σ d (k) = A(k)P d (k|k)A T (k) + Q d (k), G d (k + 1) = (E T (k + 1)Σ -1 d (k)E(k + 1)) -1 C T (k + 1) •(C(k + 1)(E T (k + 1)Σ -1 d (k)E(k + 1)) -1 •C T (k + 1) + S d (k)) -1 , L d (k + 1) = (E T (k + 1)Σ -1 d (k)E(k + 1)) -1 •E T (k + 1)Σ -1 d (k) -G d (k + 1)C(k + 1) •(E T (k + 1)Σ -1 d (k)E(k + 1)) -1 •E T (k + 1)Σ -1 d (k), xd (k + 1|k + 1) = L d (k + 1)A(k)x d (k|k) +L d (k + 1)B(k)u(k) -G d (k + 1)D(k + 1)u(k + 1) +G d (k + 1)y d (k + 1), P d (k + 1|k + 1) = L d (k + 1)A(k)P d (k|k)A T (k) •L T d (k + 1) + L d (k + 1)Q d (k)L T d (k + 1) +G d (k + 1)S d (k)G T d (k + 1). (7)
In [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]

(k +1|k +1) = C(k +1)x d (k +1|k +1)+D(k +1)u(k +1).
is not the innovation sequence in the sense of prediction error as defined in the classical Kalman filter. Usually the innovation sequence is used for likelihood evaluation in the IMM approach to hybrid system estimation. Fortunately, like the innovation sequence, this estimation error v d (k+1) is a white Gaussain sequence in [START_REF] Ali | An innovations approach to fault diagnosis in linear time-varying descriptor systems[END_REF]]. This result is important for the new IMM estimator presented in this paper for hybrid systems, at the step of likelihood computation.

HYBRID DESCRIPTOR SYSTEM IMM ESTIMATOR

The algorithm presented below for hybrid descriptor systems is quite similar to the IMM estimator for classical hybrid state space systems [START_REF] Bar-Shalom | Estimation with applications to tracking and navigation[END_REF]], the main differences reside in the Kalman filter for descriptor systems and in the evaluation of the likelihood at each iteration.

The key feature of IMM we point out is that it consists of r interacting filters operating in parallel.

The algorithm consists of the following steps.

1. Calculation of the mixing probabilities (i, j = 1, . . . , r). The probability that mode M i was in effect at k given that M j is in effect at k + 1 conditioned on Z k+1 is

µ i|j (k|k + 1) = P {M i (k) |M j (k + 1), Z k+1 } = 1 cj p i,j µ j (k),
where cj = r i=1 p i,j µ i (k).

2. Intermediate results mixing (j = 1, . . . , r).

During the last iteration, r descriptor system Kalman filters were run in parallel, each assuming a different active mode at instant k, yielding r state estimates x (k|k) and r covariance matrices P i (k|k) . Based on these results, the mixed state estimates and covariance matrices are x0,j (k|k) = r i=1

x(k|k)µ i|j (k|k),

P 0,j (k|k) = r i=1 µ i|j (k|k){P i (k|k) +[x(k|k) -x0,j (k|k)] •[x(k|k) -x0,j (k|k)] T }.
3. Mode-matched filtering (j = 1, . . . , r).

For each of the r assumed active modes at instant k + 1, say M j , a descriptor system Kalman filter delivers a state estimate xj (k + 1|k + 1) as follows:

xj (0|0) = E(x 0 ), Σ j (k) = A Mj P 0,j (k|k)A T Mj + Q Mj , G Mj (k + 1) = (E T Mj (Σ j (k)) -1 E Mj ) -1 C T Mj •(C Mj (E T Mj (Σ j (k)) -1 E Mj ) -1 C T Mj +S Mj ) -1 , L Mj (k + 1) = (E T Mj (Σ j (k)) -1 E Mj ) -1 •E T Mj (Σ j (k)) -1 -G Mj (k + 1)C Mj •(E T Mj (Σ j (k)) -1 E Mj ) -1 E T Mj (Σ j (k)) -1 , xj (k + 1|k + 1) = L Mj (k + 1)A Mj x0,j (k|k) +L Mj (k + 1)B Mj u(k) -G Mj (k + 1)D Mj u(k + 1)
+G Mj (k + 1)y(k + 1),

P j (k + 1|k + 1) = L Mj (k + 1)A Mj P 0,j (k|k)A T Mj •L T Mj (k + 1) +L Mj (k + 1)Q Mj L T Mj (k + 1) +G Mj (k + 1)S Mj G T Mj (k + 1
). The likelihood of the mode M j , given the input-output data up to instant k + 1, is evaluated through the "innovation"

vj (k + 1) = y(k + 1) -ŷj (k + 1|k + 1) = y(k + 1) -D Mj u(k + 1) -C Mj xj (k + 1|k + 1), Λ j (k + 1) = 1 (2π) p det(Σ j (k + 1)) • exp(Υ j (k + 1)),
where

Υ j (k + 1) = - 1 2 vj (k + 1) T Σ j (k + 1) -1 vj (k + 1).
4. Mode probability update.

The mode probabilities are updated as

µ j (k + 1) = 1 c Λ j (k + 1)c j with cj = r i=1 p i,j µ i (k), c = r j=1 Λ j (k + 1)c j .

Estimate and covariance combination.

Combination of the model-conditioned estimates and covariances is completed according to the mixture equations

x(k + 1|k + 1) = r j=1 xj (k + 1|k + 1)µ j (k + 1)
and

P (k + 1|k + 1) = r j=1
µ j (k + 1) P j (k + 1|k + 1)

+ xj (k + 1|k + 1) -x(k + 1|k + 1)

• xj (k + 1|k + 1) -x(k + 1|k + 1) T .

NUMERICAL EXAMPLES

Consider a hybrid descriptor stochastic system as formulated in Section 2, with (3) and ( 4), where r = 3, the mode dependent system matrices

A i , B i , C i , D i , E i , Q i , S i , Π
and the control u are as follows:

A 1 = 0 1 0 0.5 0 1 0 0 -0.1 , A 2 = 0.1 0 0.3 0 0.6 0 0 0.3 0.9 , A 3 = 0 1 0 0 0 1 -0.05 -0.1 1 , B i = 1 2 3 , D i = 0 0 ; C 1 = 0 0 1 1 1 0 , C 2 = 1 0 1 1 1 0 , C 3 = 0 1 1 1 1 0 ; E 1 = 1 0 0 0 1 1 0 0 0 , E 2 = 2 0 3 0 2 4 0 0 0 , E 3 = 3 0 0 0 4 0 0 0 0 ; Q i = 0.001 • 1 0 0 0 2 0 0 0 3 , S i = 10 • 4 0 0 1 , i = 1, 2, 3;
, u(k) =      1, if, 21 ≤ k ≤ 60; if, 101 ≤ k ≤ 140; if, 181 ≤ k ≤ 220; 0, otherwise.
We will denote x i , i = 1, 2, 3, y 1 , y 2 as the components of vector x and y, respectively in the above figures.

The simulation runs from k = 1 till k = 300. The IMM estimator for hybrid descriptor systems proposed in this paper is applied to the simulated system for state estimation. The Markov mode sequence can be seen in Fig. 1. System states and outputs are presented in Fig. 2, and the state estimates in Fig. 3. In Fig. 2, the results attained by the time varying descriptor system for Kalman filter [START_REF] Ali | An innovations approach to fault diagnosis in linear time-varying descriptor systems[END_REF]; [START_REF] Nikoukhah | Kalman filtering and riccati equations for descriptor systems[END_REF]] are also presented in the aim of comparison. This descriptor system Kalman filter assumes that the actual mode sequence is known, which is not the case of the IMM estimator.

The results shown in Fig. 3 are based on one random realization of the simulated stochastic system. In order to illustrate the statistic properties of the proposed IMM estimator, 1000 random realizations have been made, and the histograms of the state estimation errors at instants k = 50,140,200,are shown in Figs. 4,5,6.

CONCLUSIONS

In this paper the existing interacting multiple model estimator has been extended to the case of stochastic hybrid systems with modes described by descriptor equations. We have presented a numerical example to show the performance of our method. Among future research directions, joint estimation of state and parameters in descriptor hybrid systems will be studied, notably for the purpose of fault diagnosis.

Appendix A. GAIN MATRICES AND ESTIMATED STATE

In this appendix, we will check that the state estimation xd (k + 1|k + 1) computed in ( 7) is indeed the filtered state, not the predicted state, by showing that, in the particular case of E k+1 = I, it coincides with the filtered state of the classical Kalman filter. The explicit forms of the gain matrices G d (k + 1) and L d (k + 1) in ( 7) are also derived in this appendix. For simplicity of notation, we omit M i (k + 1), and set [L k+1 , K k+1 ] = L Mi(k+1) , K Mi(k+1) , etc, and B k = D k = 0. We are going to compute explicitly L k , and K k which were defined in [START_REF] Ali | An innovations approach to fault diagnosis in linear time-varying descriptor systems[END_REF]] as

[L k+1 , K k+1 ] = (H T k+1 R -1 k H k+1 ) -1 H T k+1 R -1 k , (A.1) where R k = Σ k O O S k and Σ k = A k P k A T k + Q k .
Then, we have

H T k+1 R -1 k H k+1 = E T k+1 , C T k+1 Σ -1 k O O S -1 k E k+1 C k+1 = E T k+1 Σ -1 k E k+1 + C T k+1 S -1 k C k+1 . Recall the matrix inverse formula, (A + BCD) -1 = A -1 -A -1 B •(DA -1 B + C -1 ) -1 DA -1
where A -1 , C -1 , and (DA -1 B + C -1 ) are assumed to exist.

Using the matrix inverse formula, we obtain

(H T k+1 R -1 k H k+1 ) -1 = (E T k+1 Σ -1 k E k+1 + C T k+1 S -1 k C k+1 ) -1 = (E T k+1 Σ -1 k E k+1 ) -1 -(E T k+1 Σ -1 k E k+1 ) -1 C T k+1 • C k+1 E T k+1 Σ -1 k E k+1 -1 C T k+1 + S k -1 •C k+1 (E T k+1 Σ -1 k E k+1 ) -1 . Let us define G k+1 = (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 •(C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 +S k ) -1 . (A.2) Further, (H T k+1 R -1 k H k+1 ) -1 H T k+1 R -1 k = (E T k+1 Σ -1 k E k+1 ) -1 -G k+1 C k+1 (E T k+1 Σ -1 k E k+1 ) -1 • E T k+1 , C T k+1 Σ -1 k O O S -1 k = (E T k+1 Σ -1 k E k+1 ) -1 E T k+1 Σ -1 k -G k+1 C k+1 (E T k+1 Σ -1 k E k+1 ) -1 E T k+1 Σ -1 k , (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 S -1 k -G k+1 C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 S -1 k .
Immediately, from (A.1), we get

L k+1 = (E T k+1 Σ -1 k E k+1 ) -1 E T k+1 Σ -1 k -G k+1 C k+1 (E T k+1 Σ -1 k E k+1 ) -1 E T k+1 Σ -1 k , and
K k+1 = (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 S -1 k -G k+1 C k+1 (E T k+1 Σ -1 k E k+1 ) -1 •C T k+1 S -1 k . (A.
3) Substituting (A.2) into (A.3), we have

K k+1 = (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 •(C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 +S k ) -1
, which is exactly equals to G k+1 expressed by (A.2). To obtain this result, the following equalities are used

C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 + S k -1 = S -1 k -(C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 + S k ) -1 •C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 S -1 k . (A.4) Indeed, (C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 + S k ) • (S -1 k -(C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 + S k ) -1 •C k+1 (E T k+1 Σ -1 k E k+1 ) -1 C T k+1 S -1 k = C k+1 E T k+1 Σ -1 k E k+1 -1 C T k+1 S -1 k + I -C k+1 E T k+1 Σ -1 k E k+1 -1 C T k+1 S -1 k = I.
Hence, in particular, if we set E k+1 = I, one can get

L k+1 = I -G k+1 C k+1 , G k+1 = Σ k C T k+1 (C k+1 Σ k C T k+1 + S k ) -1
, which is just the Kalman gain. At last, we have

x(k + 1) = L k+1 A k x(k) + G k+1 y(k + 1) = A k x(k) + G k+1 [y(k + 1) -C k+1 A k+1 x(k)].
The associated variance can be verified similarly.

  Fig. 1. Markov mode sequence.

Fig. 4 .

 4 Fig. 4. Histogram for k = 50, 140, 200 of the state estimation error of x 1 , respectively.

Fig. 5 .Fig. 6 .Fig. 3 .

 563 Fig. 5. Histogram for k = 50, 140, 200 of the state estimation error of x 2 , respectively.


Boyi Ni

Qinghua Zhang
email: qinghua.zhang@inria.fr

Stability of the Kalman Filter for Output Error Systems

Keywords: Stability, Kalman filter, Output Error System

Optimality and numerical efficiency are well known properties of the Kalman filter, whereas its stability property, though equally classical and important in practice, is less often mentioned in the recent literature. The stability of the Kalman filter is usually ensured by the uniform complete controllability regarding the process noise and the uniform complete observability of linear time varying systems. Such classical results cannot be applied to output error systems, in which the process noise is totally absent. It is shown in this paper that the uniform complete observability is sufficient to ensure the stability of the Kalman filter applied to time varying output error systems, regardless of the stability of the considered system itself.

INTRODUCTION

The well known Kalman filter has been extensively studied and is being applied in many different fields [START_REF] Anderson | Optimal Filtering[END_REF]; [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]; [START_REF] Zarchan | Fundamentals of Kalman filtering: a practical approach[END_REF]; [START_REF] Grewal | Kalman Filtering: Theory and Practice Using Matlab[END_REF]). The purpose of the present paper is to study the stability of the Kalman filter in a particular case rarely covered in the literature: the absence of process noise in the state equation of a linear time varying (LTV) system. Such systems are known as output error (OE) systems. Though typically process noise and output noise are both considered in Kalman filter applications, the case with no process noise is of particular interest when state equations come from physical laws that are believed accurate enough. It is also important for the application of the Kalman filter to OE system identification [START_REF] Goodman | An output error model and algorithm for electromegnetic system identification[END_REF]; [START_REF] Forssell | Identification of unstable systems using output error and box-jenkins model structures[END_REF]).

While the optimal property of the Kalman filter is frequently recalled, its stability property is less often mentioned in the recent literature. The classical stability analysis is based both on the uniform complete controllability regarding the process noise and on the uniform complete observability of the considered system [START_REF] Kalman | New methods in Wiener filtering theory[END_REF]; [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]). In the case of OE systems, there is no process noise at all in the state equation, hence the controllability regarding the process noise cannot be fulfilled, and the classical stability results are not applicable. The present paper aims at completing this missing case.

The optimal state estimation realized by the Kalman filter is usually viewed as a trade-off between the uncertainties in the state equation and in the output equation. In an OE system, the state equation is assumed noise-free. This point of view suggests that the state estimation should This work has been partly supported by the ITEA2 MODRIO project.

solely rely on the state equation, provided that the initial state of the OE system is exactly known. In practice the Kalman filter remains useful when the initial state is not exactly known or when the OE system is unstable. Of course, if the state of an unstable system diverges, so does its state estimate by the Kalman filter. Typically in practice, unstable systems are stabilized by feedback controllers so that the system state remains bounded. The Kalman filter can be applied either to the controlled system itself or to the entire closed loop system. In the latter case, the controller must be linear and completely known, excluding the saturation protection and any other nonlinearities.

The classical optimality results of the Kalman filter are also valid in the case of OE systems (Jazwinski, 1970, chapter 7). Nevertheless, it remains to complete the stability analysis, as the classical results are not applicable here.

The main results presented in this paper are as follows. Under the uniform complete observability condition, the dynamics of the Kalman filter applied to a LTV OE system is asymptotically stable, regardless of the stability of the system itself. The boundedness of the solution of the Riccati equation, which ensures the boundedness of the Kalman gain, is also proved under the same condition. These results are quite similar to the classical results [START_REF] Kalman | New methods in Wiener filtering theory[END_REF]; [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]), which exclude the case of OE systems.

For linear time invariant (LTI) systems, it is a common practice to design the Kalman filter by solving an algebraic Riccati equation (in contrast to dynamic differential Riccati equation for general LTV systems as considered in the present paper). In this case, the controllability and observability conditions can be replaced by the weaker stabilizability and detectability conditions [START_REF] Laub | A Schur method for solving algebraic Riccati equations[END_REF]; [START_REF] Arnold | Generalized eigenproblem algorithms and software for algebraic Riccati equations[END_REF]). Some preliminary results about LTI OE systems have been presented in [START_REF] Ni | Kalman filter for process noise free systems[END_REF]).

The rest of the paper is organized as follows. Some preliminary elements are introduced in Section 2. The problem considered in this paper is formulated in Section 3. The properties of the solution of the Riccati equation are analyzed in Section 4. The stability of the Kalman filter for OE systems is established in Section 5. Some numerical examples are presented in Section 6. Finally, concluding remarks are drawn in Section 7.

DEFINITIONS

Let us shortly recall some definitions about LTV systems, which are necessary for the following sections.

Let m and n be any two positive integers. For a vector x ∈ R n , x denotes its Euclidean norm. For a matrix A ∈ R m×n , A denotes the matrix norm induced by the Euclidean vector norm, which is equal to the largest singular value of A and known as the spectral norm when m = n. Then Ax ≤ A x for all A ∈ R m×n and all x ∈ R n . For two real square symmetric positive definite matrices A and B, A > B means A-B is positive definite.

Let A(t) ∈ R m×n be defined for t ∈ R. (3) 2

The last concept to be recalled here is about the observability of LTV systems, following [START_REF] Kalman | New methods in Wiener filtering theory[END_REF]). Definition 3. The matrix pair [A(t), C(t)] with A(t) ∈ R n×n and C(t) ∈ R m×n is uniformly completely observable if there exist positive constants τ , ρ 1 and ρ 2 such that, for all t ∈ R, the following inequalities hold

ρ 1 I n ≤ t t-τ Φ T (s, t)C T (s)R -1 (s)C(s)Φ(s, t)ds (4)
≤ ρ 2 I n (5) with some bounded symmetric positive definite matrix R(s) ∈ R m×m (typically the covariance matrix of the output noise in a stochastic state space system). 2

PROBLEM FORMULATION AND ASSUMPTIONS

In this section is first recalled the classical Kalman filter under its usual assumptions, before the presentation of the particular case of OE systems considered in this paper.

Kalman filter in the usual case

The Kalman filter in continuous-time is usually applied to LTV systems modeled by

dx(t) = A(t)x(t)dt + B(t)u(t)dt + Q 1 2 (t)dω(t) (6a) dy(t) = C(t)x(t)dt + R 1 2 (t)dη(t) (6b) where t ∈ R represents the time, x(t) ∈ R n is the state vector, u(t) ∈ R l the bounded input, y(t) ∈ R m the output, ω(t) ∈ R n , η(t) ∈ R m are two independent
Brownian processes with identity covariance matrices, A(t), B(t), C(t), Q(t), R(t) are real matrices of appropriate sizes. The matrix Q(t) is symmetric positive semi-definite, and R(t) is symmetric positive definite. The notations Q 1 2 (t) and R 1 2 (t) denote respectively the symmetric positive (semi)-definite matrix square roots of Q(t) and R(t). The initial state x(t 0 ) ∈ R n is a random vector following the Gaussian distribution x(t 0 ) ∼ N (x 0 , P 0 ) with x 0 ∈ R n and P 0 ∈ R n×n .

The Kalman filter for this LTV system writes

dx(t) = A(t)x(t)dt + B(t)u(t)dt + K(t)(dy(t) -C(t)x(t)dt) (7a)
K(t) = P (t)C T (t)R -1 (t) (7b) d dt P (t) = A(t)P (t) + P (t)A T (t) -P (t)C(t) T R -1 (t)C(t)P (t) + Q(t) (7c) x(t 0 ) = x 0 , P (t 0 ) = P 0 (7d)
where the solution of the Riccati equation ( 7c) is a matrix function P (t) ∈ R n×n and the Kalman gain

K(t) ∈ R n×m .
The optimal properties of the Kalman filter are well known. It is less well known, though equally classical and important, that the solution P (t) of the Riccati equation is bounded and that the dynamics of the Kalman filter is stable, provided the matrix pair [A(t), Q 1 2 (t)] is uniformly completely controllable and the matrix pair [A(t), C(t)] is uniformly completely observable [START_REF] Kalman | New methods in Wiener filtering theory[END_REF]; [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]). These are obviously crucial properties in practice for online applications.

It is worth mentioning that, in [START_REF] Kalman | New methods in Wiener filtering theory[END_REF]; [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]), these classical results are based on an important lemma, which turns out to be incorrect, as pointed out independently by the authors of [START_REF] Delyon | A note on uniform observability[END_REF]; [START_REF] Pengov | On the boundedness of the solutions of the continuous riccati equation[END_REF]). Fortunately, the mistake has been repaired in these more recent references so that the main classical results remain correct.

Output error systems and Kalman filter

In the case of OE systems, the process noise is absent from the state equation, hence the general LTV system (6) becomes

dx(t) = A(t)x(t)dt + B(t)u(t)dt (8a) dy(t) = C(t)x(t)dt + R 1 2 (t)dη(t). ( 8b
)
An OE system can be seen as a particular LTV system (6) with Q(t) ≡ 0. The Kalman filter (7) then becomes

dx(t) = A(t)x(t)dt + B(t)u(t)dt + K(t)(dy(t) -C(t)x(t)dt) (9a) K(t) = P (t)C T (t)R -1 (t) (9b) d dt P (t) = A(t)P (t) + P (t)A T (t) -P (t)C(t) T R -1 (t)C(t)P (t) (9c) x(t 0 ) = x 0 , P (t 0 ) = P 0 (9d)
In this case, the uniform complete controllability condition cannot be satisfied by the matrix pair [A(t), Q 1 2 (t)], as Q 1 2 (t) ≡ 0. Consequently, the classical results on the stability of the Kalman filter cannot be applied here. The present paper is for the purpose of studying the Kalman filter stability in this particular case.

Assumptions

The assumptions stated here are required throughout this paper.

The considered OE system (8) is defined with bounded and piecewise continuous real matrices A(t), B(t), C(t), R(t) of appropriate sizes, among which R(t) is positive definite. It is also assumed that R -1 (t) is bounded.

The initial state x(t 0 ) ∈ R n is a random vector following the Gaussian distribution

x(t 0 ) ∼ N (x 0 , P 0 ). ( 10) with some known mean vector x 0 ∈ R n and symmetric positive definite covariance matrix P 0 ∈ R n×n .

It is further assumed that the matrix pair [A(t), C(t)] is uniformly completely observable (see Definition 3).

PROPERTIES OF THE SOLUTION OF THE RICCATI EQUATION

Let us first establish the positive definiteness and the upper boundedness of P (t), which are of crucial importance in practice, before studying the stability of the Kalman filter.

It will be shown that the solution of the Riccati equation (9c) is closely related to the solution of the Lyapunov equation

dΩ(t) dt + A T (t)Ω(t) + A(t)Ω(t) = C T (t)R -1 (t)C(t) (11a)
Ω(t 0 ) = P -1 0 (11b) with the same matrices A(t), C(t), R(t), P 0 as in (9c) and (9d).

Proposition 1. The solution Ω(t) ∈ R n×n of the Lyapunov equation ( 11) is a symmetric positive definite matrix for all t ≥ t 0 .

Proof. It can be directly checked that Ω(t) = Φ T (t 0 , t)Ω(t 0 )Φ(t 0 , t)

+ t t0 Φ T (s, t)C T (s)R -1 (s)C(s)Φ(s, t)ds (12)
satisfies the Lyapunov equation ( 11). As Ω(t 0 ) = P -1 0 is positive definite, Φ(t 0 , t) is an invertible matrix, and for all t ≥ t 0 the integral in ( 12) is positive semidefinite, therefore, Ω(t) is positive definite for all t ≥ t 0 .

2

Now it is certain that, for all t ≥ t 0 , Ω(t) is positive definite, thus invertible. By applying the matrix differentiation rule

dΩ -1 (t) dt = -Ω -1 (t) dΩ(t) dt Ω -1 (t),
the following result can be directly checked.

Proposition 2. Let Ω(t) be the solution the Lyapunov equation ( 11) for t ≥ t 0 , then P (t) Ω -1 (t) solves the Riccati equation (9c) with the initial condition P (t 0 ) = P 0 . 2

It is then clear that the properties of P (t) can be studied through those of Ω(t).

Proposition 3. Under the assumptions stated in Section 3.3, for all t ≥ t 0 , the solution P (t) of the Riccati equation (9c) with the initial condition P (t 0 ) = P 0 > 0 is symmetric positive definite and is upper bounded.

2

Proof. The positive definiteness of P (t) is immediate from that of Ω(t). To show that P (t) is upper bounded, it will be shown that Ω(t) is lower bounded.

First consider the case t ∈ [t 0 , t 0 + τ ] (remind that τ is the positive constant involved in the observability condition ( 4)). In this finite time interval, the first term in ( 12) is positive definite, and the second term is positive semidefinite, then Ω(t) ≥ ρ 0 I n (13) where ρ 0 > 0 is the minimum of the smallest singular value of the matrix Φ T (t 0 , t)Ω(t 0 )Φ(t 0 , t) for t ∈ [t 0 , t 0 + τ ].

It remains to consider the case t > t

0 + τ . Let N t -t 0 τ (14)
be the largest integer smaller than or equal to (t -t 0 )/τ . For the second term of (12), decompose the integral into the sum of N integrals over intervals of size τ , possibly dropping the part of the integral at the beginning of [t 0 , t 0 + τ ] smaller than τ if (t -t 0 )/τ is not exactly an integer. It then yields

t t0 Φ T (s, t)C T (s)R -1 (s)C(s)Φ(s, t)ds ≥ N -1 k=0 t-kτ t-(k+1)τ Φ T (s, t)C T (s)R -1 (s)C(s)Φ(s, t)ds = N -1 k=0 Φ T (t -kτ, t)S k Φ(t -kτ, t) (15)
with

S k t-kτ t-(k+1)τ Φ T (s, t -kτ )C T (s)R -1 (s)C(s)
• Φ(s, t -kτ )ds ( 16) Each S k (for k = 0, 1, . . . , N-1) corresponds to the integral in the uniform complete observability condition (4), which holds for all t ∈ R. Hence S k ≥ ρ 1 I n for k = 0, 1, . . . , N-1, with the positive constant ρ 1 as in (4). Then

t t0 Φ T (s, t)C T (s)R -1 (s)C(s)Φ(s, t)ds ≥ ρ 1 N -1 k=0 Φ T (t -kτ, t)Φ(t -kτ, t) (17)
Now in the sum of ( 17) keep only the term with k = 0 (the other terms will be useful for proving other results). Then

t t0 Φ T (s, t)C T (s)R -1 (s)C(s)Φ(s, t)ds (18) ≥ ρ 1 Φ T (t, t)Φ(t, t) (19) ≥ ρ 1 I n
(20) It means that the second term of ( 12) is lower bounded by ρ 1 I n for t > t 0 + τ . By summarizing the results for the cases t ∈ [t 0 , t 0 + τ ] and t > t 0 + τ , it is thus shown that Ω(t) is lower bounded by a strictly positive definite matrix for all t ≥ t 0 . To be more specific, Ω(t) ≥ min(ρ 0 , ρ 1 )I n (21)

It is then concluded that P (t) = Ω -1 (t) is upper bounded for all t ≥ t 0 . 2

STABILITY OF THE KALMAN FILTER FOR OUTPUT ERROR SYSTEMS

For the Kalman filter ( 9), it is already shown that the matrix P (t) is upper bounded, it then follows from the assumptions about the boundedness of C(t) and R -1 (t) that the Kalman gain K(t) is also upper bounded. Equation (9a) governing x(t) can be viewed as a dynamic system driven by the "exogenous input" terms B(t)u(t)dt and K(t)dy, but its intrinsic stability property is only related to its homogeneous part, namely the LTV system

d dt z(t) = (A(t) -K(t)C(t))z(t) (22)
with z(t) ∈ R n . Like in [START_REF] Jazwinski | Stochastic Processes and Filtering Theory[END_REF]), when the stability of the Kalman filter is talked about in this paper, it is about the stability of this homogeneous LTV system.

Before studying the stability of system ( 22), let us recall a classical result. Proofs of this lemma can be found in (Anderson et al., 1986, page 38, Lemma 2.3), (Ioannou and Sun, 1996, page 221, Lemma 4.8.1) and [START_REF] Zhang | Observability conservation by output feedback and observability Gramian bounds[END_REF]).

If the state equation ( 22) was associated with a deterministic observation equation y(t) = C(t)z(t), the term K(t)C(t)z(t) in equation ( 22) would be equal to K(t)y(t) and could be seen as an output feedback. Lemma 1 means that such an output feedback preserves the observability of an LTV system. Now it is ready to present the main result of this section.

Theorem 1. Under the assumptions stated in Section 3.3, the homogeneous part of the Kalman filter state estimation equation (9a), as expressed in equation ( 22), is asymptotically stable. 2

As mentioned at the beginning of this section, the stability of the homogeneous part of the Kalman filter is an intrinsic stability property independent of the signals being processed by the filter. It implies that the mathematical expectation of the state estimation error is asymptotically stable. Moreover, the boundedness of the covariance matrix of the state estimation error, namely P (t), is ensued by Proposition 3.

Proof of Theorem 1.

Define the Lyapunov function candidate V (z(t), t) z T (t)Ω(t)z(t) (23) with the positive definite Ω(t) as defined in (11).

Compute the derivative of V (z(t), t) along the trajectory of ( 22),

dV (z(t), t) dt = z T (t) (A(t) -K(t)C(t)) T Ω(t) + Ω(t)(A(t) -K(t)C(t)) + dΩ(t) dt z(t) (24)
Remind that Ω(t) satisfies (11a), then

dV (z(t), t) dt = z T (t) (-K(t)C(t)) T Ω(t) + Ω(t)(-K(t)C(t)) + C T (t)R -1 (t)C(t) z(t)
In (9b) replace P (t) by Ω -1 (t), then

K(t) = Ω -1 (t)C T (t)R -1 (t), hence dV (z(t), t) dt = -z T (t)C T (t)R -1 (t)C(t)z(t) (25)
which is negative semidefinite. It is then clear that the value of V (z(t), t) cannot increase with the time t ≥ t 0 , thus

V (z(t), t) ≤ V (z(t 0 ), t 0 ), ∀t ≥ t 0 . (26) It then follows from the definition of V (z(t), t) that z T (t)Ω(t)z(t) ≤ z T (t 0 )Ω(t 0 )z(t 0 )
Let σ 0 be the smallest singular value of P 0 , then σ -1 0 is the largest singular value of Ω(t 0 ) = P -1 0 . Remind the lower bound of Ω(t) for all t ≥ t 0 as shown in (21), then

min(ρ 0 , ρ 1 ) z(t) 2 ≤ z T (t)Ω(t)z(t) ≤ z T (t 0 )Ω(t 0 )z(t 0 ) ≤ σ -1 0 z(t 0 ) 2 Hence z(t) ≤ 1 σ 0 min(ρ 0 , ρ 1 ) z(t 0 ) (27)
holds for all z(t 0 ) ∈ R n . Therefore the homogeneous system ( 22) is Lyapunov stable.

To show the asymptotic stability of ( 22), the limiting behavior of z(t) when t → +∞ will be analyzed.

Let the state transition matrix of the homogeneous LTV system (22) be denoted by Φ K (s, t), then z(s) = Φ K (s, t)z(t) for any s, t ∈ R. It then follows from (25) that dV (z(s), s) ds (31) This result, together with (30), leads to

= -z T (t)Φ T K (s, t)C T (s)R -1 (s)C(s)Φ K (s, t)z(t) and therefore V (z(t), t) -V (z(t -τ ), t -τ ) (28) = t t-τ dV (z(s), s) ds ds (29) = -z T (t)O K (t, t -τ )z(t) (30) with O K (t, t -τ ) t t-τ Φ T K (s, t)C T (s)R -1 (s)C(s)Φ K (s,
V (z(t), t) ≤ V (z(t -τ ), t -τ ) -ρ 3 z(t) 2 .
(32) It means that, over each time interval of size τ , the value of V (z(t), t) is decreased by ρ 3 z(t) 2 ≥ 0. Though V (z(t), t) is decreased over each time interval of size τ , it may not necessarily tend to zero, as the decrement may be infinitesimal.

It will be shown that z(t) tends to zero by means of contradiction.

For the purpose of proof by contradiction, assume that z(t) does not tend to zero when t → +∞. Then there exists a constant ε > 0, such that for any (arbitrarily large) T ∈ R, there exists t > T such that z(t) > ε, therefore

V (z(t), t) ≤ V (z(t -τ ), t -τ ) -ρ 3 ε 2 .
(33) There exist infinitely many such values of t, as T can be arbitrarily large. Moreover, it was already shown that the value of V (z(t), t) cannot increase with the time t. The inequality (33) then says that V (z(t), t) is repeatedly decreased by ρ 3 ε 2 for larger and larger values of t. Consequently, V (z(t), t) will become negative for sufficiently large values of t. This is in contradiction with the definition of V (z(t), t) in ( 23) which is positive definite. Therefore, it is proved that z(t) tends to zero when t → +∞.

The asymptotic stability of the homogeneous system ( 22) is then established. The examples of LTV systems (the A(t) and Φ(t, t 0 ) matrices) have been inspired by (Choi, 2010, Lecture 24). The numerical solutions are computed with the ode45 solver in Matlab.

λ 1 (t) λ 2 (t) z 1 (t) z 2 (t)
Fig. 1. Example 1. Top: eigenvalues of P (t) solution of (9c), bottom: components of z(t) solution of ( 22).

Example 1 -Lyapunov stable system.

Consider

A(t) = cos(0.2t) sin(0.2t) -sin(0.2t) cos(0.2t) , C(t) = 1.5 0 0 2 , R(t) = I 2 .
The corresponding state transition matrix is

Φ(t, t 0 ) = exp(ϕ(t, t 0 )) cos(ψ(t, t 0 )) sin(ψ(t, t 0 )) -sin(ψ(t, t 0 )) cos(ψ(t, t 0 ))
with ϕ(t, t 0 ) 5(sin(0.2t) -sin(0.2t 0 )) ψ(t, t 0 ) 5(cos(0.2t 0 ) -cos(0.2t)).

This system is Lyapunov stable (see Definition 1).

The Riccati equation (9c) is initialized with P (0) = I 2 and the homogenous system (22) with z(0) = [1, 2] T . In Figure 1 are shown the two eigenvalues (equal to the singular values) of P (t), which are bounded and tend to zero. In Figure 1 is also shown the state vector z(t) (solution of ( 22)), which converges to zero.

Example 2 -unstable system.

Consider

A(t) = 2-e -t e -t -2 0 2-e -t , C(t) = [1 0] , R(t) = 1.
The corresponding state transition matrix is

Φ(t, t 0 ) = exp(ψ(t, t 0 )) -ψ(t, t 0 ) exp(ψ(t, t 0 )) 0 exp(ψ(t, t 0 ))
with ψ(t, t 0 ) 2(t-t 0 )+e -t -e -t0 . This system is unstable (Φ(t, t 0 ) diverges when t → +∞).

The Riccati equation ( 9c) is initialized with P (0) = I 2 and the homogenous system (22) with z(0) = [1, 2] T . In Figure 2 are shown the two eigenvalues (equal to the singular values) of P (t), which are both upper and lower bounded. In Figure 2 is also shown the state vector z(t) (solution of ( 22)), which converges to zero.

CONCLUSION

For any recursive algorithm running continuously in real time, the boundedness of all the involved variables is obviously an important property. It is established in this paper that, when the Kalman filter is applied to LTV OE systems, the solution of the Riccati equation and the Kalman gain are both bounded, essentially under the observability condition. It is further shown that the (homogeneous part of the) Kalman filter is asymptotically stable. The case of OE systems studied in this paper is not covered by the classical results requiring a controllability condition regarding the process noise, which is totally absent in OE systems. This difference caused a technical difficulty: unlike in the classical case, here the solution of the Riccati equation does not always have a strictly positive definite lower bound, therefore the associated "natural" Lyapunov function cannot be used in the usual sense in the convergence proof. The results of this paper show that, when the Kalman filter is applied to an OE system, there is no need to artificially introduce a small process noise into the state equation that would alter the results of state estimation.

The asymptotic stability of the Kalman filter for LTV OE systems has been established in this paper, regardless of the stability of the considered systems. This result can be further refined, by characterizing the exponential or polynomial convergence rate of the Kalman filter, but such developments depend on the stability or instability properties of the considered systems. Due to space limitation of the present paper, such results will be reported elsewhere.

  This section presents some numerical examples confirming the theoretical results established in the previous sections.

Fig. 2 .

 2 Fig. 2. Example 2. Top: eigenvalues of P (t) solution of (9c), bottom: components of z(t) solution of (22).