Carlos Rengifo

Franck Plestan
email: franck.plestan@irccyn.ec-nantes.fr

Yannick Aoustin

Optimal control of a neuromusculoskeletal model: a second order sliding mode solution

The objective in neuromusculoskeletal simulation consists in finding a set of muscles excitations in order to produce the desired motion and is an exciting challenge for medical operators. In this paper, muscle excitations are computed by using a second order sliding mode controller which is optimal with respect of a functional cost based on the tracking error and the torque square. The proposed strategy was tested using a two-dimensional human anthropomorphic arm composed of two joints and six muscles.

I. INTRODUCTION

Computer simulation has become a very useful tool to investigate how the neuromusculoskeletal system interacts to produce coordinated motions and provides possibilities which are not generally possible by experimental ways. As a matter of fact, medical staff could use this simulator in order to evaluate the effect of their actions on a real patient or to analyze the behaviour of a patient who is deficient in his motion. Thus, individual muscular signals like neural excitations or forces can be perturbed in their incidence on the motion [START_REF] Komura | Simulating pathological gait using the enhanced linear inverted pendulum model[END_REF]. In this case, neuromusculoskeletal simulation is used to understand muscle coordination of human legs [START_REF]Understanding muscle coordination of the human leg with dynamical simulations[END_REF], to analyze normal and pathological gait [START_REF] Piazza | Muscle-driven forward dynamic simulations for the study of normal and pathological gait[END_REF] and to estimate muscle forces [START_REF] Buchanan | Neuromusculoskeletal modeling: Estimation of muscle forces and joint moments and movements from measurements of neural command[END_REF]. The role of closed loop control techniques in neuromusculoskeletal model simulation is crucial, as the following survey will display. One of the formerly works in neuromusculoskeletal simulation [START_REF] Chow | Studies of human locomotion via optimal programming[END_REF] (forward dynamics simulation) consists in finding a set of muscles excitations which, when applied to forward dynamics equations, produces desired motions. In [START_REF] Anderson | Dynamic optimization of human walking[END_REF] the muscle excitations are taken as the optimization variables for minimizing the metabolic energy expended per traveled distance unit. This approach unfortunately is computationaly expensive: for example, muscle excitations computations for a 3D human walking model require more than 10000 hours of computation with a supercomputer [START_REF] Anderson | Dynamic optimization of human walking[END_REF]. Another proposition, called inverse dynamics [START_REF] Crowninshield | A physiologically based criterion of muscle force prediction in locomotion[END_REF], estimates muscle forces from joint torques, which are computed as a function of joint positions, velocities and accelerations. As muscles number is greater than joints number, different set of muscular forces could give the same joint torques. Several criteria have been proposed in order to solve the redundancy problem: among them, the minimum sum of muscle efforts [START_REF] Seireg | Biomechanical Analysis of the Musculoskeletal Structure for Medicine and Sports[END_REF], the minimum sum of muscle stresses [START_REF] Crowninshield | A physiologically based criterion of muscle force prediction in locomotion[END_REF] and the minimization of muscle work. Although inverse dynamics approach is a computationally efficient technique, its main drawback is its lack of robustness [START_REF] Silva | Sensitivity of the results produced by the inverse dynamic analysis of a human stride to perturbed input data[END_REF]. The second generation of neuromusculoskeletal simulation approaches is using closed-loop control strategies in order to compute muscular excitations. The controller input is the difference between desired and simulated motions, the output being muscle excitations vector. As stated in [START_REF] Seth | A neuromusculoskeletal tracking method for estimating individual muscle forces in human movement[END_REF], the benefit of this tracking formulation is to constraint the model to produce more realistic results; however, the question on how neural central system is solving the redundancy problem still remains as an open problem. Three main strategies have been proposed, each one solving this redundancy problem by a different manner. In [START_REF] Thelen | Generating dynamic simulations of movement using computed muscle control[END_REF], computed muscle control is an extension of the computed torque control idea [START_REF] Lewis | Control of Robot Manipulators[END_REF]: the minimization criteria is the square of euclidean norm of steady state muscular activations vector. Neuromusculoskeletal tracking [START_REF] Seth | A neuromusculoskeletal tracking method for estimating individual muscle forces in human movement[END_REF] solves the redundancy problem using feedback linearization and linear quadratic control. Recently, sliding mode control has been also applied to neuromusculoskeletal simulation. In [START_REF] Lister | Simulation of human gait and associated muscle activation strategies using sliding-mode control techniques[END_REF], a strategy based on inverse dynamic and first order sliding mode control is used to track desired motion. The problem with first order sliding mode is that, before the system reachs the sliding surface, the insensitivity to variations of system parameters and external disturbances can not be ensured [START_REF] Utkin | Integral sliding mode in systems operating under uncertainty conditions[END_REF] and the convergence to origin is only asymptotic. The three aforementioned strategies solve in an effective way the tracking motion problem but use very artificial optimization criteria from physiological pointof-view. When desired motions are derived from experimental data, another shortcoming arises because in these works the desired movement is specified using positions, velocities and accelerations: however, it is well-known that second order time derivatives estimation of data signals is highly sensitive to noise. The purpose of this paper is to present an application of second order sliding-mode techniques to neuromusculoskeletal simulation. As the used controller [START_REF] Plestan | A new algorithm for high order sliding mode control[END_REF] is based on a high-order sliding mode approach [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF], the robustness is ensured during the entire response as the sliding manifold is designed such that the system is evolving on this manifold early from initial time. Furthermore, the convergence time is finite and well-known in advance, and the reference trajectories are specified using only joint positions and velocities (and not accelerations, which is a positive point with respect to previous approaches), and the minimization criteria for solving the redundancy problem has a strong physiological support. The paper is organized as follows: The human arm model equations are given in Section II, as desired motions design. Section III describes the control strategy. Section IV presents the simulation results. The final section V contains conclusion and perspectives.

II. ANTHROPOMORPHIC ARM MODEL

The driven system is depicted in Figure 1. The arm model is composed of 6 muscles and 2 joints and it is assumed the arm is evolving in a vertical plane.

A. Skeletal dynamics

Using the second Lagrange method, skeletal dynamics equations are obtained

D (q 2 ) q + C (q, q) q + G(q) + τ f ( q) = Γ (1)
with q = q 1 q 2 T the generalized coordinates vector 1 , D(q 2 )(2 × 2) the symmetric positive inertia matrix, C(q, q)(2 × 2) the Coriolis and centrifugal effects matrix and G(q)(2 × 1) the gravity effects vector. Γ = Γ 1 Γ 2 T is the torques vector applied at the links joints and τ f ( q) = τ f1 τ f2 T is the friction effects vector (j = {1, 2}) such that τ fj = B j qj + F sj sign ( qj ). B j and F sj are Coulomb and viscous friction coefficients for the considered links joint. Arm parameters are given in [START_REF] Stroeve | Impedance characteristics of a neuromusculoskeletal model of the human arm i: Posture control[END_REF].

B. Torque-force relation

As no muscle is considered on the second link (Figure 1), relation between torques and forces is given by a constant 1 Notation T denote matrix transposition.

matrix

Γ 1 Γ 2 = r 11 . . . r 16 r 21 . . . r 26 R •    F 1 . . . F 6    (2)
with F i the force developped by muscle i. Coefficients of R matrix are given in [START_REF] Stroeve | Impedance characteristics of a neuromusculoskeletal model of the human arm i: Posture control[END_REF].

C. Muscular forces

In the sequel, each muscle is supposed to be an instantaneously force generator, which means that the muscular force is a nonlinear algebraic function of the muscular activation a(t). Furthermore, each tendon is supposed to be rigid, which means that its length is constant. Then, each muscle of the arm is characterized through the active muscular force which it is able to engender, this force being given from (i {1, • • • , 6})

F ai = a i f li (l i ) f vi li , v max (a i , l i ) F maxi (3)
with l i the length of muscle i, F max the maximal value of muscular effort. The functions f l (.) and f v (.) and the parameters of the muscles are described in [START_REF] Stroeve | Impedance characteristics of a neuromusculoskeletal model of the human arm i: Posture control[END_REF].

D. Musculo-skeletal geometry

Muscles are now linked to external supports. Of course, it means that their lengths depend on the articular positions.

The relation between muscular fiber length and joint positions reads as

L = L r -R T (q -Q r ) (4)
L = l 1 . . . l 6 T is the vector of muscular fiber lengths, L r = l r1 . . . l r6 T is the vector of muscle length in rest position, and and Q r = Q r1 Q r2 T is the vector of rest position for every articular joint i. As the muscular force depends also on fiber length contraction velocity, the relation between these variables and joint velocities reads as L = -R T q.

E. Activation dynamics

The activations dynamics displays the relation between the muscular excitation u i and the mechanical activation a i and reads as

ȧi =    (u i -a i ) [u i /τ ai + (1 -u i )/τ di ] , u i ≥ a i , (u i -a i )/τ di , u i < a i , (5)
with τ ai and τ di respectively the time constants for the activation and desactivation dynamics. Excitation and activation levels are allowed to vary continuously between 0 and 1 [START_REF] Zajac | Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control[END_REF].

F. Simulation model

Simulation model is displayed in Figure 2 and is based on previous equations. Let x (resp. U ) denote the state (resp. control input) vector:

x s = x s1 . . . x s10 T = q 1 q1 q 2 q2 a 1 . . . a 6 T , U = u 1 . . . u 6 T .
Then, the arm model reads as a nonlinear system which is not affine in control input,

ẋs = f (x s , U ) (6)
As as matter of fact, the system can not be written as a nonlinear model affine in U (which would be more convenient for control design) because activation dynamics (5) depends on u i and u2 i . As the objectice consists in forcing the arm to track articular positions trajectories, the system output is defined as y = x 1 x 3 T .

Dynamics

U = [u 1 • • • u 6 ] T
is the muscular excitations vector (and then the control input vector

). A = [a 1 • • • a 6 ] T is the muscular activations vector. F = [Fa 1 • • • Fa 6 ] T is the active muscular forces cvector. L = [l 1 • • • l 6 ] T is the muscular fiber length vector. Γ = [Γ 1 Γ 2 ]
T is the joint torques vector. q = [q 1 q 2 ] T is the joint positions vector.

G. Control synthesis model

Muscle dynamics are stable 2 and faster than arm motions, which yields that muscle dynamics can be neglected. Then, the model, which will be used for control design, is built by supposing that there is no dynamics between A and U (see Figure 2), i.e. A = U . Then, the system can be written as

ẋ = f (x) + g(x)U (7)
with x = x 1 . . . x 4 T = q 1 q1 q 2 q2 T , U = u 1 . . . u 6 T .

III. CONTROL STRATEGY

The proposed control strategy is a hierarchical one composed by four levels (Figure 3), the output of each level being a reference signal for the next one. It is recalled that the control objective consists in tracking articular reference trajectoire q ri (t) by acting through the excitation of the six muscles (through input U ). An other objective consists also in designing robust controller with respect to parameters variations (for example, mass which would be adding on the forarm as an object taking in the hand).

Excitations Activations Forces Desired Desired Torques Desired A r F r U Γ r q, q q r , qr
Eq.( 10) Eq.( 8) Eq.( 12)

Eq.( 13) Fig. 3. Proposed control scheme.

A. Control scheme

The four levels of the controller scheme can be detailed as follows:

First level. This block computes the desired torques Γ r using computed-torque formulation [START_REF] Lewis | Control of Robot Manipulators[END_REF] and robust control (9)

Γ r = D (q 2 ) η + Ĉ (q, q) q + Ĝ(q) + τf ( q) (8)

