M Rodrigues
email: mickael.rodrigues@univ-lyon1.fr

H Hamdi

D Theilliol
email: didier.theilliol@univ-lorraine.fr

C Mechmeche

N Benhadj Braiek

Actuator Faut Estimation based Adaptive Polytopic Observer for a Class of LPV Descriptor Systems

Keywords: Adaptive Polytopic Observer, LPV descriptor systems, Actuator Fault Estimation, LMI

This paper presents an Adaptive Polytopic Observer (APO) design in order to develop an actuator fault estimation method dedicated to polytopic Linear Parameter Varying (LPV) descriptor systems. This paper extends a Fault Diagnosis (FD) method developed for regular LTI systems to polytopic LPV descriptor systems. Here, time-varying actuator faults are also considered whereas in many papers actuator faults are generally assumed to be constant. The design and convergence conditions of this APO are provided. The design is formulated through Linear Matrix Inequalities (LMI) techniques under equality constraints. The performances of the proposed actuator fault estimation scheme are illustrated using an electrical circuit.

INTRODUCTION

During the last three decades, both theoretical and experimental researches on fault detection have been intensively developed. For obvious economical reasons, systems in the real world have to work perfectly at all time under all conditions. That's why it is crucial to be able to detect and identify the possible faults that may affect a system as early as possible in order to prevent significant performance degradations or damages of the system [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF], [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF]. Fault Detection and Diagnosis (FDD) has been well studied for physical plants which were modeled by T-S fuzzy models. Recently Such efforts studded by Jiang B. et al. [START_REF] Jiang | Integrated fault estimation and accommodation design for discrete-time Takagi-Sugeno fuzzy systems with actuator faults[END_REF], have led to the development of a method for integrated robust fault estimation and accommodation for a class of discrete-time nonlinear systems described by a T-S fuzzy model. Also, dynamic output feedback fault-tolerant controllers have been developed by Zhang K. et al. [START_REF] Zhang | Dynamic output feedback fault tolerant controller design for Takagi-Sugeno fuzzy systems with actuator faults[END_REF] for T-S fuzzy systems with actuator faults in which a fuzzy augmented fault observer is proposed to yield fault estimates and, based on the information of on-line fault estimates, observer based output feedback fault-tolerant controllers are designed. Nevertheless for nonlinear systems it still remain a challenge [START_REF] Puig | Passive Robust Fault Detection of Dynamic Processes Using Interval Models[END_REF].

By the way, some techniques have been developed to approximate nonlinear systems by Linear Parameter Varying (LPV) models that can be used efficiently to represent nonlinear systems as in [START_REF] Bamieh | Identification of linear parameter varying models[END_REF], [START_REF] Rodrigues | Sensor Fault Detection and Isolation Filter for Polytopic LPV Systems: A Winding Machine Application[END_REF] for a winding machine, [START_REF] Rotondo | Quasi-LPV modeling, identification and control of a twin rotor MIMO system[END_REF] for a twin rotor system. These methods allow to apply powerful linear design tools to complex nonlinear models. Based on this LPV representation of nonlinear systems, some researchers from the FDI community have developed model-based methods using LPV ordinary models as in [START_REF] Armeni | Robust Fault detection and isolation for LPV systems under a sensitivity constraint[END_REF], [START_REF] Bokor | Detection filter design for LPV systems-a geometric approach[END_REF], and [START_REF] Lira | LPV observer design for PEM fuel cell system: Application to fault detection[END_REF]. Similar result is found in [START_REF] Wei | LMI solutions to the mixed H -/H∞ fault detection observer design for linear parametervarying systems[END_REF] for linear parameter-varying (LPV) systems the presented method is based on multi-objective H -/H ∞ fault detection observer design. LPV sliding mode observers have been largely introduced by [START_REF] Alwi | Robust Fault Reconstruction for Linear Parameter Varying Systems Using Sliding Mode Observers[END_REF] so as to reconstruct the states of the system. A particular class of LPV systems is the polytopic LPV form which allows to describe the system as a convex combination of sub-models defined by the vertices of a convex polytope [START_REF] Grenaille | A method for designing Fault Diagnosis Filters for LPV polytopic systems[END_REF], [START_REF] Hamdi | Robust Fault Detection and Estimation for Descriptor Systems Based on Multi-Models Concept[END_REF], [START_REF] De | Fault estimation and virtual sensor FTC approach for LPV systems[END_REF]. These sub-models are then combined by convex weighting functions that yield to a global model.

Generally speaking, most of control and fault diagnosis methods for physical systems, use normal (or regular) models i.e. there is no algebraic relations between the system variables. However, Differential-Algebraic Equations (DAE) or implicit systems or singular systems or descriptor systems are of quite importance for the physical representation of some systems [START_REF] Lewis | A Survey of Linear Singular Systems Circuits Systems and Signal Processing[END_REF], [START_REF] Darouach | Design of Observers for Descriptor Systems[END_REF]. Such systems appear for example in electrical circuits, mechanical systems with holonomic or non holonomic constraints, robotic systems with kinematical constraints and chemical systems [START_REF] Mattson | Physical system modeling with modelica[END_REF]. Some practical problems must take into account physical constraints or algebraic relations and more generally impulsive behaviors caused by an improper transfer matrix: see the following books on singular systems [START_REF] Dai | Singular Control Systems[END_REF], [START_REF] Duan | Analysis and Design of descriptor Linear Systems Advances in Mechanics and Mathematics[END_REF].