η = η 1 η 2 T := f c (q r , qr , q, q) (9)
with D, Ĉ, Ĝ, τf the estimated values of D, C, G and τ f (1) (knowing that these terms have uncertainties). In Section III-B, a solution for η based on second order sliding mode controller is proposed and compared (in simulation) with a classical proportional-derivative control. Note that these both controllers are computed in order to minimize a criteria such that the torques are minimal. Second level. Desired muscular forces F r = F a1 r . . . F a6 r T are computed as a function of desired torques using the linear relation (from equation ( 2))

Γ r = RF r . (10)
This system having six unknowns and two equations, there exists an infinity of exact solutions. As a consequence, vector F r is chosen in order to minimize muscle efforts. Then, at each simulation step, the following constrained optimization problem is solved using linear programming

min Fa i r >0 6 i=1 F ai r F maxi 2 (11)
with F maxi being the maximum value allowed for muscle i (for numerical values, see [START_REF] Stroeve | Impedance characteristics of a neuromusculoskeletal model of the human arm i: Posture control[END_REF]). Third level. Desired activations A r = a 1r . . . a 6r T are computed in order to obtain desired forces vector F r . For each muscle, a ri (i = {1 . . . 6}) is obtained as a function of F ai r using contraction dynamics (3). For a sake of clarity, equation ( 3) is rewritten for each muscle i by using the reference activations (with i = {1 . . . 6})

F ai r = a ri F l (l i ) F vi li , v max (a ri , l i ) F maxi (12)
In this case, one gets a system with six unknowns and six independent nonlinear equations. Find an explicit solution for a ri is a hard task: solutions of system ( 12) are numerically calculated by using a root solver. Fourth level. Muscle excitations input vector U = u 1 . . . u 6 T is computed as a function of the desired activations. Given the activation dynamics equation ( 5), for a constant excitation u(t) = ū, the corresponding steady-state activation ā equals ū. Thus, excitation inputs are computed by using the trivial equation (with i = {1 . . . 6})

u i = a ri . (13)

B. Robust controller

In this section, two different controllers are presented. The first one is based on a second order sliding mode control approach, the second one using a proportional-derivative control and being designed in order to offer performances comparisons with respect to the previous one. The both controllers are synthetized from model displayed in Section II-G whereas they are evaluated on a simulator based on the model displayed in Section II-F. Controller parameters are determined by minimization of the following performance index

J = t f 0 (q -q r ) T (q -q r ) e t + Γ T QΓ dt (14)
with

Γ T QΓ = Γ 1 Γ ss1 2 + Γ 2 Γ ss2 2
The exponential term e t gives a strong weight to tracking errors in latter time interval. The weight matrix Q was chosen to penalize large control signals. Q allows to compare each joint torque with respect to its steady-state value. Γ ss = Γ ss1 Γ ss2 T = G (q f ) is obtained by evaluating skeletal dynamics at t = t f (q (t f ) = q f , q (t f ) = 0, q (t f ) = 0).

C. Sliding mode control

The feedback control ( 9) is defined by the second-order integral sliding control law proposed by [START_REF] Plestan | A new algorithm for high order sliding mode control[END_REF]. Let s i (i = {1, 2}) denote the so-called sliding variable defined as s i = q i -q ri (t).

The objective of a second order sliding mode controller consists in forcing the system to track trajectories such that, in finite time and in spite of uncertainties, one has s i = ṡi = 0 for t > t c .

1) Some recalls: Define s = [s 1 s 2 ] T . From ( 1) and ( 8) one gets, s = q -qr (t) = D -1 (q 2 ) D (q 2 ) η + Ĉ (q, q) q + Ĝ(q)+ τf ( q) -C (q, q) q -G(q) -τ f ( q)] -qri (t)

:= χ(•) + β(•)η (15)
H1. The solutions of [START_REF] Stroeve | Impedance characteristics of a neuromusculoskeletal model of the human arm i: Posture control[END_REF] are understood in the Filippov sense [START_REF] Filippov | Differential equations with discontinuous right-hand side[END_REF], and system trajectories are supposed to be infinitely extendible in time for any bounded Lebesgue measurable input. H2. β is a dominant diagonal matrix for x ∈ X ⊂ IR 4 , X being a bounded open subset of IR 4 within which the boundedness of the system dynamics is ensured.

Last assumption implies that system (15) is "sufficiently" decoupled which allows to consider that input η i is acting on sliding variable s i and then to design by a separate way each component of input vector. Then, for a sake of clarity, one supposes only the case of monovariable system. Consider now single input equation

si = χ i (•) + β i (•)η i (16)
H3. Functions χ i (•) and β i (•) are bounded uncertain functions, and, without loss of generality, let also the sign of the control gain β i be constant and strictly positive. Thus, there exist

K m ∈ IR + * , K M ∈ IR + * , C 0 ∈ IR + such that 0 < K m < β i < K M |χ i | ≤ C 0 . (17)
for x ∈ X.

The synthesis of a high order sliding mode controller for ( 16) is made through the following idea: switching variable is defined such that the system evolves, early from t = 0, on a switching manifold. Furthermore the sliding variable and its time derivatives reach the origin in finite time in spite of uncertainties thanks to discontinuous control input. The design of the controller consists in two steps

• Design of the switching variable for system (16),

• Design of a discontinuous control input η i maintaining the system trajectories on the switching manifold which ensures the establishment of a second order sliding mode, in spite of uncertainties.

Switching variable.

Let S i denote the switching variable defined as

S i = ṡi -Ḟi (t) + λ 0i [s i -F i (t)] , (18)
with λ 0i > 0. The function F i (t), called pre-computed trajectory, is at least a C 2 -one defined such that S i (t = 0) = 0 and s

(k) i (t = t c ) -F (k) i (t c ) = 0 (k = {0, 1}).
Then, from initial and final conditions, the problem consists in finding the function F i (t) such that

s i (t = 0) = F i (0), s i (t = t c ) = F i (t c ) = 0, ṡi (t = 0) = Ḟi (0), ṡi (t = t c ) = Ḟi (t c ) = 0 (19)
A solution for F i (t) reads as

F i (t) = a 3,i t 3 + a 2,i t 2 + a 1,i t + a 0,i (20)
H4. There exists a finite positive constant

Θ i ∈ IR + such that, for 0 ≤ t ≤ t c , si -Fi (t) + λ 0i ṡi -Ḟi (t) < Θ (21)
Equation S i = 0 describes the desired dynamics which satisfy the finite time stabilization of vector [ ṡi s i ] T to zero. Then, the switching manifold on which system ( 16) is forced to slide on via a discontinuous control η i , is defined as

S i = {x | S i = 0} (22)
Given equation [START_REF]Understanding muscle coordination of the human leg with dynamical simulations[END_REF], one gets S i (t = 0) = 0: at the initial time, the system still evolves on the switching manifold.

Controller design.

The attention is now focused on the design of the discontinuous control law u which forces the system trajectories of ( 16) to slide on S i , to reach in finite time the origin and to maintain the system at the origin.

Theorem 1: [START_REF] Plestan | A new algorithm for high order sliding mode control[END_REF] Consider the nonlinear system ( 16) with a relative degree 2 with respect to s i (x, t). Suppose that it is minimum phase and that hypotheses H 1 , H 2 , and H 3 are fulfilled. Let 2 be the sliding mode order and 0 < t c < ∞ the desired convergence time. Define S i ∈ IR by

S i = ṡi -Ḟi (t) + λ 0i [s i -F i (t)] 0 ≤ t < t c ṡi + λ 0i s i t ≥ t c (23
) with pre-computed trajectory derived from [START_REF]Understanding muscle coordination of the human leg with dynamical simulations[END_REF] such that assumption H4 is fulfilled. The control input η i defined by

η i = -k i sign(S i ) (24)
with

k i ≥ C 0 + Θ + ǫ K m , (25)
C 0 , K m defined by [START_REF] Utkin | Integral sliding mode in systems operating under uncertainty conditions[END_REF], Θ defined by (21), ǫ > 0, leads to the establishment of a second order sliding mode with respect to s i . The convergence time is t c .

Sketch of proof.

Condition (25) allows to satisfy the ηattractivity condition ṠS ≤ -ǫ|S|. For more details, see [START_REF] Plestan | A new algorithm for high order sliding mode control[END_REF].

2) Application:

F i (t) = a 3,i t 3 + a 2,i t 2 + a 1,i t + a 0,i is computed such that F i (0) = e i (0), F i (t c ) = 0, Ḟi (0) = ėi (0), Ḟi (t c

D. Proportional-derivative control

In order to evaluate performances (and its interest versus a "classical" controller) of previous second order sliding mode controller, the feeedback control law (9) reads as a proportional-derivative control equation (with s i (t) = q i (t) -

IV. SIMULATIONS

In this section simulation results are displayed. In order to evaluate control performance and robustness versus parameters variations, two sets of simulation have been made: the first one is made by supposing that all parameters are fully well-known (called "nominal case"), whereas the second one supposes an error in the upper-arm mass of 50% (called "uncertain case"). For the second order sliding mode control, positions, torques and excitations for both nominal and uncertain cases are respectively displayed in Figures 4, and5, whereas, for PD control, results are compared in Figures 6 and7. These simulations clearly show that second order sliding mode controller is more robust to mass variations than PD one: furthermore, the finite time feature of second order sliding mode control clearly appears.

V. CONCLUSION AND PERSPECTIVES

A new neuromusculoskeletal simulation technique based on second-order sliding mode control has been proposed. It allows to evaluate, under optimal criteria as minimum sum of muscle efforts, muscle activations in order to move an humanequivalent arm under a prescribed planar motion. From the control point-of-view, simulations have shown that the closedloop system is robust versus mass variations. Furthermore, in spite of difference between control synthesis model and simulation model (muscle activation dynamics are neglected for control design), tracking error converges towards zero in an arbitrary fixed time t c . This latter time is limited only by physical constraints [START_REF] Zajac | Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control[END_REF] as muscle activation u i (t) must be such that (0 ≤ u i (t) ≤ 1). Futher works consist in applying a such strategy on a biped model in order to simulate human gait and to obtain medical staffs simulators usable for pathology diagnostic or operation consequences.

Gravity

4 Fig. 1 .

 41 Fig. 1. Schematic representation of the human arm. The arm is moved by 6 muscles: 4 monoarticular ones (m 1 , m 2 , m 3 m 4 ) and 2 biarticular ones (m 5 ,m 6 ).

Fig. 2 .

 2 Fig. 2. Elements of human arm musclen driven simulation.U = [u 1 • • • u 6 ]T is the muscular excitations vector (and then the control input vector). A = [a 1 • • • a 6 ] T is the muscular activations vector. F = [Fa 1 • • • Fa 6 ] T is the active muscular forces cvector. L = [l 1 • • • l 6 ]T is the muscular fiber length vector. Γ = [Γ 1 Γ 2 ] T is the joint torques vector. q = [q 1 q 2 ] T is the joint positions vector.

) = 0 .

 0 Parameters t c , λ 01 and λ 02 are determined by minimization of the performance index (14) after setting k 1 = 20 and k 2 = 30. The obtained values are λ 0i = 17.81, t c = 4.8, a 3,i = 0.0142, a 2,i = -0.10230, a 1,i = 0, a 0,i = 0.7854 (i = 1, 2).

Fig. 4 .

 4 Fig. 4. Second order sliding mode control. Joint positions (deg) versus time (sec): trajectory reference (solid line), nominal case (dashed line) and uncertain case (+50% mass error) (dotted line).

q

  ri (t), i = {1, 2}) η i = k p s i + k d ṡi . Parameters k p and k d are determined by minimization of the performance index (14) which gives k p = 18.2070, k d = 0.1774.