Various observer-based fault diagnosis methods for descriptor systems have been proposed but generally for Linear Time Invariant models [START_REF] Wang | Observer design for discrete-time descriptor systems: An LMI approach[END_REF]. The author in [START_REF] Koenig | Unknown Input Proportional Multiple-Integral Observer Design for Linear Descriptor Systems: Application to State and Fault Estimation[END_REF] studies an observer coming from a continuous nonlinear descriptor systems via a convex optimization. Recently, a fault diagnosis method to detect and estimate actuator faults for multi-models descriptor systems represented by using unknown input observers [START_REF] Hamdi | Robust Fault Detection and Estimation for Descriptor Systems Based on Multi-Models Concept[END_REF]. Therefore, in [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor Linear Parameter Varying systems[END_REF] the authors have presented a method of fault estimation for a particular class of discrete-time LPV descriptor systems. The authors in [START_REF] Hamdi | State Estimation for Polytopic LPV Descriptor Systems: Application to Fault Diagnosis[END_REF] and [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF], have proposed a polytopic unknown inputs and proportional integral observers for LPV descriptor systems respectively. These observers are used to detect and isolate actuator faults. However, all these observers can only estimate constant faults which is a restrictive condition.

In order to estimate both constant and time varying faults, an Adaptive Polytopic Observer is proposed in this paper. In many applications where the parameters are unknown and states are not accessible, adaptive observer appears to be a valuable method in estimation of both parameters and states of the system. In [START_REF] Wang | Actuator fault diagnosis: An adaptive observer-based technique[END_REF] the authors have presented an adaptive fault diagnosis observer approach dedicated to regular LTI systems which can detect and estimate constant faults. In [START_REF] Zhang | Adaptive Observer-based Fast Fault Estimation[END_REF], the authors have designed an adaptive observer so as to estimate time-varying faults but only for regular LTI systems.

Thus, this paper proposes an Adaptive Polytopic Observer (APO) for a class of LPV descriptor systems. Fault Diagnosis for descriptor LPV system still remains a challenge and as far as authors knows there are very few papers for such systems that consider time-varying actuator fault estimation under disturbances. In literature, most of fault diagnosis approaches are designed only for regular systems and not for descriptor systems. To tackle such problem, the entire development of this paper is a real contribution for descriptor LPV systems with significant new results in terms of time-varying fault estimation in spite of disturbances. Both methods presented in [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF] or in [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor Linear Parameter Varying systems[END_REF] can only estimate constant fault for descriptor LPV Systems; the real novelty in the paper consists in the ability of the approach to estimate time-varying actuator fault despite disturbances which is a real huge problem and a key problem to solve for some potential extensions to Fault Tolerant Control.

In comparison with the polytopic LPV proportional integral observer studied in [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF] which can only estimate constant actuator fault, the main goals of this paper are: to extend existing results about fault estimation of LTI regular systems to polytopic LPV descriptor systems, and also to consider time-varying actuator faults estimation. Sufficient convergence conditions for the APO are given by solving a set of Linear matrix inequalities (LMI).

The outline of this paper is as follows. The structure of the LPV descriptor systems is formulated in Section 2. The Adaptive Polytopic Observer-based actuator fault estimation is presented in Section 3. Finally, in Section 4, a numerical example that considers a LPV descriptor systems is used to assess the validity of the proposed approach.

In this paper, the notations are standard. R denotes the set of real numbers, C is the complex plane, A + denotes the generalized inverse of A; Q > 0 or (Q < 0) indicates the symmetric matrix Q is positive (or negative) definite; . denotes the standard norm symbol. Also, ∀ means "for all". An asterisk * denotes the transposed element in the symmetric position.

PROBLEM STATEMENT

Let us consider the following LPV descriptor system in fault-free case subject to disturbance:

E ẋ(t) = A(θ(t))x(t) + B(θ(t))(u(t) + f (t)) + R(θ(t))d(t) y(t) = Cx(t) + Hd(t) (1)
where x(t) ∈ R n is the state vector, u(t) ∈ R p is the inputs vector, d(t) ∈ R q is a bounded unknown input vector i.e., d(t) ≤ β which represents disturbances and y(t) ∈ R m represents the measured outputs vector. f (t) ∈ R p represents the actuator fault vector by an additive external signal as in [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF], [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF]. Sometimes in FDI, these malfunctions can also be represented by a multiplicative external signal as in [START_REF] Zhang | Design of Integrated Fault Detection, Diagnosis and Reconfigurable Control Systems[END_REF], [START_REF] Rodrigues | Fault Tolerant Control Design for Polytopic LPV Systems[END_REF] by the following faulty control input u f (t) = (I pγ)u(t) which is equivalent to an additive signal such that:

u(t) + f (t) where f (t) = -γu(t) with γ diag[γ 1 , γ 2 , . . . , γ p ], 0 ≤ γ k ≤ 1 such that γ k = 1 → a total failure of the k th actuator k ∈ [1, . . . , p] γ k = 0 → the k th healthy actuator (2)
where the term γ k represents the loss of effectiveness of k th actuator. Disturbances and faults are supposed to be unknown. Matrix E ∈ R n×n may be singular and rank(E) = r < n. A(•), B(•) and R(•) are continuous functions which depend affinely on the time-varying parameter vector θ(t) ∈ R l . Matrices C and H are fixed. It is also assumed that each component θ i (t), i ∈ [1, . . . , l] of this time-varying parameter vector θ(t) is bounded and lies into a hypercube such that [START_REF] Wu | Control of Parameter Varying Systems[END_REF]:

θ(t) ∈ Υ = {θ i | θ i (t) ≤ θ i (t) ≤ θ i (t)}, ∀i ∈ [1, . . . , l], ∀t ≥ 0 (3)
The matrices A(θ(t)), B(θ(t)), R(θ(t)) of the LPV descriptor system (1) with the affine parameter dependence (3) are represented such that ∀θ(t) ∈ Υ:

M (θ(t)) = M 0 + l i=1 θ i (t) M i (4)
where M stands for matrices A, B and R. The system (1) can be transformed into a convex interpolation of the vertices of Υ where the vertices S i of the polytope are defined such that [START_REF] Rodrigues | Design of an Active Fault Tolerant Control for Nonlinear Systems described by a Multi-Model Representation[END_REF]: . . , h] where h = 2 l . The polytopic coordinates are denoted ρ(θ(t)) and vary into the convex set Ω:

S i = A i , B i , R i , C, H , ∀i ∈ [1, .
Ω = ρ(θ(t)) ∈ R h , ρ(θ(t)) = [ρ 1 (θ(t)), ..., ρ h (θ(t))] T , ρ i (θ(t)) ≥ 0, ∀ i, h i=1 ρ i (θ(t)) = 1 (5)
Then, the usual assumptions related to LPV systems are: ρ(θ(t)) is bounded, it is assumed to be fault-free and it is available. To ease the presentation, it is assumed that the matrices A(•), B(•) and R(•) are given by convex combinations ∀t ≥ 0 and the polytopic LPV descriptor system with the time-varying parameter vector ρ(θ(t)) ∈ Ω given by: where

   E ẋ(t) = h i=1 ρ i (θ(t))(A i x(t)+B i u(t)+F i f (t)+ R i d(t) ) y(t) = Cx(t) + Hd(t) (6
A i ∈ R n×n , B i = F i ∈ R n×p , R i ∈ R n×q , C ∈ R m×n
and H ∈ R m×q are time invariant matrices defined for the i th model. Before starting the design of the APO, we assume that:

A1. rank(CF i ) = rank(F i ) = p ∀i = [1, . . . , h]
A2. The triple matrix (E, A i , C) is observable, i.e. [START_REF] Dai | Singular Control Systems[END_REF], [START_REF] Duan | Analysis and Design of descriptor Linear Systems Advances in Mechanics and Mathematics[END_REF] rank sE

-A i C = n, ∀s ∈ C, ∀i = [1, . . . , h] (7)
and

rank E C = n (8)
A3. The fault f (t) satisfies f (t) ≤ α 1 and the derivative of f (t) with respect to time is norm bounded i.e. ḟ (t) ≤ α 2 and 0 ≤ α 1 , α 2 < ∞.

A4.

p + q ≤ m A5. rank E 0 C H = n + rank(H)
Note that A3 is a quite general assumption in the literature [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF], [START_REF] Jiang | Fault Accommodation for Nonlinear Dynamic Systems[END_REF]. Assumption A4 has to be verified to be able to provide actuator fault estimation.

The main goal of the paper is to estimate time-varying actuator fault for polytopic LPV descriptor systems. In [START_REF] Zhang | Adaptive Observer-based Fast Fault Estimation[END_REF], the authors have performed an adaptive observer so as to estimate time-varying faults but only for regular LTI systems. The authors in [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF] have developed a proportional integral observer for actuator fault estimation but with the following constraint ḟ (t) = 0. Moreover, the same constraint has been considered in [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor Linear Parameter Varying systems[END_REF] for discrete-time descriptor systems where the main goal was also to estimate faults. So, in this paper, the main contribution consists in designing an Adaptive Polytopic Observer for LPV descriptor systems which is able to deal with time-varying actuator faults and by the way to tackle the previous restrictive constraints of the above mentioned papers.

The following section is dedicated to the design of the Adaptive Polytopic Observer.

ACTUATOR FAULT ESTIMATION FOR POLYTOPIC LPV DESCRIPTOR SYSTEM

Adaptive Polytopic Observer Design

The proposed Adaptive Polytopic Observer (APO) has the following structure:

                   ż(t) = h i=1 ρ i (θ(t))(N i z(t) + G i u(t) + L i y(t) + F i f (t)) x(t) = z(t) + T 2 y(t) ŷ(t) = C x(t) r(t) = S(y(t) -ŷ(t)) ḟ (t) = Γ h i=1 ρ i (θ(t))U i ( ṙ(t) + σr(t)) (9)
where z(t) ∈ R n is the observer state vector, x(t) ∈ R n is the estimated state vector, r(t) ∈ R m is the residual vector, ŷ(t) ∈ R n is the estimated output vector and f (t) ∈ R p is the estimated actuator fault of f (t). N i , G i , L i , T 2 , U i , S are unknown matrices of appropriate dimensions to be determined, Γ ∈ R p×p is a learning rate symmetric positive definite matrix and σ ∈ R is a positive scalar. The state estimation error e(t) is defined as:

e(t) = x(t) -x(t)
Then, it follows from ( 6) and ( 9) that:

e(t) = (I n -T 2 C)x(t) -z(t) -T 2 Hd(t) (10)
Under assumption that rank E C = n, there exists nonsingular matrices T 1 ∈ R n×n and T 2 ∈ R n×m such that [START_REF] Darouach | Design of Observers for Descriptor Systems[END_REF]:

T 1 E + T 2 C = I n (11)
Furthermore, the fault estimation error e f (t) can be expressed as:

e f (t) = f (t) -f (t) (12)
From ( 10), the state estimation error dynamic with the relation ( 11) is given by:

ė(t) = T 1 E ẋ(t) -ż(t) -T 2 H ḋ(t) (13)
Using ( 6) and ( 9), the residual r(t) can be expressed as:

r(t) = SCe(t) + SHd(t) (14)
Then, if the following conditions hold true ∀ i = 1, ..., h:

T 1 A i -L i C -N i T 1 E = 0 (15) T 1 B i -G i = 0 (16) T 1 E = I n -T 2 C (17) T 1 R i + N i T 2 H -L i H = 0 (18) T 2 H = 0 ( 19
) SH = 0 (20)
By taking into account ( 6), ( 9) and ( 13), the estimation error dynamic ė(t) and the residual r(t) can be written such that:

ė(t) = h i=1 ρ i (θ(t))[N i e(t) + M i f (t) + F i e f (t)] (21)
r(t) = SCe(t) (22)
with

M i = (T 1 -I n )F i (23)
The substitution of ( 17) into (15) yields to:

N i = T 1 A i + (N i T 2 -L i )C ( 24
) N i = T 1 A i + K i C ( 25
)
where

K i = N i T 2 -L i (26)
To find simultaneously matrices T 1 and T 2 from equations ( 17) and [START_REF] Koenig | Unknown Input Proportional Multiple-Integral Observer Design for Linear Descriptor Systems: Application to State and Fault Estimation[END_REF], one can define the following augmented matrix: Under assumption A5 , [START_REF] Koenig | Unknown Input Proportional Multiple-Integral Observer Design for Linear Descriptor Systems: Application to State and Fault Estimation[END_REF], the solution of ( 27) can be expressed such that:

T 1 T 2 E 0 C H = I n 0 (27
T 1 T 2 = I n 0 E 0 C H + (28)
Matrix S can be determined by solving equation [START_REF] Lewis | A Survey of Linear Singular Systems Circuits Systems and Signal Processing[END_REF] by considering assumption A4 i.e. the number of disturbances and faults can not be larger than the number of measurements, a general solution for ( 20) is given by (See Chapter 4 of [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF]):

S = Π 1 [I m -H(H T H) -1 H T ] (29)
where Π 1 ∈ R m×m is an arbitrary design matrix.

Convergence analysis

Let us consider the following Lemma:

Lemma 1 Given a scalar µ > 0 and a symmetric positive definite matrix P 1 , the following inequality holds [START_REF] Zhang | Adaptive Observer-based Fast Fault Estimation[END_REF]:

2x T y ≤ 1 µ x T P 1 x + µy T P -1 1 y x, y ∈ R n (30)
In contrast to the proportional integral observer [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF], here time-varying faults are considered. Then, it follows that ḟ (t) = 0 and by consequence, the dynamic of fault estimation error is expressed as follows:

ėf (t) = ḟ (t) -ḟ (t) (31)
By using [START_REF] Rodrigues | Sensor Fault Detection and Isolation Filter for Polytopic LPV Systems: A Winding Machine Application[END_REF], the dynamic of the state estimation error (21) becomes:

ė(t) = h i=1 ρ i (θ(t))[(T 1 A i + K i C) Ni e(t) + M i f (t) + F i e f (t)] (32)
The convergence of the state estimation error (32) can be verified by the following Theorem.

Theorem 1

Under Assumptions A 1 to A 5 , given scalars σ, µ > 0, if there exists symmetric positive definite matrices Q, P 1 , P 2 , P 3 and matrices W i = QK i and U i such that, ∀ i ∈ [1, . . . , h]:

Ψ i * -1 σ F T i (QT 1 A i + W i C) -2 σ F T i QF i + 1 σµ P 2 + 1 σµ P 3 < 0 (33) s.t.
W i H + QT 1 R i = 0 (34) F T i Q -U i SC = 0 (35)
where

Ψ i = (T 1 A i ) T Q + Q(T 1 A i ) + 1 µ P 1 + C T W T i + W i C
Then, the following adaptive fault estimation algorithm

ḟ (t) = Γ h i=1 ρ i (θ(t))U i ( ṙ(t) + σr(t)) (36)
can realize e(t) and e f (t) uniformly bounded, where Γ ∈ R p×p is a symmetric positive definite learning rate matrix.

Convergence condition for the estimation error yields that the time derivative of the Lyapunov function (37) should be negative definite. The derivative of V (e(t), e f (t)) with respect to time evaluated on trajectories of the estimation error equation ( 21) is:

V (e(t), e f (t)) = h i=1 ρ i (θ(t)) e T (t)[N T i Q + QN i ]e(t) + 2e T (t)QM i f (t) +2e T (t)QF i e f (t) + 1 σ ėT f (t)Γ -1 e f (t) + 1 σ e T f (t)Γ -1 ėf (t) (38)
From ( 22) and ( 9), the equation ( 36) can be expressed as:

ḟ (t) = Γ h i=1 ρ i (θ(t))U i S ẏ(t) -ẏ(t) + σ(y(t) -ŷ(t)) (39)
and under e(t) = x(t) -x(t), (39) is equivalent to:

ḟ (t) = Γ h i=1 ρ i (θ(t))U i SC ė(t) + σe(t) (40)
By using [START_REF] Zhang | Dynamic output feedback fault tolerant controller design for Takagi-Sugeno fuzzy systems with actuator faults[END_REF] and substituting (40) into equation ( 38), one can obtain:

V (e(t), e f (t)) = h i=1 ρ i (θ(t)) e T (t)[N T i Q + QN i ]e(t) + 2e T (t)QM i f (t) + 2e T (t)QF i e f (t) + 1 σ ėT f (t)Γ -1 e f (t) + 1 σ e T f (t)Γ -1 [ ḟ (t) -ΓU i SC( ė(t) + σe(t))]
(41) Substituting ( 21) into (41) and using (35), it follows that:

V (e(t), e f (t)) = h i=1 ρ i (θ(t)) e T (t)[N T i Q + QN i ]e(t) -2 σ e T f (t)F T i QN i e(t) -2 σ e T f (t)F T i QF i e f (t) -2 σ e T f (t)F T i QM i f (t) + 2e T (t)QM i f (t) + 2 σ e T f Γ -1 ḟ (t) (42
) By taking into account assumption A3, for a positive scalar µ and symmetric matrices P 1 > 0, P 2 > 0 and P 3 > 0 and by using Lemma 1, we can obtain:

2e T (t)QM i f (t) ≤ 1 µ e T (t)P 1 e(t) + µf T (t)M T i QP -1 1 QM i f (t) ≤ 1 µ e T (t)P 1 e(t) + µα 2 1 λ max (M T i Q T P -1 1 QM i ) 2 σ e T f (t)Γ -1 ḟ (t) ≤ 1 µσ e T f (t)P 2 e f (t) + µ σ ḟ T (t)Γ -T P -1 2 Γ -1 ḟ (t) ≤ 1 µσ e T f (t)P 2 e f (t) + µ σ α 2 2 λ max (Γ -1 P -1 2 Γ -1 ) -2 σ e T f (t)F T i QM i f (t) ≤ 1 µσ e T f (t)P 3 e f (t) + µ σ f T (t)M T i QF i P -1 3 F T i QM i f (t) ≤ 1 µσ e T f (t)P 3 e f (t) + + µ σ α 2 1 λ max (M T i Q T F i P -1 3 F T i QM i ) (43)
Using (43), the time derivative of (42) can be bounded as follows:

V (e(t), e f (t))

≤ h i=1 ρ i (θ(t)) e T (t)[N T i Q + QN i ]e(t)
i µα 2 1 λ max (M T i Q T P -1 1 QM i )+ µ σ α 2 2 λ max (Γ -1 P -1 2 Γ -1 )+ µ σ α 2 1 λ max (M T i Q T F i P -1 3 F T i QM i )
Using ( 25) and for W i = QK i , the inequality (44) can be written as:

V (e(t), e f (t)) ≤ h i=1 ρ i (θ(t)){ e T (t) e T f (t) Ξ i e(t) e f (t) + δ} (45)
with

Ξ i = Ψ i * -1 σ F T i (QT 1 A i + W i C) -2 σ F T i QF i + 1 σµ P 2 + 1 σµ P 3 and Ψ i = (T 1 A i ) T Q + Q(T 1 A i ) + 1 µ P 1 + C T W T i + W i C Considering that h i=1 ρ i (θ(t)) = 1 and ρ i (θ(t)) ≥ 0 and because F i is of full column rank ∀ i ∈ [1, ..., h]
(see assumption A1), when Ξ i < 0 one can obtain that:

V (e(t), e f (t)) ≤ -ε e(t) e f (t) 2 + δ ( 46
)
where ε = min i (λ min (-Ξ i )). Then, V (e(t), e f (t)) < 0 for ε e(t) e f (t)

2

> δ, ∀t ≥ 0 which means that e(t), e f (t) converges to a small set according to Lyapunov stability theory. Now, the obtained gains matrices K i must satisfied also the constraint (18

) ∀ i ∈ [1, . . . , h].
Then, the equation ( 18) can be written as follows:

(N i T 2 -L i )H + T 1 R i = 0 (47)
and under equation [START_REF] Rodrigues | Sensor Fault Detection and Isolation Filter for Polytopic LPV Systems: A Winding Machine Application[END_REF], it can be expressed as:

K i H + T 1 R i = 0 ( 48
)
For K i = Q -1 W i ∀ i ∈ [1, . . . , h],
equation (48) becomes:

W i H + QT 1 R i = 0 (49)
Then, gains matrices K i will be obtained by solving LMIs [START_REF] Zhang | Design of Integrated Fault Detection, Diagnosis and Reconfigurable Control Systems[END_REF] under constraints [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF] and (35). ✷ Therefore, Theorem (1) implies that error (e(t), e f (t)) are uniformly bounded. Now, its easy to show that the estimated fault can be deduced from the expression (36) as follows :

f (t) = Γ h i=1 ρ i (θ(t))U i   r(t) + σ t t f r(τ )dτ    ( 50
)
where t f is the time of fault occurring. The estimated fault expression (50) combines a proportional term with an integral one. The proportional term have an influence on the rapidity of fault estimation [START_REF] Zhang | Adaptive Observer-based Fast Fault Estimation[END_REF]. The proposed APO makes possible to reconstruct the state of the system whatever the presence of the unknown inputs and actuator faults and generate a residual signal which can indicates that a fault occurs.

Remark 1: It is easy to solve inequality (33) of Theorem 1, but the solving difficulty is added due to the presence of equality constraints [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF] and ( 35). However, it is possible to transform [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF] and (35) in Theorem 1 into the following optimization problem [START_REF] Corless | State and Input Estimation for a Class of Uncertain Systems[END_REF]:

Int. J.
ηI W i H + QT 1 R i * ηI < 0 (52) ηI F T i Q -U i SC * ηI < 0 ( (51)
) 53
Remark 2:Note that, after solving the LMI problem under equalities constraints given in Theorem 1, the input-to-state convergence condition given in (46) is satisfied. Thus, in the case of timevarying faults with a bounded first time derivative, the state estimation error e(t) and the fault estimation error e f (t) converge to a ball, centered at the origin, defined by the terms δ and ε. The radius of the ball in which e(t) converges can be minimized by a choice of the parameter Γ that minimizes δ without changing ε (that does not depend on Γ). It thus improves the accuracy of the estimation.

Moreover, the APO studied in this paper can be considered as an improvement of the classical PI observer, in the sense that the convergence of the state estimation error and fault estimation error is proved (in a ball centered at the origin) even in a nonconstant fault case, whereas the assumption of a constant fault is needed to prove the convergence of the state estimation error when using a PI observer [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF], [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor Linear Parameter Varying systems[END_REF].

ILLUSTRATIVE EXAMPLE

The proposed example, considers an electrical network as shown in Figure [START_REF] Alwi | Robust Fault Reconstruction for Linear Parameter Varying Systems Using Sliding Mode Observers[END_REF] where R j , j = 1, . . . , 8 and L 1 and L 2 stand for the resistors and inductors, respectively. e 1 (t) and e 2 (t) are the voltage sources which are taken as the control inputs. We denote by i 1 (t), i 2 (t), i 3 (t) and i 4 (t) the amperage of the currents. According to the basic circuit theory and the Kirchoff's laws, we get the following Differential-Algebraic Equations (DAE) which describe the system:

L 1 L 2 R 1 R 2 R 4 R 5 R 3 e (t) 2 e (t) 1 R 8 R 7 R 6 i (t) 1 i (t) 3 i (t) 4 i (t) 2
L 1 di1(t) dt = -(R 1 + R 3 + R 5 )i 1 (t) + R 3 i 3 (t) + R 5 i 4 (t) L 2 di2(t) dt = -(R 4 + R 6 + R 7 )i 2 (t) + R 4 i 3 (t) + R 7 i 4 (t) 0 = -(R 2 + R 3 + R 4 )i 3 (t) + R 4 i 2 (t) + R 3 i 1 (t) + e 1 (t) 0 = -(R 5 + R 7 + R 8 )i 4 (t) + R 7 i 2 (t) + R 5 i 1 (t) + e 2 (t) y 1 (t) = i 1 (t) + i 4 (t) y 2 (t) = i 2 (t) + i 3 (t) y 3 (t) = i 4 (t) Let denote [ x 1 (t) x 2 (t) x 3 (t) x 4 (t) ] T = [ i 1 (t) i 2 (t) i 3 (t) i 4 (t) ] T the state vector, u(t) = [e1(t) e2(t)]
T the control inputs, y 1 (t), y 2 (t) and y 3 (t) the outputs signals, R 1 and R 6 two variable resistors. The previous set of DAE can be represented as a LPV descriptor systems as follows:

E ẋ(t) = A(θ(t))x(t) + B(θ(t))u(t) + R(θ(t))d(t) + F (θ(t))f (t) y(t) = Cx(t) + Hd(t) (54)
The resulting matrices of the LPV descriptor system are given by:

E=    1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0   , B(θ(t))=B =    0 0 0 0 1 0 0 1   , F (θ(t))= F =B, C =   1 0 0 1 0 1 1 0 0 0 0 1   , H =   2 1 1   Ã(θ(t)) =     -R11+θ1(t) L1 0 R13 L1 R14 L1 0 -R22+θ2(t) L2 R23 L2 R24 L2 R 31 R 32 -R 33 0 R 41 R 42 0 -R 44     and R(θ(t)) =    0 0 0.7 + θ 1 (t) 0.2 + θ 2 (t)    where R 11 = R 1 + R 3 + R 5 = 15Ω, R 22 = R 4 + R 6 + R 7 = 40Ω, R 33 = R 2 + R 3 + R 4 = 25Ω, R 44 = R 5 + R 7 + R 8 = 20Ω, R 13 = R 31 = R 3 = 3Ω, R 14 = R 41 = R 5 = 2Ω, R 23 = R 32 = R 4 = 5Ω, R 24 = R 42 = R 7 = 8Ω and L 1 = 0.3H, L 2 = 0.65H.
R is a disturbance distribution matrix that affects the system where d(t) is a random vector uniformly distributed in [-1, 1] which acts like additive noise. H is also a disturbance distribution matrix that affects the outputs. θ 1 (t) and θ 2 (t) are two time-varying parameters which vary according to θ 1 (t) ∈ -0.5, 0.5 and θ 2 (t) ∈ -1, 1 . The inputs control are described by

u 1 (t) = 12 sin(2.5t) u 2 (t) = 5
It is assumed that the LPV descriptor system (54) is affected by two actuator faults which followed the assumption A 3 and that are defined by:

f 1 (t) = 0.35u 1 (t) occurs in u 1 (t) when 8 ≤ t ≤ 15s f 2 (t) = 0.25u 2 (t) occurs in u 2 (t) when 20 ≤ t ≤ 30s
Note that different actuator fault can appear in voltage source: an offset for bias or temporal variation for transient fault. Such fault can appear on real source voltage with various temporal duration and different shapes.

The parameters of the polytopic LPV descriptor system evolve in a hyper-rectangle. Consequently, the matrices A i and R i of this system are given by:

A 1 =    -51
   R 1 =    0 0 0.2 -0.8    , R 2 =    0 0 0.2 1.2    , R 3 =    0 0 1.2 -0.8    and R 4 =    0 0 1.2 1.2   
The weighting functions ρ i (θ(t)) are defined as combinations of θ j and are given by [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF]:

ρ 1 (θ(t)) = θ 1 (t) -θ 1 θ 1 -θ 1 θ 2 (t) -θ 2 θ 2 -θ 2 = (θ 1 (t) + 0.5)(θ 2 (t) + 1) 2 ρ 2 (θ(t)) = θ 1 (t) -θ 1 θ 1 -θ 1 θ 2 -θ 2 (t) θ 2 -θ 2 = (θ 1 (t) + 0.5)(1 -θ 2 (t)) 2 ρ 3 (θ(t)) = θ 1 -θ 1 (t) θ 1 -θ 1 θ 2 (t) -θ 2 θ 2 -θ 2 = (0.5 -θ 1 (t))(θ 2 (t) + 1) 2 ρ 4 (θ(t)) = θ 1 -θ 1 (t) θ 1 -θ 1 θ 2 -θ 2 (t) θ 2 -θ 2 = (0.5 -θ 1 (t))(1 -θ 2 (t)) 2

APO design and residual generation

To design the presented APO according to the proposed methodology, let's verify a necessary assumption A 4 to be able to provide actuator fault estimation with equation (50). After that, a matrix S (solution of and SH = 0) can be given by [START_REF] Wei | LMI solutions to the mixed H -/H∞ fault detection observer design for linear parametervarying systems[END_REF]. Then, let's check the necessary assumption A 5 . Then, matrices T 1 and T 2 can be computed from equation [START_REF] Wang | Actuator fault diagnosis: An adaptive observer-based technique[END_REF]. After by checking A 1 and A 2 , the gains matrices can be obtained by solving the LMIs (33) under constraints [START_REF] Zhang | Bibliographical review on reconfigurable fault-tolerant control systems[END_REF] and (35) via numerical approach within the LMI Toolbox. By choosing the scalars σ = 2, µ = 0.5, these inequalities are fulfilled with: and choosing the learning rate matrix Γ = diag(1, 2), the simulation results for time-varying actuator faults estimation using the adaptive algorithm given by (50) are mentioned in the following Figures. Figures [START_REF] Bamieh | Identification of linear parameter varying models[END_REF] and [START_REF] Bokor | Detection filter design for LPV systems-a geometric approach[END_REF] show that estimated actuator faults by APO can converge to their real values which is more powerful than Proportional Integral Polytopic Observer (PIPO) presented in [START_REF] Hamdi | Fault Detection and Isolation in Linear Parameter Varying Descriptor systems via Proportional Integral Observer[END_REF]. Then, it is shown that the PIPO makes possible to decouple the disturbances while it can only detect and estimate constant actuator faults as in Figure [START_REF] Bokor | Detection filter design for LPV systems-a geometric approach[END_REF]. Indeed, its actuator fault estimation does not match correctly the real one for the time-varying case in Figure [START_REF] Bamieh | Identification of linear parameter varying models[END_REF]. In contrast, the APO shows very good results for both time-varying and constant actuator faults estimations despite the presence of disturbance: it is a good improvement for such LPV descriptor systems.

K 1 =   

CONCLUSION

In this paper, an actuator fault estimation scheme based on an Adaptive Polytopic Observer for LPV descriptor systems has been proposed. The developed strategy allows to consider not only constant faults but also time-varying actuator faults. The convergence conditions of this observer has been formulated and solved within a set of LMI under equalities constraints. The developed scheme has been applied in an electrical circuit so as to estimate both time-varying and constant actuator

)

  Int. J. Robust. Nonlinear Control (2014) Prepared using rncauth.cls DOI: 10.1002/rnc

)

  Int. J. Robust. Nonlinear Control (2014) Prepared using rncauth.cls DOI: 10.1002/rnc

Figure 1 .

 1 Figure 1. An electrical circuit

Figure 3 .

 3 Figure 3. Residual signal r 2 (t) when a fault occurs in u 2 (t)

  Int. J. Robust. Nonlinear Control (2014) Prepared using rncauth.cls DOI: 10.1002/rnc

Figure 4 .

 4 Figure 4. Actuator fault f 1 (t) and its estimated f1 (t)

Figure 5 .

 5 Figure 5. Actuator fault f 2 (t) and its estimated f2 (t)

  Consider the Lyapunov function with the following quadratic form:

	V (e(t), e f (t)) = e T (t)Qe(t) +	1 σ	e T f (t)Γ -1 e f (t)

Int. J. Robust. Nonlinear Control (2014) Prepared using rncauth.cls DOI: 10.1002/rnc Proof 1

  Figure 2. Residual r 1 (t) signal when a fault occurs in u 1 (t)

			0.2
											r 1 (t)
		0.15
			0.1
	residual signal r 1 (t)	0 0.05
		-0.05
		-0.1
		-0.15	0	5	10	15		20	25	30	35	40
										t(s)
			0.18
											r 2 (t)
			0.16
			0.14
			0.12
	residual signal r 2 (t)	0.06 0.08 0.1
			0.04
			0.02
			0
		-0.02	0	5	10	15		20	25	30	35	40
										t(s)
			21.1564 -17.1376 -25.1753			1.7542	1.8862	-5.3945	
			-16.1420 -0.9267 38.9566 -18.4702 -59.4431 33.2106	  , K 2 =	 	-3.7885 13.0178 -5.4409 5.4732 -11.0272 0.0808	 
			23.4193 -16.5291 -30.3096		1.7234	1.8499	-5.2967
			1.4095		2.0748	-4.8937			1.2491	2.5999	-5.0981
	K 3 =	 	-2.9091 12.5797 -6.7615 4.2535 -10.4436 1.9365	  , K 4 =	 	-2.7452 13.0423 -7.5518 3.9211 -10.3360 2.4938
			1.2742		2.0773	-4.6257			1.1244	2.5655	-4.8143

  

The simulation results for time varying actuator faults detection are shown in Figures

[START_REF] Armeni | Robust Fault detection and isolation for LPV systems under a sensitivity constraint[END_REF]

and (3). From figures (2) and (3), one can see that the residuals are almost zero throughout the time simulation in fault-free case. The residual signals increase in magnitude considerably when actuators faults occur for t ∈

[START_REF] Dai | Singular Control Systems[END_REF][START_REF] Hamdi | Robust Fault Detection and Estimation for Descriptor Systems Based on Multi-Models Concept[END_REF]

s and t ∈

[START_REF] Lewis | A Survey of Linear Singular Systems Circuits Systems and Signal Processing[END_REF][START_REF] Wu | Control of Parameter Varying Systems[END_REF]

s. These two faults can be estimated by using the adaptive fault estimation algorithm as in (50). Int. J. Robust. Nonlinear Control (2014) Prepared using rncauth.cls DOI: 10.1002/rnc

  Prepared using rncauth.cls DOI: 10.1002/rnc faults for polytopic LPV descriptor systems. A comparison with previous works like PI observers, underlines the theoretical and practical improvements for time-varying actuator faut estimation for such LPV descriptor systems.

Int. J. Robust. Nonlinear

Control (2014)

Int. J. Robust. Nonlinear Control (2014) Prepared using rncauth.cls DOI: 10.1002/rnc


Alexandre Bilger

Éric Bardinet

Sara Fernandez-Vidal

Christian Duriez

Pierre Jannin

Stéphane Cotin

Intra-operative Registration for Stereotactic Procedures driven by a combined Biomechanical Brain and CSF Model

During stereotactic neurosurgery, the brain shift could affect the accuracy of the procedure. However, this deformation of the brain is not often considered in the pre-operative planning step or intraoperatively, and may lead to surgical complications, side effects or ineffectiveness. In this paper, we present a method to update the preoperative planning based on a physical simulation of the brain shift. Because the simulation requires unknown input parameters, the method relies on a parameter estimation process to compute the intracranial state that matches the partial data taken from intra-operative modalities. The simulation is based on a biomechanical model of the brain and the cerebro-spinal fluid. In this paper, we show on an anatomical atlas that the method is numerically sound.

Introduction

Stereotactic neurosurgery requires a high precision for target location and target definition. For example, one of the main targets in deep brain stimulation, the subthalamic nucleus, is about 9 millimeters long. However, a combination of deformation and motion of the brain, known as brain shift, can occur during the procedure, depending on the surgical technique and the patient's brain anatomy and pathology. In stereotactic procedures, it may cause a displacement of the target or other structures. Elias et al. [START_REF] Elias | Cortical and subcortical brain shift during stereotactic procedures[END_REF] reported an anterior commissure shift up to 5.67 mm with a mean of 0.98 mm, and higher shift values for the frontal pole.

Because brain shift is a major source of errors when applying the planned strategy intra-operatively, several groups have studied the problem of brain shift compensation using a deformable model, instead of using a fully intensity-based method. Skrinjar et al. [START_REF] Skrinjar | Model-driven brain shift compensation[END_REF] presented a method to deform pre-operative data according to a partial intra-operative brain surface, captured by a stereo camera, while Chen et al. [START_REF] Chen | Intraoperative brain shift compensation: accounting for dural septa[END_REF] and Audette et al. [START_REF] Audette | An integrated range-sensing, segmentation and registration framework for the characterization of intra-surgical brain deformations in image-guided surgery[END_REF] used a laser range scanner. However the equipment is not standard and requires a large working space without occlusion, which is not always the case in neurosurgery. The work of Wittek et al. [START_REF] Wittek | Patient-specific model of brain deformation: application to medical image registration[END_REF] Fig. 1: Pre-operative (left) MR image without brain shift, and post-operative (right) MRI showing asymmetric brain shift. The green and pink contours is a segmentation of the pneumocephalus (air inside skull). The pink line represents the surface of the brain in contact with air. is very similar to [START_REF] Skrinjar | Model-driven brain shift compensation[END_REF] but they used intra-operative MRI to guide the deformable model and took extra care on the complex model of the brain tissue deformation. However [START_REF] Skrinjar | Model-driven brain shift compensation[END_REF] and [START_REF] Wittek | Patient-specific model of brain deformation: application to medical image registration[END_REF] do not model physically the brain shift phenomenon, because they include artificial forces by adding virtual springs or constraints between pre-and intra-operative control points. None of these methods accounts for the effects of gravity, interaction of CSF and brain tissue and loss of CSF. On the contrary, Chen et al. [START_REF] Chen | Intraoperative brain shift compensation: accounting for dural septa[END_REF] used a computational model which accounts for CSF loss and gravity. They pre-operatively built a statistical atlas of deformation to solve the inverse problem intra-operatively.

We propose a new method to take into account brain shift during a stereotactic neurosurgical procedure. It relies on high resolution pre-operative MR image and intra-operative image (CT scan) of the patient's brain, acquired once the brain has shifted, but the method could be extended to other modalities such as intra-operative MRI. The pre-operative data are deformed according to a physical model of brain shift, mainly based on a cerebro-spinal fluid (CSF) interaction with the brain tissue, gravity and brain-skull interactions. To best fit the intra-operative configuration, the input parameters of the brain shift model are estimated intra-operatively with a measure of similarity with a series of precomputed deformed configurations. The output of this inverse problem gives parameters used to compute the intra-operative configuration, as detailed as the pre-operative data. With this information, the surgeon is able to update the pre-operative planning data (target coordinates, trajectory angles etc).

The following section describes the brain shift simulation and how we use it to solve the intra-operative registration based on a physical parameters estimation. Section 3 presents the results and the numerical validation on an anatomical atlas. Section 4 concludes and addresses future steps for the method.

The overall pipeline of our method is detailed in the figure 2. This section presents the model of brain shift used in the simulations and the parameters estimation (the highlighted steps in the figure 2). The other aspects are not treated in this paper as we focus on the numerical validation. Fig. 2: Steps of our method (read from left to right, then from top to bottom): the acquired pre-operative MRI enables to generate a simulation. A set of deformation is pre-computed. During the surgery, the pneumocephalus is extracted from the intra-operative CT, rigidly registered to the pre-operative data to be in the same space and compared to the set of pre-computed deformation. Finally, the displacement of the structures of interest is deduced. Among these steps, the most time consuming processes are done before the surgery.

Presentation of the model

The brain shift simulation relies on a physics-based model of the brain tissue deformation, the contact response with the skull and the falx cerebri, and the interaction with the CSF. It is based on the work of [START_REF] Bilger | Computation and visualization of risk assessment in deep brain stimulation planning[END_REF].

The brain deformation is computed using a non-linear geometric finite element method, with a linear constitutive law [START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. Theory[END_REF]. This allows for rotation in the model, while relying on a linear expression of the stress-strain relationship. The viscous part of the brain behavior is omitted, because we focus on the static equilibrium after the brain has shifted. The equation relating the external forces f to the nodal displacements u can be written as f = K(u)u, with K the stiffness matrix depending on u. Both hemispheres are meshed independently with linear tetrahedrons.

When the brain deforms, the simulation algorithm solves the contacts with the skull and the falx cerebri, considered as rigid. Unilateral interpenetration constraints are created when a collision is detected. The area near the brainsterm is assigned fixed Dirichlet boundary conditions to impose a null displacement constraint. The contacts are solved with Lagrangian multipliers. Figure 3 shows the boundary conditions of the model.

The main cause of brain shift is a CSF loss [START_REF] Slotty | The impact of brain shift in deep brain stimulation surgery: observation and obviation[END_REF]. The resulting force of the surrounding fluid acting on the brain tissue acts against the weight, so that the brain is balanced. A CSF loss causes a fluid forces decrease, breaking the balance and leading to a brain shift. In order to cause a brain shift in our simulation, we include a model of CSF forces:

f CSF = S ρ g h(P ) dS
where ρ is the density of CSF (1000 kg/m 3 ), g is the norm of the gravity and h is the distance between a point P on the mesh and the fluid surface. The force f CSF is applied onto every immersed triangle S of the brain surface mesh. Asymmetric brain shift is observed on post-operative images (see Fig. 1). To handle this property, two independent models of fluid forces are present in the simulation, each acting independently on each hemispheres. With this model, the deformation depends also on the patient's head orientation compared to gravity direction, mechanical properties of the tissue and patient variability.

To recover the true rest configuration of the brain before applying gravity and external forces on the finite element model, we solve an inverse problem using an iterative geometric algorithm [START_REF] Sellier | An iterative method for the inverse elasto-static problem[END_REF]. The algorithm is used with a static solver. The iterations are smoothed to handle the geometrical non-linearities of the deformation model.

Non-rigid registration based on physical parameters estimation

Parameters to estimate As explained in section 2.1, the brain shift model depends on the following parameters: patient's head orientation compared to gravity direction, patient variability (geometry), mechanical properties of the tissue and the amount of CSF lost during the procedure. The two first parameters are known: we can measure the orientation of the patient's head and the geometry is acquired by the segmentation of patient images. The mechanical parameters of the brain tissue have been estimated by several groups using different techniques, but there is no consensus and there is no direct measurement technique personalized for a patient. Although it would be interesting to add Fig. 3: Schematic representation of the brain model and the set of constraints describing its interactions with the environment. In this pre-operative configuration, the patient is in the supine position (see the gravity vector direction -→ g ).

The cerebro-spinal fluid surrounds the brain tissue and acts on it with pressure forces. The resultant force balances the weight (see the ---→ f CSF compared to the weight force m -→ g ). The brain is under a null displacement constraint near the brainstem area. The illustration also shows the contacts between the brain and the skull.

the mechanical properties to the parameters to estimate, the deformation is also controlled by the CSF volume lost, meaning that these two parameters cannot be dissociated. That is why the mechanical properties are fixed to the value used in [START_REF] Ferrant | Serial registration of intraoperative MR images of the brain[END_REF]. Moreover, the mechanical parameters are not the physical cause of the brain shift phenomenon. We will show that an acceptable error in the mechanical parameters does affect the parameter estimation, but not the output of our algorithm. Finally, only the CSF volume lost is estimated. As we model both hemispheres independently, two quantities have to be estimated.

Error minimization In our algorithm, the only unknown input parameters are two volumes of CSF lost during the surgery. To estimate these parameters, we introduce a measure of similarity, to compare the final deformed geometry of the brain, obtained at equilibrium after applying a CSF loss in a simulation, and the visible surface of the brain, extracted from an intra-operative modality. The similarity between these two objects is defined with a least-squares approach:

d(S, M ) = x S ∈S x S -p(x S , M ) 2 (1)
with x S a point on the surface S of the intra-operative brain, p(x S , M ) is the projection of the point x S on the surface M of the simulated brain. The measure is normalized by dividing by the number of points in S, and represents the average distance between both surfaces. The figure 4 shows a representation of such a projection in image space.

Our model is consistent with the correlation between the CSF loss and the amount of deformation, studied by [START_REF] Elias | Cortical and subcortical brain shift during stereotactic procedures[END_REF]. The displacement norm of a point in the brain is an increasing function of the CSF volume lost. This property guarantees the convexity of the error as a function of the two CSF volume losses (see an example in Fig. 5): a local minimum is a global minimum.

Even an efficient optimization method for our problem would require at least several tens of iterations to converge to an estimation of the parameters. Each iteration requires to compute a direct simulation, which takes time (more than 5 minutes) and therefore is not compatible with an application during a neurosurgical procedure. To speed up the intra-operative computation, we pre-compute before the surgery a set of brain configurations, deformed by a simulated brain shift for a regular sampling of the parameters domain. Then, during the procedure, a mesh is extracted from an intra-operative image, and compared to all the pre-computed deformed configurations. The configuration with the smallest error corresponds to the registered brain.

Finally, the displacement field enables to compute rapidly the displacement and deformation of other structures such as the ventricles, the target or the blood vessels. To validate our method, we first apply it on a template containing high-fidelity models of the brain, skull, skin, ventricles and blood vessels. In our tests, both hemispheres are meshed with approximately 10,000 tetrahedrons. The Young's modulus E is set to 3000 Pa and Poisson's ratio to 0.45 according to [START_REF] Ferrant | Serial registration of intraoperative MR images of the brain[END_REF].

In this section, we show that our method is numerically sound. A virtual brain shift, obtained with fixed chosen parameters, is applied on the template. The obtained configuration will be used as a synthetic of 3D intra-operative data. These data are then used in the parameters estimation process to prove the consistency of the method: the estimated parameters are compared to the parameters used to generate the data.

Finally, we also show that if the Young's modulus differs slightly from the mechanical parameters used to generate the synthetic intra-operative mesh, our algorithm still register the intra-operative brain with a minimal error.

Data Generation With given input parameters V L and V R (CSF loss volume for both sides), a deformation is computed with the brain shift simulation: the pre-operative configuration undergoes a CSF leak leading to a brain shift. After the brain shift, when the deformed brain is at equilibrium, the surface in contact with air is extracted, i.e. the same data we extract from intra-operative patient data (see pink line in Fig. 1). The surface is uniformly remeshed with a marching cube algorithm to mimic data coming from 3D images, and noise is added, to meet as much as possible the intra-operative conditions. Finally, the parameters leading to this surface are estimated.

Validation protocol Multiple intra-operative data were generated in order to validate different scenarios: symmetric brain shift(#1), asymmetric brain shift (#2), unknown mechanical parameters for both symmetric (#3) and asymmetric brain shifts (#4). The parameters (CSF loss volumes and Young's modulus) of each scenario are summarized in the table 1.

Input parameters

Optimization output Volumes(cm 3 ) Young's modulus(Pa) Volumes(cm Fig. 5: The CSF loss volume are varying (X and Y axis). The colors represents the average distance between the equilibrium state of the deformed brain caused by a CSF loss (x,y). The minimum is located at (1 × 10 -5 , 2 × 10 -5 ), which is the set of parameters used to generate the intra-operative data (scenario #2).

Error measurement For each scenario, a set of deformation is computed, and compared to the synthetic intra-operative model, in order to find the minimum distance. The figure 5 shows the average distance between data generated with the scenario #2 and the simulated brain shifts, with CSF losses varying between 0 and 30cm 3 .

The table 1 summarizes the error between the values used to generate the synthetic intra-operative model, and the values estimated by the minimization process. These values are comparable only if the mechanical parameters used for generation and optimization process. To measure errors if the mechanical parameters are different, we compare directly the brain geometry with the measure of similarity presented in the section 2.2.

The data in the table 1 show that, in the scenarios #1 and #2, the volume parameters are the same in input and output (the error in the measure of similarity comes from the remeshing and noise in the intra-operative data). In scenarios #3 and #4, where the brain is stiffer in the parameters estimation, the volumes estimated are greater, as expected, but the brain deformed with these parameters still matches the intra-operative data. We conclude that the Young's modulus does not have an effect on the registration output.

The figure 6 shows an example of the difference between the pre-operative configuration of the brain and the deformed brain fitting the intra-operative data. With our method, we are able to estimate the CSF volume lost, which gives a registration of the brain if used in the brain shift model. The entire process is based on a physical model, and no artificial forces are needed. The registration algorithm also provides the deformation and displacement of brain structures thanks to the displacement field.

Conclusion

In this paper, we presented a method to register the pre-operative configuration of the brain onto the intra-operative state. The registration is entirely based on physics and handle patient-specific geometry, patient's head orientation and asymmetry in the deformation. We showed the numerical consistency of the method on 4 scenarios involving symmetry/asymmetry and mechanical parameters changes. The algorithm provides also the physical displacement and deformation of other structures such as the blood vessels or the target. The method is independent from the constitutive law of the brain, and could be improved using a more accurate law. The next step of this work is to validate it on a series of patient.

  (a) Measure of similarity: the intra-operative triangular surface (pink) is projected on the simulated brain surface (yellow). The measure is the average distance of projection. (b) At each time step of a simulation, the average distance between the simulated brain and the intra-operative data is measured. In this example, the simulation is parameterized with the same values used to generate the intra-operative data, that is why the distance converges to a small value.

Fig. 4 :

 4 Fig. 4: Mesure of similarity between the simulated brain state and intra-operative data.

Fig. 6 :

 6 Fig. 6: 3D view of the mechanical model of the brain (in blue, only cut for visualization purpose). In pink, the mesh corresponding to intra-operative data. Left: the mechanical model is in its pre-operative configuration. Right: The brain fits the intra-operative data after a deformation

Table 1 :

 1 Different scenarios of brain shift used to validate the method. The scenarios #2 and #4 mimic an asymmetric brain shift. #3 and #4 use data generated with mechanical parameters different from the minimization process.

	3 ) Average distance (mm)

Acknowledgment

The authors would like to thank the French Research Agency (ANR) for funding this study through the ACouStiC project.