Fig. 5 .

 5 Fig. 5.Second order sliding mode control. Torques (N • m) versus time (sec): nominal case (dashed line) and uncertain case (+50% mass error) (dotted line).

Fig. 6 .

 6 Fig.6.PD control. Joint positions (deg) versus time (sec): trajectory reference (solid line), nominal case (dashed line) and uncertain case (+50% mass error) (dotted line).

Fig. 7 .

 7 Fig. 7. PD control. Torques (N •m) versus time (sec): nominal case (dashed line) and uncertain case (+50% mass error) (dotted line).

Equations (5) are linear in a i . Then, it is easy to find an analytic solution for a i and to show that it is a BIBO (bounded input bounded output) stable system.

ACKNOWLEDGMENT

Carlos Rengifo would like to thank Cauca University (Colombia) for his financial support.


Yang Tian
email: yang.tian@inria.fr

Thierry Floquet
email: thierry.floquet@ec-lille.fr

Lotfi Belkoura
email: lotfi.belkoura@univ-lille1.fr

Wilfrid Perruquetti
email: wilfrid.perruquetti@ec-lille.fr

Algebraic switching time identification for a class of linear hybrid systems

Keywords: Linear systems, switched systems, switching time estimation, distribution theory

In this paper, a method for the finite time estimation of the switching times in linear switched systems is proposed. The approach is based on algebraic tools and distribution theory. Switching time estimates are given by explicit algebraic formulae that can be implemented in a straightforward manner using standard tools from computational mathematics. Simulations illustrate the proposed techniques.

Introduction

Many systems encountered in practice exhibit switchings between several subsystems, inherently by nature, such as when a physical plant has the capability of undergoing several operational modes, or as a result of the controller design, such as in switching supervisory control.

Switched systems can be seen as higher-level abstractions of hybrid systems, obtained by neglecting the details of the discrete behavior. A switched system is composed of a family of dynamical (linear or nonlinear) subsystems and a rule, called the switching law, that orchestrates the switching between them (see [START_REF] Liberzon | Switching in systems and control[END_REF][START_REF] Sun | Switched linear systems: Control and design[END_REF] for surveys). In the recent years, there has been an increasing interest in the control problems of switched systems due to their significance from both a theoretical and practical point of view. Important results for switched systems have been achieved for problems including stability [START_REF] Agrachev | Lie-algebraic stability criteria for switched systems[END_REF][START_REF] Boscain | Stability of planar switched systems: The linear single input case[END_REF][START_REF] Branicky | Multiple Lyapunov functions and other analysis tools for switched and hybrid systems[END_REF][START_REF] Mancilla-Aguilar | A converse lyapunov theorem for nonlinear switched systems[END_REF][START_REF] Vu | Common Lyapunov functions for families of commuting nonlinear systems[END_REF], stabilization [START_REF] Moulay | Stabilization of nonlinear switched systems using control lyapunov functions[END_REF][START_REF] Persis | Supervisory control with statedependent dwell-time logic and constraints[END_REF][START_REF] Wicks | Switched controller synthesis for the quadratic stabilization of a pair of unstable linear systems[END_REF][START_REF] Zhai | Quadratic stabilizability of switched linear systems with polytopic uncertainties[END_REF], tracking [START_REF] Bourdais | Towards a model-free output tracking of switched nonlinear systems[END_REF] or controllability [START_REF] Sun | Controllability and reachability criteria for switched linear control[END_REF][START_REF] Xie | Controllability of switched linear systems[END_REF].

Observability and state estimation is a key problem for such systems because both the active mode and the continuous state have to be estimated and this during a finite time interval. The notion of state estimation for switched systems was introduced in [START_REF] Ackerson | On state estimation in switching environnements[END_REF]. Observability notions for some classes of hybrid systems such as switched linear systems has been discussed and characterized in recent works such as [START_REF] Babaali | Observability of switched linear systems in continuous time[END_REF], [START_REF] Santis | On observability and detectability of continuous-time linear switching systems[END_REF], [START_REF] Vidal | Observability of linear hybrid systems[END_REF].

The problem is to recover from available measurements the state of the system and/or the switching signal, and eventually the switching time. Different observation and identification methods have been performed during the last years ( [START_REF] Alessandri | Design of luenberger observers for a class of hybrid linear systems[END_REF][START_REF] Balluchi | Design of observers for hybrid systems[END_REF][START_REF] Basseville | Detection of abrupt changes -theory and application[END_REF][START_REF] Dasgupta | A polynomial-time algorithm for checking equivalence under certain semiring congruences motivated by the state-space isomorphism problem for hybrid systems[END_REF][START_REF] Davila | Finite-time observation of the continuous and discrete state for a class of nonlinear switched dynamics[END_REF][START_REF] Domlan | Switching systems: Active mode recognition, identification of the switching law[END_REF][START_REF] Fantuzzi | Identification of piecewise affine models in noisy environment[END_REF][START_REF] Hashambhoy | Recursive identification of switched arx models with unknown number of models and unknown orders[END_REF][START_REF] Heemels | Equivalence of hybrid dynamical models[END_REF][START_REF] Huang | Identification of hybrid linear timeinvariant systems via subspace embedding and segmentation (ses)[END_REF][START_REF] Orani | Robust reconstruction of the discrete state for a class of nonlinear uncertain switched systems[END_REF][START_REF] Paoletti | Identification of hybrid systems: a tutorial[END_REF][START_REF] Saadaoui | Exact differentiation and sliding mode observers for switched lagrangian systems[END_REF]). Usually, the hybrid observer consists of two parts: an index estimator of the current active sub-model and a continuous observer that estimates, asymptotically in most cases, the continuous state of the hybrid system.

The aim of this paper is to estimate in "real-time" the switching time sequence of some class of switched linear systems with the knowledge (full or partial) of the continuous state only. The possibility to have finite time estimate for that kind of systems is clearly important, not to say crucial. The approach considered here takes root in recent works developed in [START_REF] Fliess | An algebraic framework for linear identification[END_REF] for parameter identification of linear time-invariant systems. This method is based on algebraic tools (differential algebra, module theory and operational calculus) and results in finite time estimates given by explicit algebraic formula that can be implemented in a straightforward manner using standard tools from computational mathematics. Those results have been extended to the problems of closed-loop parametric estimation for continuous-time linear systems in [START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF], parametric estimation of systems with delayed and structured entries in [START_REF] Belkoura | Parameters estimation of systems with delayed and structured entries[END_REF], state estimation of linear systems with time-varying parameters in [START_REF] Tian | Fast state estimation in linear timevariant systems: an algebraic approach[END_REF] or with delays in [START_REF] Belkoura | On-line identification of systems with delayed inputs[END_REF], fault diagnosis in [START_REF] Fliess | Robust residual generation for linear fault diagnosis: an algebraic setting with examples[END_REF], nonlinear systems with unknown inputs in [START_REF] Barbot | An algebraic framework for the design of nonlinear observers with unknown inputs[END_REF] or some class of linear infinite dimensional systems governed by partial differential equations in [START_REF] Rudolph | An algebraic approach to parameter identification in linear infinite dimensional systems[END_REF]. This approach was also applied in [START_REF] Fliess | Real-time estimation for switched linear systems[END_REF] for the estimation of the index corresponding to the current active subsystem, and the state variable of this subsystem. In this paper, the algebraic approach is extended to the problem of the finite time identification of the switching occurrences. Using some properties of the distribution theory, the switching time estimation is given by an explicit formula, as a function of the integral of the output, in order to attenuate the influence of measurement noises.

This paper is organized as follows: Section 2 gives the problem formulation. The main result is derived in Section 3. First, the switching time identification of one commutation between two subsystems is analyzed. Then, the result is extended to the case of commutations among an arbitrary number of subsystems. The last section provides a simulation that illustrates the proposed approach.

PROBLEM STATEMENT

In this work, we study a class of linear hybrid systems called switched linear systems (see [START_REF] Vidal | Observability of linear hybrid systems[END_REF] for a definition), i.e. systems whose evolution is determined by a collection of linear models (Q subsystems) with continuous state x ∈ R n connected by switches among a number of discrete states q ∈ I Q {1, . . . , Q}, modeled by linear ordinary differential equations of the following form:

ẋ = A q x, (1)
where A q ∈ R n×n are constant matrices. For the sake of convenience and without loss of generality, it is assumed that at each time t ∈ R only one discrete event can act on the system.

The objective of this paper is to estimate on line the switching times of the system (1) using the dynamics of the continuous state only. The method that will be developed will involve high order time derivatives of x. Because of the presence of non smooth dynamics, derivation has to be understood in the distribution sense.

Distribution Framework

We recall here some standard definitions and results from distribution theory developed in [START_REF] Schwartz | Théorie des distributions[END_REF], and fix the notations to be used in the sequel. The space of C ∞ -functions having compact support in an open subset Ω of R is denoted by D(Ω), and D ′ (Ω) is the space of distributions on Ω, i.e., the space of continuous linear functionals on D(Ω).

When concentrated at a point {τ}, the Dirac distribution δ (t -τ) is written δ τ . Functions are considered through the distributions they define and are therefore indefinitely differentiable. Hence, if y is a continuous function except at a point a with a finite jump σ a , its derivative writes

ẏ = dy/dt + σ a δ a , (2)
where dy/dt is the distribution stemming from the usual derivative of y.

A distribution is said to be of order r if it acts continuously on C r -functions but not on C r-1 -functions. Measures and functions are of order 0. By virtue of Schwartz Theorem (see [START_REF] Schwartz | Théorie des distributions[END_REF]), α × T = 0 for any smooth function α whose derivatives of appropriate order vanish on the support of a given distribution T ∈ D ′ (Ω). In particular, one has for the Dirac distribution 1 :

f (t) • δ t i = f (t i ) • δ t i , ∀ f ∈ C ∞ (R) if f is continuous at t i .
(3)

Switching time identification

The goal of this work is to obtain an algebraic relation of the measured variables which involves only known parameters in order to explicitly obtain the unknown switching times. To do so, first, a time-varying state transformation is used to get an intermediate differential algebraic relation which parameters are not depending anymore on the switching times but on a Dirac distribution at the switching instant. Then, an annihilating algebraic manipulation based on (3) is provided to get the desired differential algebraic relation where the unknown switching times explicitly appear. Finally, using iterative integral operators, the final expression of the switching instant is obtained in terms of time integrals of the state variables.

The case of two subsystems and one commutation

In this case, q ∈ {1, 2}. Assume that the system (1) switches from a subsystem to the other one at time t c . The aim is to get a finite time estimation of the switching time t c . For the sake of simplicity, the method is first detailed assuming that all the state is measured (some hints for the case of partial state measurement are given in Section 3.4).

The dynamical behavior of the system can be written as follows:

ẋ = Γ(t)x Γ(t) ∈ {A 1 , A 2 } . (4)
Under the change of variable z = e Gt x, where G will be defined later, the system ( 4) is transformed into:

ż = M(t)z with M(t) = G + e Gt Γ(t)e -Gt . (5)
Then, the matrix G is chosen such that:

M 1 (t) + M 2 (t) = 0 (6)
where

M i (t) = G + e Gt A i e -Gt , i = 1, 2
Since G and e Gt commute, equation ( 6) implies that

G = - A 1 + A 2 2 (7)
and

ż = σ (t)M 1 (t) z (8)
with σ (t) ∈ {-1, 1} and, without loss of generality,

σ (0+) = 1 if Γ(0+) = A 1 .

Remark 1. When the matrices A 1 and

A 2 commute, M 1 (t) = -M 2 (t) = A 1 -A 2 2
are constant. This case was treated in [START_REF] Tian | Algebraic switching time identification for a class of linear hybrid systems[END_REF].

The determinant ∆ 1 (t) of M 1 (t) satisfies the following property:

2 n ∆ 1 (t) = det(2M 1 (t)) = det(M 1 -M 2 ) = det(A 1 -A 2 ) = cte. ( 9
)
Assume henceforth that the matrix (A 1 -A 2 ) is full rank. Then, M 1 (t) is invertible and one can define W 1 (t) := Ad(M 1 (t)) = ∆ 1 (t)M -1 1 (t). Equation ( 8) becomes:

W 1 (t)ż = σ (t) ∆ 1 (t) z. ( 10
)
Note that every term in [START_REF] Boscain | Stability of planar switched systems: The linear single input case[END_REF] is known but the evolution of σ (t). Since σ (t) is a constant outside the set {t c }, the time derivative of the product σ (t)g(t) is well-defined as soon as g(t) is a smooth function for t = t c . Denoting σ c = σ (t c +) -σ (t c -) = ±2, one has:

σ (t) = σ c δ t c (σ (t)g(t)) = σ (t)g(t) + σ (t) ġ(t) = σ c g(t c )δ t c + σ (t) ġ(t). ( 11
)
By derivation of [START_REF] Boscain | Stability of planar switched systems: The linear single input case[END_REF], one obtains:

Ẇ1 (t)ż + W 1 (t)z = σ (t) ∆ 1 (t) z + σ (t) ∆ 1 (t) ż = σ (t) ∆ 1 (t) z + ∆ 1 (t) M 1 (t) z W 1 (t)z + Ẇ1 (t)ż -∆ 1 (t) M 1 (t) z = σ (t) ∆ 1 (t) z. ( 12
)
The equation ( 12),using property [START_REF] Bourdais | Towards a model-free output tracking of switched nonlinear systems[END_REF], can be expressed as:

2 ∑ i=0 K i (t) z (i) = γ c δ t c ( 13
)
K 0 (t) = -∆ 1 (t) M 1 (t) K 1 (t) = Ẇ1 (t) K 2 (t) = W 1 (t) γ c = σ c ∆ 1 (t c )z(t c )
Thus, using the change of variables z = e Gt x, one obtains a differential systems with a left-hand side that only contains known quantities. The right-hand side involves a Dirac distribution that will be annihilated using property (3), as shown in the next section.

Explicit computation of the switching instant

Take any function with the following properties:

(i) f (t,t c )δ t c = 0, (ii) f (0,t c ) = ḟ (0,t c ) = 0.
Multiplying (13) by f (t,t c ), one obtains:

f (t,t c ) 2 ∑ i=0 K i (t) z (i) = 0. ( 14
)
Integrating ( 14) from 0 to t > t c leads to

t 0 f (τ,t c ) W 1 (τ)z + Ẇ1 (τ)ż -∆ 1 (τ) M 1 (τ)z dτ = 0. ( 15
)
Integration by parts gives:

t 0 f (τ,t c )W 1 (τ)z(τ)dτ = [ f (τ,t c )W 1 (τ)ż(τ)] τ=t τ=0 - t 0 ḟ (τ,t c )W 1 (τ) + f (τ,t c ) Ẇ1 (τ) ż(τ)dτ = f (t,t c )W 1 (t)ż(t) - t 0 ḟ (τ,t c )W 1 (τ)ż(τ)dτ - t 0 f (τ,t c ) Ẇ1 (τ)ż(τ)dτ and t 0 ḟ (τ,t c )W 1 (τ)ż(τ)dτ = ḟ (τ,t c )W 1 (τ)z(τ) τ=t τ=0 - t 0 f (τ,t c )W 1 (τ) + ḟ (τ,t c ) Ẇ1 (τ) z(τ)dτ = ḟ (t,t c )W 1 (t)z(t) - t 0 f (τ,t c )W 1 (τ)z(τ)dτ - t 0 ḟ (τ,t c ) Ẇ1 (τ)z(τ)dτ. Hence t 0 f (τ,t c ) W 1 (τ)z + Ẇ1 (τ)ż dτ = f (t,t c )W 1 (t)ż(t) -ḟ (t,t c )W 1 (t)z(t) + t 0 f (τ,t c )W 1 (τ)z(τ)dτ + t 0 ḟ (τ,t c ) Ẇ1 (τ)z(τ)dτ (16)
Using the property

t 0 t ν-1 0 • • • t 1 0 x(τ)dt ν-1 • • • dt 1 dτ = t 0 (t -τ) ν-1 (ν -1)! x(τ)dτ,
integrating one more time (15) from 0 to t > t c and using ( 16), one obtains:

f (t,t c )W 1 (t)z(t) -2 t 0 ḟ (τ,t c )W 1 (τ)z(τ)dτ - t 0 f (τ,t c ) Ẇ1 (τ)z(τ)dτ + t 0 (t -τ) f (τ,t c )W 1 (τ) + ḟ (τ,t c ) Ẇ1 (τ) -∆ 1 (τ) M 1 (τ) z(τ)dτ = 0. ( 17
)
An estimate of t c can be obtained from [START_REF] Fliess | Real-time estimation for switched linear systems[END_REF]. Take for instance the function f (t,t c ) = t 2 (tt c ) that satisfies the properties (i)-(ii) and

ḟ (t,t c ) = t(3t -2t c ), f (t,t c ) = 6t -2t c .
Using the result of ( 17), the estimate of t c is given by the following formula:

D(t, z, M 1 )t c = N(t, z, M 1 )
with

N(t, z, M 1 ) = t 3 W 1 (t)z(t) - t 0 (6τ 2 W 1 (τ) + τ 3 Ẇ1 (τ))z(τ)dτ + t 0 (t -τ)(6τW 1 (τ) + 3τ 2 Ẇ1 (τ) -∆ 1 (τ) M 1 (τ))z(τ)dτ D(t, z, M 1 ) = t 2 W 1 (t)z(t) - t 0 (4τW 1 (τ) + τ 2 Ẇ1 (τ))z(τ)dτ + t 0 (t -τ)(2W 1 (τ) + 2τ Ẇ1 (τ))z(τ)dτ (18)
Note that N(t, z, M 1 ) and D(t, z, M 1 ) are column vectors of dimension n. So

t c = N g (t, z, M 1 ) D g (t, z, M 1 ) , (19)
for any 1 ≤ g ≤ n, where N g (t, z, M 1 ) and D g (t, z, M 1 ) are the components of N(t, z, M 1 ) and D(t, z, M 1 ), respectively.

Remark 2. When the measured signal z(t) is perturbed by noise, ( 17) can be integrated one or several more times in order to reduce the noise effect. In this case, the expression of N(t, z, M 1 ) and D(t, z, M 1 ) in [START_REF] Fliess | Robust residual generation for linear fault diagnosis: an algebraic setting with examples[END_REF] becomes

N(t, z, M 1 ) = t 0 τ 3 W 1 (τ)z(τ)dτ - t 0 (t -τ)(6τ 2 W 1 (τ) + τ 3 Ẇ1 (τ))z(τ)dτ + t 0 1 2 (t -τ) 2 (6τW 1 (τ) + 3τ 2 Ẇ1 (τ) -∆ 1 (τ) M 1 (τ))z(τ)dτ D(t, z, M 1 ) = t 0 τ 2 W 1 (τ)z(τ)dτ - t 0 (t -τ)(4τW 1 (τ) + τ 2 Ẇ1 (τ))z(τ)dτ + t 0 1 2 (t -τ) 2 (2W 1 (τ) + 2τ Ẇ1 (τ))z(τ)dτ (20)
Hereafter, the method is extended to the problem of the identification of an arbitrary number of switching times between an arbitrary number of subsystems.

The case of Q subsystems and S commutations

Assume that every pair of matrices A i and A j (i, j ∈ I Q and i = j) in system ( 1) is such that (A i -A j ) is full rank. It has been seen that an estimator E i, j that computes the quantities D(t, z, M i, j ) and N(t, z, M i, j ), with M i, j = G i, j + e G i, j t A i e -G i, j t , G i, j = -

A i +A j 2
can be used to determine a switching time t i, j that occurs between the two subsystems i and j (either from mode i to mode j or from mode j to mode i).

Hence, in order to identify all the switches among Q modes, one can use

C 2 Q = Q(Q-1) 2
estimators in parallel. Then, the output signals of each estimator can be analyzed as follows to determine the occurrence of a switch and its associated mode. Indeed, assume that the system (1) is in the mode i for t ∈ [t 0 ,t i, j [ and in mode j for t ∈ [t i, j , T [, where times t 0 and T stands for other switch occurrences. Then, one has:

1. D g (t, z, M i, j ) = 0 and N g (t, z, M i, j ) = 0 for t ∈ [t 0 ,t i, j [: in practice, this property is checked on a few sampling periods before the switching instant; 2. D g (t, z, M i, j ) and N g (t, z, M i, j ) are straight lines for t ∈ [t i, j , T [: in practice, several points have to be used to detect that the slope of the line is constant.

3. for t ∈ [t i, j , T [ the ratio of

N g (t,z,M i, j ) D g (t,z,M i, j
) is constant and equal to t i, j : to detect when this ratio is constant, successive values within a boundary of a given thickness are needed. Properties 1. and 2. arise from the fact that Dirac distributions in equation ( 13) are integrated twice. Thus, the estimator that fulfills those three conditions provides the switching instant as well as the index of the subsystems between which the switch occurs. Furthermore, with the knowledge of the first active mode, the sequence of all active models can be estimated.

The implementation of the numerical algorithm which allows to identify all the switching instants is done on a sliding window as follows: function detection of the switching time t c read x(kT e ) where y ∈ R p is the output vector of available measurements. Assume it is possible to find an output relation given by:

Ẏ = Γ(t)Y + Θ δ t c , δt c , . . . , δ (l) t c , Γ(t) ∈ Ā1 , Ā2 (21)
where Y = y T , . . . , y (r) T T is a vector composed of the output and a finite number of its time derivatives, Ā1 , Ā2 are square matrices, and where Θ is a square matrix that depends on the Dirac distribution and a finite number of its time derivatives. Using the same procedure as in Section 3.1, it is possible to obtain a differential system similar to (13):

2 ∑ i=0 Ki (t) z(i) = Θ δ t c , δt c , . . . , δ ( l) t c
where z = e Ḡt Y . Then, setting ω = r z leads to:

2 ∑ i=0 Ki (t) ω (i+r) = Θ δ t c , δt c , . . . , δ ( l) t c
.

Following the lines of Section 3.2, with a suitable choice of the function f (t,t c ) (that allows for the annihilation of Θ) and a sufficient number of integrations, one can obtain an algebraic relation of the form:

D(t, ω, M1 )t c = N(t, ω, M1 ),
that provides an on-line estimation of t c , since ω is available for measurement (it depends on the output y and a finite number of its integrals).

To completely solve the problem of switching time identification of systems with partial state measurement using the proposed approach, one has to study the possibility to obtain the relation [START_REF] Hashambhoy | Recursive identification of switched arx models with unknown number of models and unknown orders[END_REF]. For this, observability and identifiability properties of the system (4) have to be investigated. This problem is out of the scope of this paper (whose aim is to introduce an algebraic methodology for fast identification of switching times) and will be studied in future work.

EXAMPLE

Let us consider the buck-boost converter shown in Figure 1, the switched model consisting of two second-order linear time invariant systems is given by: The mode 1 represents the transistor T closed and the diode D blocked (T = 1, D = 0), while the mode 2 correspond with T OFF and D ON (T = 0, D = 1). This model is well being to the class defined in the section 3.4. Note that mode 1 is unstable, since the current increases linearly with time when T is closed, while mode 2 (T open) is stable at the origin.

       ẋ = A 1 x + B = -1 RC 0 0 0 x + 0 U L ẋ = A 2 x = -1 RC -1 C 1 L 0 x ( 22
)
where x = [V i] T = [x 1 x 2 ]
One can rewrite this system as follows:

       Ẏ = 0 0 0 0 Y + U L U L δ t c Ẏ = 0 1 -1 LC -1 RC Y (23)
with Y = [y ẏ] T corresponds with the relation [START_REF] Hashambhoy | Recursive identification of switched arx models with unknown number of models and unknown orders[END_REF].

One has Ḡ = -

0 0 0 0 + 0 1 -1 LC -1 RC 2 = 0 1 1 LC 1 RC 2
to change the variable z = e Ḡt Y . Setting ω = t 0 z(τ)dτ, then ω is available for measurement (it is as a function of the output y and its integral). In order to annihilate Θ (δ t c ), here one chooses the function f (t,t c ) = t 2 (tt c ) 2 satisfying with the properties:

(i) f (t,t c )δ t c = f (t,t c ) δt c = 0, (ii) f (0,t c ) = ḟ (0,t c ) = 0.
After twice integrations one obtains:

D 1 (t, ω, M1 )t 2 c -D 2 (t, ω, M1 )t c + N(t, ω, M1 ) = 0 (24)
with

N(t, ω, M1 ) = t 4 W 1 (t)ω(t) - t 0 (8τ 3 W 1 (τ) + τ 4 Ẇ1 (τ))ω(τ)dτ + t 0 (t -τ)(12τ 2 W 1 (τ) + 4τ 3 Ẇ1 (τ) -∆ 1 (τ) M1 (τ))ω(τ)dτ D 1 (t, ω, M1 ) = t 2 W 1 (t)ω(t) - t 0 (4τW 1 (τ) + τ 2 Ẇ1 (τ))ω(τ)dτ + t 0 (t -τ)(2W 1 (τ) + 2τ Ẇ1 (τ))ω(τ)dτ D 2 (t, ω, M1 ) = 2t 3 W 1 (t)ω(t) - t 0 (12τ 2 W 1 (τ) + 2τ 3 Ẇ1 (τ))ω(τ)dτ + t 0 (t -τ)(12τW 1 (τ) + 6τ 2 Ẇ1 (τ))ω(τ)dτ
Hence, one can calculate t c explicitly from the equation ( 24). The Buck-Boost circuit has the following characteristics: U = 24(V ), L = 100(µH), C = 220(µF), R = 20(Ω). In order to get a good simulation, since the above parameters L and C are of order 10 -4 , we use a time scaling T = 1 ε t with ε = 10 -4 , that means we simulate dx dT = Ax + b instead of dx dt = 1 ε (Ax + b). The switching instants of this system are detected using the proposed algorithm with T e = 0.001s with respect to the new time scaling. The switching times are assumed to occur at: t 1 = 80(µs), t 2 = 150(µs), t 3 = 220(µs), t 4 = 280(µs) (real time) (see Fig. 3) and the first active mode is 1. Note that in the simulation, one chooses 30 sampling periods to identify the switching time, so the dwell time of the system has to be at least larger than 3(µs) with respect to the real time. In Fig. 2 is reported the behavior of the state x(t). Fig. 3 shows that the estimator scheme can detect all switching times t c accurately. Note that the methodology of this algebraic approach can be applied to solve some practical problems, such as for instance fault detection, fault tolerant control issues [START_REF] Rhee | Issues of fail-over switching for fault-tolerant ethernet implementation[END_REF].

Conclusion

In this paper, an algebraic approach for switching time estimation of a class of hybrid systems has been introduced. Using an approach based on algebraic tools and distribution theory, an explicit algorithm which computes on-line the switching time instants in a fast way has been derived for systems with full and partial state measurement. Future works are concerned with observability and identifiability issues.

Figure 1 :

 1 Figure 1: Buck-boost converter circuit.

Figure 2 Figure 3 :

 23 Figure 2: State x

The existence of the product f (t) • δ t i is ensured as soon as f is a C ∞ -function.

T e : sample time; h: number of chosen points to verify the constraints with h.T e less than the dwell time of the system; t s = 0: left bound of the integration window; m = length(x): length of the sampled state x; z i, j (kT e ) = e - A i +A j 2 x(kT e ): new state for k = h : m Reset all the components of N(t, z, M i, j )(k) and D(t, z, M i, j )(k) to zero for i = t s : k Approximate numerically the integrals as sums using the trapezoidal rule end Compute the value of N(t, z, M i, j )(k) and D(t, z, M i, j )(k) in an integration window with [START_REF] Fliess | Robust residual generation for linear fault diagnosis: an algebraic setting with examples[END_REF] if an estimator satisfies the three detection criterions update the integration window: t s ←kh; print this estimator and the active model. else break; end print the switching time end

The numerical implementation of the above obtained algorithm has also to address the following points:

• integration window: when none of the

2 estimator detect the commutation, the algorithm will calculate the value of N(t, z, M i, j )(k), D(t, z, M i, j )(k) in a new integration window (k ←k + 1) whose width is one sample time larger than the old one (the lower bound of the integration is not changed while the upper bound is increased of one sampling), until a switching instant is detected with one of the estimator.

• update: when t c is detected, one has to update the scheme by resetting the integral computation and replacing the lower bound of the integration with t c in order to identify the next switching time.

A hint for an extension to systems with partial state measurement

Consider the following system:


Yang Tian

Thierry Floquet

Lotfi Belkoura

Wilfrid Perruquetti

Switching time estimation for linear switched systems: an algebraic approach

Keywords: Linear systems, switched systems, hybrid systems, switching time estimation, distribution theory

This paper aims at estimating the switching time for linear switched systems, i.e. the time instant when a sub-model is switched on while another one is switched off. Assuming that the state of the active sub-models is known, a distribution point of view is adopted to get a real-time estimation of these switching times. Real-time means that an explicit algorithm computes on-line these time instants in a fast and efficient way. Simulations illustrate the proposed techniques which is easily implementable.

I. INTRODUCTION

Many systems encountered in practice exhibit switchings between several subsystems, inherently by nature, such as when a physical plant has the capability of undergoing several operational modes, or as a result of the controller design, such as in switching supervisory control.

Switched systems may be viewed as higher-level abstractions of hybrid systems, obtained by neglecting the details of the discrete behavior. Roughly speaking, a switched system is composed of a family of dynamical (linear or nonlinear) subsystems and a rule, called the switching law, that orchestrates the switching between them. In the recent years, there has been increasing interest in the control problems of switched systems due to their significance from both a theoretical and practical point of view. Important results for switched systems have been achieved for problems including stability [START_REF] Agrachev | Lie-algebraic stability criteria for switched systems[END_REF], [START_REF] Boscain | Stability of planar switched systems: The linear single input case[END_REF], [START_REF] Vu | Common Lyapunov functions for families of commuting nonlinear systems[END_REF], controllability [START_REF] Sun | Controllability and reachability criteria for switched linear control[END_REF], [START_REF] Xie | Controllability of switched linear systems[END_REF], stabilization [START_REF] Moulay | Stabilization of nonlinear switched systems using control lyapunov functions[END_REF], [START_REF] Zhai | Quadratic stabilizability of switched linear systems with polytopic uncertainties[END_REF], tracking [START_REF] Bourdais | Towards a model-free output tracking of switched nonlinear systems[END_REF], . . . . See also [START_REF] Liberzon | Switching in systems and control[END_REF], [START_REF] Sun | Switched linear systems: Control and design[END_REF] for surveys.

Observability and state estimation is a key problem for such systems, where discrete and continuous parts are mixed. The vector-valued output is a function of the continuous state. The observation problem is to determine the continuous state and the active ODE driving the current continuous state evolution. In [START_REF] Ackerson | On state estimation in switching environnements[END_REF], the notion of state estimation for switched systems is introduced. Observability notions for some classes of hybrid systems such as switched linear systems has been discussed and characterized in recent works, see e.g. [START_REF] Babaali | Observability of switched linear systems in continuous time[END_REF], [START_REF] Santis | On observability and detectability of continuous-time linear switching systems[END_REF], [START_REF] Vidal | Observability of linear hybrid systems[END_REF]. The problem is to recover from available measurements the state of the system and/or the switching signal, Y. and eventually the switching time. Different observation and identification methods have been proposed during the last years ([3], [START_REF] Balluchi | Design of observers for hybrid systems[END_REF], [START_REF] Fantuzzi | Identification of piecewise affine models in noisy environment[END_REF], [START_REF] Hashambhoy | Recursive identification of switched arx models with unknown number of models and unknown orders[END_REF], [START_REF] Huang | Identification of hybrid linear timeinvariant systems via subspace embedding and segmentation[END_REF], [START_REF] Paoletti | Identification of hybrid systems: a tutorial[END_REF], [START_REF] Saadaoui | Exact differentiation and sliding mode observers for switched lagrangian systems[END_REF]). In most of the cases, the hybrid observer consists of two parts: an index estimator of the current active sub-model and a continuous observer that estimates, asymptotically in most cases, the continuous state of the hybrid system.

The aim of this paper is to estimate in "real-time" the switching time sequence with the knowledge of the continuous state of some class of switched linear systems. The possibility to have finite time estimate for that kind of systems is clearly important, not to say crucial. The approach considered here is based on an idea recently published by M. Fliess and H.-S. Ramírez [START_REF] Fliess | An algebraic framework for linear identification[END_REF] for parameter identification of linear time-invariant systems. This method is based on tools which are of algebraic flavor: differential algebra, module theory and operational calculus. The estimators are non asymptotic: solutions are provided in finite time by explicit algebraic formulae and result in relatively simple and fast algorithms, with straightforward implementation using standard tools from computational mathematics. Those results have been extended to the problems of closed-loop parametric estimation [START_REF] Fliess | Closed-loop parametric identification for continuous-time linear systems via new algebraic techniques[END_REF], state estimation of linear systems with time-varying parameters [START_REF] Tian | Fast state estimation in linear time-variant systems: an algebraic approach[END_REF] or with delays [START_REF] Belkoura | On-line identification of systems with delayed inputs[END_REF], fault diagnosis [START_REF] Fliess | Robust residual generation for linear fault diagnosis: an algebraic setting with examples[END_REF], nonlinear systems with unknown inputs [START_REF] Barbot | An algebraic framework for the design of nonlinear observers with unknown inputs[END_REF] or some class of linear infinite dimensional systems governed by partial differential equations [START_REF] Rudolph | An algebraic approach to parameter identification in linear infinite dimensional systems[END_REF]. This approach was also applied in [START_REF] Fliess | Real-time estimation for switched linear systems[END_REF] for the estimation of the index corresponding to the current active subsystem, and the state variable of this subsystem, via quite robust methods with respect to corrupting noises. Here, explicit computation of the switching times is obtained: the switching time estimation is given by an explicit formula, as a function of the integral of the output, in order to attenuate the influence of measurement noises. The proposed method exhibits the following features:

• the switching times can be efficiently identified on-line and in finite time, • computations can be carried out by a computer and in a very fast manner. The rest of the paper is organized as follows: the problem formulation is given in Section II. Then, the main result is derived in Section III last section provides a simulation example in order to highlight the efficiency of the proposed approach.

II. PROBLEM STATEMENT

Consider the class of linear switched system made of subsystems modeled linear ordinary differential equations (LODE) of the following form:

ẋ = A q x, (1)
where q ∈ I Q {1, . . . , Q} is the discrete state, x ∈ R n is the continuous state and A q ∈ R n×n . For the sake of convenience and without loss of generality, it is assumed that at each time t ∈ R only one discrete event can act on the system. It is also assumed that there is no state jump, i.e. that if t k is a switching instant lim

t→t - k x(t) = lim t→t + k x(t)
The switching function driving the change of mode is defined by:

σ (t) : R → I Q x → σ (t) (2)
where σ (t) ∈ I Q corresponds to the index associated with the current active LODE. The general problem to be solved is: "Is it possible to get in real time (when they occur) the switching times t i using the available past and current measurements?". In this paper, the study is restricted to the case when Q = 2 and when the whole continuous states are measured. In that case, system (1) can be rewritten as follows:

ẋ = Γ(t)x, (3)
Γ(t) = A 1 H(t) + A 2 (1 -H(t)) ( 4
)
where A 1 , A 2 ∈ R n×n are constant matrices. Without loss of generality 1 , assume that q = 1 within some given interval t ∈ [t i-1 ,t i [ and q = 2 on t ∈ [t i ,t i+1 [ . Then, the function H is given by:

H(t) = 1, t i-1 ≤ t < t i 0, t i ≤ t < t i+1
Suppose that the previous switching time was found: thus t i-1 is known. The objective is to have a finite time estimation of the switching time t i when the system commutes from the mode q = 1 to the mode q = 2. This will be done by using higher order time derivatives of x. Because of the presence of non smooth dynamics, derivation has to be understood in the distribution sense.

Distribution Framework

We recall here some standard definitions and results from distribution theory [START_REF] Schwartz | Théorie des distributions[END_REF], and fix the notations to be used in the sequel. The space of C ∞ -functions having compact support in an open subset Ω of R is denoted by D(Ω), and D ′ (Ω) is the space of distributions on Ω, i.e., the space of continuous linear functionals on D(Ω).

When concentrated at a point {τ}, the Dirac distribution δ (t -τ) is written δ τ .

Functions are considered through the distributions they define and are therefore indefinitely differentiable. Hence, if y is a continuous function except at a point a with a finite jump σ a , its derivative writes ẏ = dy/dt + σ a δ a ,

where dy/dt is the distribution stemming from the usual derivative of y.

A distribution is said to be of order r if it acts continuously on C r -functions but not on C r-1 -functions. Measures and functions are of order 0. By virtue of Schwartz Theorem [START_REF] Schwartz | Théorie des distributions[END_REF], α × T = 0 for any smooth function α whose derivatives of appropriate order vanish on the support of a given distribution T ∈ D ′ (Ω). In particular, one has for the Dirac distribution 2 :

f (t) • δ t i = 0, ∀ f ∈ C ∞ (R) with f (t i ) = 0 f (t) • δ t i = f (t i ) • δ t i , ∀ f ∈ C ∞ (R) with f (t i ) = 0 III. MAIN RESULT

A. Change of variable

Assume that A 1 and A 2 are two square commuting matrices, and consider the following change of variable:

z = e Gt x ( 6
)
where G = -A 1 +A 2 2 . One has: ż = Ge Gt x + e Gt ẋ = Gz + e Gt Γ(t)x = (G + e Gt Γ(t)e -Gt )z which leads to ż = M(t)z [START_REF] Belkoura | On-line identification of systems with delayed inputs[END_REF] where

M(t) = (G + e Gt A 1 e -Gt )H(t) + (G + e Gt A 2 e -Gt )(1 -H(t))
Since A 1 and A 2 are two square commuting matrices, one gets:

M(t) = (G + A 1 )H(t) + (G + A 2 )(1 -H(t)) = A 1 -A 2 2 H(t) + A 2 -A 1 2 (1 -H(t)) = -AH(t) + A(1 -H(t)) with A = A 2 -A 1 2 .
The choice of the change of coordinates (6) was motivated by the fact that the subsequent computations involve M 2 = A 2 which is a known quantity independent of the switching time.

Differentiating [START_REF] Belkoura | On-line identification of systems with delayed inputs[END_REF] , one has: 2 The existence of the product f (t) • δ t i is ensured as soon as f is a C ∞ -function which can be rewritten in the following compact form:

z -A 2 z = -2Aδ t i-1 z(t i-1 ), t i-1 ≤ t < t i z -A 2 z = -2Aδ t i-1 z(t i-1 ) + 2Aδ t i z(t i ), t i ≤ t < t i+1
z -A 2 z = i ∑ g=i-1 γ g δ t g z(t g ), t i-1 ≤ t < t i+1 (8)
where γ g is obviously defined. Note that a similar expression is obtained if the initial system is q = 2, but with the γ g having opposite signum. As it will be seen hereafter, this has no influence because the right hand side of (8) will be canceled by suitable algebraic manipulations. Thus the proposed method does not need the knowledge of the current mode.

B. Explicit computation of the switching instant

Recall that (tt i )δ t i = 0 and multiply (8) by

f (t,t i-1 ,t i ) = (t -t i-1 ) 2 (t -t i ). (9)
One obtains:

(t -t i-1 ) 2 (t -t i )(z -A 2 z) = 0
Integrating from t i-1 to t > t i leads to

t t i-1 (τ -t i-1 ) 2 (τ -t i )(z(τ) -A 2 z)dτ = 0 ( 10
)
The function ( 9) has also been chosen in such a way that f

(t i-1 ,t i-1 ,t i ) = f ′ (t i-1 ,t i-1 ,t i ) = 0.
Then, integration by parts leads to:

t t i-1 f (τ,t i-1 ,t i )z(τ)dτ = [ f (τ,t i-1 ,t i )ż] τ=t τ=t i-1 - t t i-1 f ′ (τ,t i-1 ,t i )ż(τ)dτ = f (t,t i-1 ,t i )ż - t t i-1 f ′ (τ,t i-1 ,t i )ż(τ)dτ and t t i-1 f ′ (τ,t i-1 ,t i )ż(τ)dτ = f ′ (τ,t i-1 ,t i )z τ=t τ=t i-1 - t t i-1 f ′′ (τ,t i-1 ,t i )z(τ)dτ = f ′ (t,t i-1 ,t i )z - t t i-1 f ′′ (τ,t i-1 ,t i )z(τ)dτ Hence t t i-1 f (τ,t i-1 ,t i )z(τ)dτ = f (t,t i-1 ,t i )ż -f ′ (t,t i-1 ,t i )z + t t i-1 f ′′ (τ,t i-1 ,t i )z(τ)dτ (11)
Taking one more time the integral from t i-1 to t > t i leads to the following relation:

t t i-1 [ f (τ,t i-1 ,t i )ż(τ) -f ′ (τ,t i-1 ,t i-1 ,t i )z(τ)]dτ + t t i-1 τ 1 t i-1 ( f ′′ (τ,t i-1 ,t i ) -A 2 f (τ,t i-1 ,t i ))z(τ)dτdτ 1 = 0 (12)
Using the properties

t 0 t ν-1 0 ••• t 1 0 x(τ)dt ν-1 ••• dt 1 dτ = t 0 (t -τ) ν-1 (ν -1)! x(τ)dτ t-t i 0 t ν-1 -t i 0 ••• t 1 -t i 0 x(τ)dt ν-1 ••• dt 1 dτ = t-t i 0 (t -τ) ν-1 (ν -1)! x(τ)dτ yields t t i t ν-1 t i • • • t 1 t i x(τ)dt ν-1 • • • dt 1 dτ = t t i (t -τ) ν-1 (ν -1)! x(τ)dτ.
Thus, one has:

f (t,t i-1 ,t i )z -2 t t i-1 f ′ (τ,t i-1 ,t i )z(τ)dτ + t t i-1 (t -τ)( f ′′ (τ,t i-1 ,t i ) -A 2 f (τ,t i-1 ,t i ))zdτ = 0 (13) Since f (t,t i-1 ,t i ) = (t -t i-1 ) 2 (t -t i ) = t(t -t i-1 ) 2 -(t -t i-1 ) 2 t i f ′ (t,t i-1 ,t i ) = (t -t i-1 ) 2 + 2(t -t i-1 )(t -t i ) f ′′ (t,t i-1 ,t i ) = 4(t -t i-1 ) + 2(t -t i )
the estimate of t i is given by the following formula:

D(t,t i-1 , z)t i = N(t,t i-1 , z) with N(t,t i-1 ,z) = t(t -t i-1 ) 2 z(t) - t t i-1 2(τ -t i-1 ) 2 + 4τ(τ -t i-1 ) z(τ)dτ + t t i-1 (t -τ) 4(τ -t i-1 )I + 2τI -A 2 τ(τ -t i-1 ) 2 z(τ)dτ D(t,t i-1 ,z) = (t -t i-1 ) 2 z(t) - t t i-1 4(τ -t i-1 )z(τ)dτ + t t i-1 (t -τ)(2I -A 2 (τ -t i-1 ) 2 )z(τ)dτ (14)
where I ∈ R n×n is the unit matrix. Note that D(t,t i-1 , z) and N(t,t i-1 , z) are column vectors of dimension n. So

t i = N j (t,t i-1 , z) D j (t,t i-1 , z) ,
for any 1 ≤ j ≤ n, where N j (t,t i-1 , z) and D j (t,t i-1 , z) are the components of N(t,t i-1 , z) and D(t,t i-1 , z), respectively. In the above expression, t i-1 has been estimated with a similar procedure and is known. The first switching time t 1 is obtained with t 0 = 0 and the following relation:

D(t, z)t 1 = N(t, z) N(t, z) = t 3 z(t) - t 0 6τ 2 z(τ)dτ + t 0 (t -τ)(6τI -A 2 τ 3 )z(τ)dτ D(t, z) = t 2 z(t) - t 0 4τz(τ)dτ + t 0 (t -τ)(2I -A 2 τ 2 )z(τ)dτ (15)
Thus, one can obtain the exact formula for the estimation of the switching times t i . Remark: one could also choose any function f (t,t i-1 ,t i ) with the following properties:

(i)

f (t,t i-1 ,t i )δ t i = 0, (ii) f (t i-1 ,t i-1 ,t i ) = ḟ (t i-1 ,t i-1 ,t i ) = 0.

C. Numerical implementation of the Integration Window

The numerical implementation of the above obtained result has to address the following problems: 1) sampling: for computation purpose the continuous signal is sampled (the sampling period is T e and quantized (not taken into account here)). It is clear that the choice of T e is of importance and should be less than the minimum dwell time (i.e. the minimum difference between two successive switching instants).

2) numerical integration: for the implementation of the algorithm, since it involves integrals (see ( 14)), one should choose an efficient numerical algorithm in order to compute the integrals. Here, a trapezoidal rule was chosen.

3) singularity: t i is given by ( 14). As a result, a singularity (0 divided by 0) occurs at the beginning of the computation. In order to overcome this problem, one has to wait a very small amount of time before the computation is taken as valid. In practice, two samplings are sufficient.

4) initialization: let us recall that this algorithm needs to be initialized with a first computation of the first switching time t 1 . For this, one has to use (15) instead of ( 14).

5) detection: the ratio that defines t i by ( 14) is fluctuating when t ∈ [t i-1 ,t i [ and constant as soon as t > t i : to detect when this ratio is constant, successive values within a boundary of a given thickness are needed. In practice 2 or 3 successive values have to be within the interval of thickness.

6) update: when t i is detected, one has to update the scheme by resetting the integral computation and changing the lower bound of the integration in order to compute the next switching time (for this in the integral approximation.

The advantage of the implementation of the algorithm with an integration window is that the estimate can be done at each integration step.

IV. EXAMPLE

Consider a linear switched system consisting of two second-order linear time invariant systems given by:

A 1 = 1 0 0 1 A 2 = -3 0 5 0
The switching instants of this hybrid system are detected using the algorithm proposed in section III with T e = 0.001.

The switching times are assumed to occur at: t 1 = 0.2(s), t 2 = 0.8(s), t 3 = 0.85(s), t 4 = 1.8(s) and t 5 = 3(s) and the first active model is subsystem 1.Since one needs three samplings to identify the switching time, the dwell time has to be at least larger than 0.003(s).

Figure 1 shows that the estimator performs well. Note that t 3 occurs shortly (0.05(s)) after t 2 but that the algorithm can detect it accurately.

Figure 2 gives the behavior of the state z.

V. CONCLUSION AND PERSPECTIVES

In this paper, an algebraic approach for switching time estimation of hybrid systems has been introduced. An explicit algorithm which computes on-line the switching time instants in a fast way has been derived. In future works, this approach will be extended to more general cases when the switching system is composed of more than two subsystems (using a sufficient number of estimators in parallel) and with only partial state measurement outputs. Such techniques could also be used to estimate on-line and in real time both the switching instants and the switching function.

Fig. 1 .

 1 Fig. 1. Estimation of the switching times.

Fig. 2 .

 2 Fig. 2. State z

  Tian, T. Floquet, and W. Perruquetti are with LAGIS (CNRS, UMR 8146) & ALIEN Project (INRIA Lille-Nord Europe), École Centrale de Lille, BP 48, Cité Scientifique, 59650

Villeneuve d'Ascq, France. E-mail: yang.tian@inria.fr, thierry.floquet@ec-lille.fr, wilfrid.perruquetti@ec-lille.fr L. Belkoura is with LAGIS (CNRS, UMR 8146) & ALIEN Project (INRIA Lille-Nord Europe), Université des Sciences et Technologies de Lille, 59650 Villeneuve d'Ascq, France. E-mail: lotfi.belkoura@univ-lille1.fr

This will be explained later (see subsection III-A).


Philippe Ben-Abdallah

Karl Joulain
email: karl.joulain@univ-poitiers.fr

Fundamental limits for non-contact transfers between two bodies

Keywords: number(s) : 44.40.+a, 42.50.Lc, 12.20.-m, 68.35.Af

We investigate energy and momentum non-contact exchanges between two arbitrary flat media separated by a gap. This problem is revisited as a transmission problem of individual system eigenmodes weighted by a transmission probability obtained either from fluctuational electrodynamics or quantum field theory. An upper limit for energy and momentum flux is derived using a general variational approach. The corresponding optimal reflectivity coefficients are given both for identical and different media in interaction.

2

Two arbitrary media in relative motion or at rest and separated by a gap continually exchange in permanence energy and momentum [1] throughout the thermally and quantum fluctuating electromagnetic field they radiate in their surrounding. At long separation distance compared to the Wien wavelength (

) T k /( c B T     
), energy radiative transfer is maximal when both media behave like blackbodies [START_REF] Planck | The Theory of Heat Radiation[END_REF]. In this situation, the transfer is driven by the famous Stefan-Boltzmann law and exchanges only depends on the difference of media temperatures power four. When two bodies are in relative motion, the momentum exchanges through Doppler-shifted photons give rise to a van der Waals interaction [START_REF] Dzyaloshinskii | [END_REF][4][5] also called the van der Waals friction stress which is opposed to the motion. At long separation distances, this transfer is maximal when both media are perfect absorbers [5]. At subwavelength scale, the situation radically changes for both energy and momentum exchanges. Indeed, in this case, the presence of evanescent modes gives rise to wave phenomena such as tunnelling and interferences which drastically affect the transfers. Hence, energy and momentum exchanges can respectively exceed by several orders of magnitude the blackbody limit [6][7][8][9][10] for heat transfer and the van der Waals friction intensity between two perfect absorbers [4,5,[11][12][13]. However, these results raised new questions. Are there fundamental limits for these exaltation mechanisms for energy [14] and momentum transfer between two bodies ? If they exist, what are these ultimate limits [15] and what are the media required to reach such values? In this paper we use calculus of variations principles in order to bring a general answer to these questions. In addition, we introduce a general Landauer-type formulation [16] for the transfer, to describe non-contact exchanges as a simple transmission problem. Relating this formulation to the results predicted by the fluctuational electrodynamic [START_REF] Rytov | Theory of Electrical Fluctuation and Thermal Radiation (Academy of Science of USSR[END_REF] and the quantum field [18] theory we derive energy and momentum transmission probabilities for both radiative and non-radiative photons. Moreover we derive and give a physical interpretation for the optimal conditions maximizing exchanges between two identical or distinct bodies. Finally, we examine near-field heat transfer between two media which support surface polaritons and compare it with the fundamental limits of transfer. and layered, separated by a vacuum gap of thickness  . These resonators are reservoirs of radiative and non-radiative electromagnetic modes that are entirely defined by a couple ) , ( q  ,  being the mode angular frequency and q the parallel component of its wavevector. Each of these modes carries a quantum of energy    e and a quantum of momentum q   m (we henceforth denote generically  any of these quantities). Whatever the separation distance between the two reservoirs, a  -flux is exchanged. This exchange occur through radiative modes coupling at long separation distances and both through non-radiative and radiative modes at short distances. If these reservoirs are maintained in nonequilibrium situations at two distinct temperatures A T and B T and are animated by a relative parallel motion with a velocity r v then a certain amount of  can be transmitted from A to B throughout a coupling channel with a transmission probability 1 ) , , (

B A       q
. Then, using a Landauer-type formulation for this  -transfer [16] we have

          q q q v , B A r A B A ) , , ( ) , , T , ( f  (1)
where f is a function which depends on the nature of exchanges we are dealing with. Suppose now that the transfer occur during a time t  . Mode angular frequency are quantified as well as parallel mode if we suppose that the system is bounded by an arbitrary large square box area 2 L S 

. Then, we have

t n t     and y y x x L 2 n L 2 n e e q     , (2) where ) y
x e , (e denote two orthogonal vectors in the plane parallel to the interfaces, t n is a positive integer, x n and y n are relative integers respectively. This quantification of mode specifies the number of channels available in the (q,)-space. Examining the  exchange between A and B (resp. B and A) from this surface during a time interval t  we get , after transforming the discrete summation in (1) into an integration over a continuum

               ) , q , ( ) q , v , T , ( f ) L / 2 ( ) t / ( d 2 1 B A r A 2 0 B A  dq . ( 3
)
and time interval t  are macroscopic. It follows that the net flux of  exchanged between both media per unit surface reads

)] , , ( ) , , T , ( f ) , , ( ) , , T , ( f [ d 8 1 A B r B B A r A 0 3 total   q q v q q v dq                    . ( 4
)
At equilibrium

0 total    and ) , , T , ( f ) , , T , ( f r B r A q v q v    so that from (4) we have        A B B A .
On the other hand, from fluctuational electrodynamics [1,[START_REF] Rytov | Theory of Electrical Fluctuation and Thermal Radiation (Academy of Science of USSR[END_REF] or quantum field [18] theory any  -flux can be casted into a general form

) ( ) v , T , T , ( F d r B A 0 total max           (5)
Let us separate propagative and evanescent modes contributions.

We have

] ), ( R ), ( R [ L ) ( B A prop q q q dq /c q       and ] ), ( R ), ( R [ L ) ( B A eva q q q dq /c q      
. In this generic expression  F depends on the type of transfer that is considered and prop L and eva L are functionals of the reflectivity B , A R of exchanging media. Due to the exponential damping of evanescent waves, q is limited, at a given frequency, to a value below the cutt-off wavector

2 2 c ) c / ( / 4 q     .
According to the fundamental lemma of calculus of variations, when the monochromatic flux is extremal, the real and imaginary parts of reflectivity coefficients satisfy the so called Euler-Lagrange's

(EL) equations 0 ] R Re[ L B , A eva , prop    and 0 ] R Im[ L B , A eva , prop    (6)
Now let us consider the transfer of heat between two bodies in non-equilibrium thermal situation.

According to the fluctuational electrodynamic theory [1,[6][7] and using the azimuthal symmetry of

problem (i.e.     qdq 2 dq
) we can see that

   max , 2 i 2 B A 2 B 2 A 2 prop e R R 1 ) R 1 )( R 1 ( q 4 1 L         ,              2 2 2 B A B A 2 eva e e R R 1 ) R Im( ) R Im( q 1 L and ) T , ( ) T , ( F B A       , ] 1 - ) T k / /[exp( ≡ ) T , ( B      
being the mean energy of a Planck oscillator at equilibrium and

2 2 q ) c / (     (with 0 Im      
) the normal component of the wave vector in the intracavity space. Then, given the reflectivity of one of two interacting media, let say B R , the EL's equations ( 6) lead, after a straightforward calculation, to the solution in the non-radiative range (i.e. (7) while in the particular case of identical media

c / q   )      2 ) R ( iArg opt A e e R B
2 / 0 any for e e R 2 i opt A           . ( 8
)
Note that for both geometrical configurations 0 R opt A  is the optimal reflectivity in the radiative range.

It follows from these expressions that

1 e R R 1 ) R 1 )( R 1 ( 2 i 2 B opt A 2 B 2 opt A        and 4 1 e e R R 1 ) R Im( ) R Im( 2 2 2 B opt A B opt A           
. Therefore, by identifying (4) and ( 5) we obtain

                               c / q , e e R R 1 ) R Im( ) R Im( 4 c / q , e R R 1 ) R 1 )( R 1 ( ) , q , ( 2 2 2 B A B A 2 i 2 B A 2 B 2 A e B A     and qF 4 1 ) q , v , T , ( f 2 r A    for any q. ( 9
) For c / q   , 1 e B A   
corresponds to an energy exchange between blackbodies that is a transfer of radiative waves between two perfect emitters. In this case, by performing integration of flux over all the spectrum it is direct to see from (4), after summation over the two polarizations states, that we recover the Stefan-Boltzmann law [START_REF] Planck | The Theory of Heat Radiation[END_REF] (i.e.

) T T ( 4 B 4 A e B A      ). When c / q   , the condition 1 e B A   
corresponds to a perfect tunnelling of non-radiative photons. Also, we see from (4) that this transfer is maximal at a given frequency when the number of coupled (evanescent) modes per unit

surface dq ) q , ( 4 q ) ( N c q c / e B A 2        
becomes maximum. This precisely occurs when

1 e B A    . In this case 2 2 max 2 1 N ) ( N     
so that, by taking into account the two polarization states of nonradiative photons, the upper limit for the near-field heat transfer between two media reads

] T T [ 6 k )] T , ( ) T , ( [ d N 2 2 B 2 A 2 2 B B A 0 max max , e B A               . ( 10
)
If we assume that

T T T B A    with 1 T / T A   we can introduce a heat transfer coefficient from T h e e B A    
. From (10) after linearization, we see that

2 0 max e g 2 h    where h 3 / T k g A 2 B 2 0  
is the quantum of thermal conductance at A T . Note that the separation distance  cannot go below the scales where non local effects appear. Thus, for metals it is limited by the Thomas-Fermi screening length while for dieletrics it is the interatomic distance a which defines the lower limit [19]. For silicon, interatomic distance is roughly 0.24 nm so that the ultimate conductance at 300 K between two silicon samples is approximately

    at 300 K),
we remark that it can hardly beat the classical thermal conduction of bulk materials.

If we turn out to the non-contact friction problem at zero temperature between two reservoirs in parallel relative motion at non-relativistic velocity r v , then according to the quantum field theory

[4,17] r x max v q   , 0 L prop  and    q 2 2 q 2 B A B A x 3 eva e e R R 1 ) R Im( ) R Im( q 4 L        ,
where

) v q ( R R r x B B     and 2 y 2 x q q q  
. In that case, frictional stress optimization

           r x v
                 . ( 12
)
Now let us compare energy we can exchange between two semi-infinite media with the ultimate values predicted above. Assuming these media are identical with a dielectric permittivity given by a Lorentz-Drude model

                             (13)
It appears from this expression and from numerical calculation of heat transfer at T=300K when that it reaches about 25% of the ultimate heat transfer. Polaritons are thus excellent candidates to maximize non contact heat transfer between materials where they could be used in next generation thermophotovoltaic devices to enhance conversion of radiation into electricity. These results are in full agreement with the works of Wang et al. [20] on the magnification of transfers between two optimized dielectrics.

We have theoretically derived the fundamental limits, for energy and momentum exchanges between two parallel layered media separated by a vacuum gap at rest in nonequilibrium thermal situation and in relative motion at zero temperature respectively. The corresponding optimal reflection coefficients have been precisely determined using basic principles of calculus of variations. Our approach is general provided that non-local effects can be neglected and the sliding velocity is nonrelativistic. It could be extended to optimize exchanges between any structure shapes using the concept of generalized scattering operators rather than that of reflectivities. Quantum friction at finite temperature could also be investigated in the same way. In addition, a Landauer-type formulation of non-contact exchanges has been presented both for radiative and non-radiative modes. We have shown that while the near-field heat transfer at a given frequency is maximum when the number of coupled evanescent modes per unit surface is maximum, the momentum transfer mediated by shearing depends on the energy density of electromagnetic modes on the surfaces in interaction. These results should provide a guidance for the design of composite materials dedicated to exalted transfers in near-field technology.

Captions list

  follows from (5) that the upper limit for the frictional stress between two media in relative

  like optical phonon pulsations while  is a damping factor. These media support a surface polariton at

.

  In the neighbourhood of this curve the transfer of heat by tunnelling is very efficient. At the surface polariton frequency and far from frequency as shown in Fig.1. In this figure, note that  goes to high q values and contributes to high heat transfer values. Comparing the ratio between the heat transfer coefficient e h due to polaritons and the

Fig. 1

 1 Fig. 1 Transmission probability in p-polarization of quantum of heat between two (massive) samples

Figure 1

 1 Figure 1


Gael Moneron

Albert-Claude Boccara

Arnaud Dubois
email: arnaud.dubois@institutoptique.fr

Polarization-sensitive full-field optical coherence tomography

We present a polarization-sensitive full-field optical coherence tomography system that can produce highresolution images of the linear retardance and reflectivity properties of biological media. En face images can be delivered at a frame rate of 3.5 Hz by combination of interferometric images acquired by two CCD cameras in an interference microscope illuminated with a tungsten halogen lamp. Isotropic spatial resolution of ϳ1.0 m is achieved. The technique is demonstrated on ex vivo muscle tissues.

Optical coherence tomography (OCT) is a technique of choice for noninvasive three-dimensional imaging of biological tissues with micrometer-scale resolution [1,2]. OCT measures the amplitude of the light backscattered by the object illuminated by a broadband source. Polarization-sensitive OCT (PS-OCT) is an extension of OCT that provides additional information on the polarization properties of the tissues. PS-OCT can be performed in the time domain [3][4][5][6][7] or in the frequency domain [8,9]. The sensitivity advantage of frequency-domain OCT over time-domain OCT can be used either to shorten the acquisition time or to improve the signal-to-noise ratio. Many biological tissues exhibit birefringence related to various components such as collagen, myelin, and elastic fibers. PS-OCT can offer enhanced contrasts and specificity in structure identification. Conventional PS-OCT experimental setups require scanning the sample transversally to build up images. Fullfield OCT, also sometimes referred to as full-field optical coherence microscopy, was proposed a few years ago as an alternative method to perform ultrahighresolution OCT using a CCD camera to avoid lateral scanning and using a simple white-light source such as a tungsten halogen lamp [10][11][12]. We present in this paper a new full-field OCT technique that can also provide linear retardance measurements in real time.

The schematic of the experimental setup is represented in Fig. 1. It is based on a bulk Michelson interferometer with identical water-immersion microscope objectives in both arms. This configuration is referred to as a Linnik interference microscope [13]. A reference mirror with 2% reflectivity is realized by the interface between the polished surface of a Y 3 Al 5 O 12 (YAG) crystal and water. Ultrastable highpower spatially and temporally incoherent illumination is achieved by using a 12 V, 100 W Köhler halogen illuminator with a stabilized DC power supply. Light is linearly polarized by a polarizer (P) oriented at a 0°angle and split into the reference and sample arms by a broadband polarization-insensitive beam splitter (BS). Light in the reference arm passes through an achromatic quarter-wave plate (AQWP) oriented at a 22.5°angle to the incident polarization, is then reflected by a reference mirror, and passes through the AQWP again. Reference light is thus linearly polarized at a 45°angle to the incident polarization. Light in the sample arm passes through an identical AQWP oriented at a 45°angle, so that the sample is illuminated with circularly polarized light. Returning from the sample and passing through the AQWP, light has an arbitrary (elliptical) polarization state due to the optical anisotropy properties of the sample. After recombination by the beam splitter, the sample and reference images interfere. The interferometric signal is split into two orthogonal polarization components by a polarizing beam splitter (PBS) oriented at 0°( ʈ ) and 90°(Ќ) angles to the incident polarization. The two resulting images are projected onto two CCD arrays (Dalsa 1M15, 1024 ϫ 1024 pixels, 12 bits) by an achromatic doublet lens with a 300 mm focal length. 2 ϫ 2 binning is applied to increase the effective dynamic range of the CCD cameras. Because of the spectral response of the silicon-based CCDs and the use of a long-wave-pass filter, the effective spectrum of the illumination light has a width of ⌬ = 300 nm (FWHM) centered at = 750 nm, which fits the bandwidth acceptance of the achromatic wave plates. The reference mirror is attached to a piezoelectric translation stage (PZT) to make it oscillate at the frequency f = 3.5 Hz. This oscillation generates a sinusoidal phase modulation in the interferometer. The intensity on each pixel ͑x , y͒ of the CCDs is then time modulated at the frequency f. The CCD cameras are triggered at 4f = 14 Hz and synchronized with the PZT oscillation. The CCD cameras capture four images per modulation period. An adjustable number N of series of four images are accumulated in four buffers on the computer. By combining the four frames

E 1 Ќ , E 2 Ќ , E 3 Ќ , E 4 Ќ , and E 1 ʈ , E 2 ʈ , E 3 ʈ , E 4
ʈ , acquired by CCD 1 (Ќ) and CCD 2 ( ʈ ) respectively, according to

⌫ Ќ = ͑E 1 Ќ -E 2 Ќ ͒ 2 + ͑E 3 Ќ -E 4 Ќ ͒ 2 , ⌫ ʈ = ͑E 1 ʈ -E 2 ʈ ͒ 2 + ͑E 3 ʈ -E 4 ʈ ͒ 2 ,
we can show, using the Jones matrix formalism, that

⌫ Ќ ϰ R cos 2 ͑͒, ⌫ ʈ ϰ R sin 2 ͑͒,
where denotes the single-pass phase retardation due to the object birefringence and R is the object reflectivity. Phase-retardation images ͑x , y͒ are then obtained by calculating the ratio ⌫ ʈ / ⌫ Ќ . Polarizationindependent tissue reflectivity images R͑x , y͒ can also be obtained from the sum ⌫ Ќ + ⌫ ʈ . Data acquisition, computation, and display of the computed images are performed by home-written Visual C + + software. Images are produced at a frame rate of several hertz, depending on the accumulated image number N (3.5 Hz for N =1). Conventional PS-OCT requires a large depth of field to avoid focus adjusting while acquiring the axial dimension of the image. Very low numerical aperture (NA) objectives are therefore used, which limits the transverse resolution ͑10-20 m͒. Our fullfield OCT technique produces en face tomographic images, which eliminates the limitation to such low NA values. We use here water-immersion microscope objectives with a moderate NA value of 0.3, which yields a theoretical transverse resolution of ϳ1 m without introducing detectable polarization artifacts that could be present with a higher NA. The axial resolution of the images is determined here by the effective coherence length of the illumination source and not by the objective NA. In a medium with refractive index n = 1.33 (water), the theoretical axial resolution is ⌬z =1 m. Since biological tissues consist mainly of water, the use of identical waterimmersion microscope objectives minimizes dispersion mismatch when imaging deep inside the tissues. Nevertheless, residual dispersion mismatches and optical aberrations due to refractive index inhomogeneities within the tissues inevitably cause a slight degradation of the axial and transverse resolutions compared with the theoretical value of 1 m.

The instrument accuracy for linear retardance measurements was tested by measuring the phase retardation introduced by a Babinet compensator placed into the sample arm between the the quarter wave-plate (AQWP) and the microscope objective (MO). A glass plate with appropriate thickness was placed in the other arm of the interferometer to compensate for the dispersion. A phase retardation varying between 0°and 225°was generated by the calibrated Babinet compensator. Figure 2 compares the measured phase retardation with the true phase retardation introduced by the compensator. The RMS error was 3°(7% average error over the whole scale). The error may result from the calculations, which assume a monochromatic illumination. A rotation of the Babinet compensator in the plan orthogonal to the optical axis had no effect on the measured value of the phase retardation, which was a confirmation of the insensitivity of the setup to the object orientation. When imaging a birefringent object with a microscope objective, the angle between the optical rays and the birefringence optic axis is not as well defined. It has recently been shown that a substantial discrepancy between measured retardance and true retardance can be induced by variations in this angle [14]. We also point out that the phase retardation at a depth z under the object surface does not represent a strictly depth-resolved measurement [15], but the accumulated linear retardance, i.e., the linear retardance integrated along the depth z and over the optical ray angle distribution. To illustrate the system performances for biomedical applications, we imaged fresh muscle tissues extracted from a shrimp tail. The tissues were placed in a tank filled with phosphate-buffered saline (PBS, pH 7) for preservation. En face oriented images are shown in Fig. 3 using 256 gray levels. Linear retardance varies from 0°(black) to 90°(white). Intensity, on a logarithmic scale, varies from -90 dB (black) to -45 dB (white). Couples of images were acquired at two different depths within the tissues, separated by 70 m. Intensity images and phase images exhibit different and uncorrelated contrasts. One can notice a dark area on Fig. 3(c) that becomes bright on Fig 3(d), indicating a variation of ϳ70°in phase retardation, whereas the gray levels in some other areas remain almost unchanged. This clearly reveals optical anisotropy inhomogeneities in the tissues.

In summary, we have developed a polarizationsensitive full-field OCT system that produces images of biological media with both intensity-based and linear-retardance-based contrasts. Compared with conventional full-field OCT, polarization-sensitive full-field OCT provides additional information related to the object polarization properties. With a moderate increase in experimental setup complexity and price, polarization full-field OCT offers larger potentialities for biomedical applications.

Fig. 1 .

 1 Fig. 1. Experimental arrangement. P, polarizer; BS, beam splitter (broadband, nonpolarizing); PBS, polarizing beam splitter (broadband); MO, microscope objective (water immersion, ϫ10, 0.3 NA); L, achromatic doublet lens (300 mm focal length). The reference mirror (2% reflectivity) oscillates at 3.5 Hz. The synchronized CCD cameras are triggered at 14 Hz to acquire four images per modulation period.

Fig. 2 .

 2 Fig. 2. Measured versus actual phase retardation introduced by a calibrated Babinet compensator.

July 15, 2007 / Vol. 32, No. 14 / OPTICS LETTERS