Chouaib Afri

Madiha Nadri Wolf

Pascal Dufour

Observer Design for a Ternary Distillation Column with Side Stream

This paper investigates the problem of observer design to provide temperature and molar fraction spatial profiles in a ternary distillation column with side stream. Indeed, both temperature and molar fraction spatial profiles are required to be used in advanced controller design. Such an estimation problem is usually solved using an implicit observer which requires an optimization algorithm when a differential algebraic system is used. In this paper, we consider a realistic nonlinear differential algebraic equation based model validated experimentally. We propose an approach to design an explicit observer that reconstructs all the state variables with an asymptotic convergence. The observer is build as a coupling of two nonlinear observers in full interconnection. Finally, the performance of the proposed observer is illustrated using a model of a ternary distillation column with side stream.

I. INTRODUCTION

The industry is still looking for better strategies for improving purity control in distillation columns, while optimizing installation costs and energy. The constraints imposed on such system are: the high cost of the measures and the nonlinear coupling between physical quantities. The most common controllers in the industrial world are PID, which tuning may be based on a simple model to describe the process. Their decentralized control loop feature results in high energy consumption and poor control performances. Hence, advanced controllers are needed, using a model based state estimation. The model used here is based on the total and partial material balances and the vapor liquid equilibrium (VLE) [START_REF] Gmehling | Chemical thermodynamics for process simulation[END_REF]. In this work, we do not consider the classical assumption of ideal mixture: the components to be separated can have a completely different chemical structures (such as oil and water) resulting in forces of attraction or repulsion between adjacent molecules of different types. This phenomenon is reflected by the activity coefficients which are functions of the component concentrations of the liquid mixture. A few works have considered the nonideal case [START_REF] Baratti | A composition estimator for multicomponent distillation columns-development and experimental test on ternary mixtures[END_REF]. Unlike to the classical models based on LEWIS assumptions [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF], which can be represented as an ordinary differential equations (ODE) system, the model considered here is closer to the real behaviour and is a set of nonlinear differential algebraic equations (DAE). Using algebraic methods, observer design for linear DAE systems has widely been investigated in the literature (see for instance [START_REF] Muller | On the observer design for descriptor systems[END_REF], [5], [START_REF] Darouach | Reduced order observer design for descriptor systems with unknown inputs[END_REF], [START_REF] Boukhobza | Observability of structured linear systems in descriptor form: A graph-theoretic approach[END_REF] and the references herein). [START_REF] Nikoukhah | Kalman filtering and Riccati equations for descriptor systems[END_REF] extends the Kalman filter for discrete time systems (the proposed Kalman filter is also an implicit system). The observer synthesis for nonlinear DAE systems is still an open problem. The most C. Afri, M. Nadri Wolf and P. Dufour are with the Université de Lyon, F-69622, Lyon, France -Université Lyon 1, Villeurbanne, France -CNRS, UMR 5007, LAGEP, France. {nadri}@lagep.univ-lyon1.fr common approaches consist in characterizing this DAE class which can be transformed by a diffeomorphism to an ODE system and for which an observer can be designed [START_REF] Krener | Nonlinear observers with linear error dynamics[END_REF], [START_REF] Hammouri | Bilinearization up to output injection[END_REF], [START_REF] Gauthier | Observability and observers for nonlinear systems[END_REF]. The observer for the original system is an implicit system obtained by using the pullback map.

Beside these geometric approaches, other approaches gave implicit observers whose practical implementation usually requires the use of optimization techniques to solve at each time t the algebraic constraints associated to the state manifold together with an ODE numerical method (see for instance [START_REF] Tjoa | Simultaneous solution and optimization strategies for parameter estimation of differentialalgebraic equation systems[END_REF], [START_REF] Becerra | Applying the extended Kalman filter to systems described by nonlinear differential algebraic equations[END_REF], [START_REF] Nikoukhah | A new methodology for observer design and implementation[END_REF] and [START_REF] Aslund | An observer for non-linear differentialalgebraic systems[END_REF]). Such methods have been shown to lack robustness when solving large scale complex systems or systems with fast dynamics. Another interesting approach is presented in [START_REF] Zimmer | On observing nonlinear descriptor systems[END_REF]. In their work, the authors used the algebraic constraints as additional output and gave an explicit observer of higher dimension. This work aims to design an explicit state observer based on a realistic experimentally validated nonlinear DAE model. The outline of the paper is as follows. In section 2, we summarize the model development for a ternary distillation column with side stream. In section 3, we synthesize an explicit observer for a class of nonlinear DAE systems which contains most models of distillation columns. Finally, the usefulness of this approach is illustrated using only four molar fraction measurements to estimate both temperature and molar fractions spatial profiles.

II. STATE MODEL OF A TERNARY DISTILLATION COLUMN WITH SIDE STREAM

The considered process is a ternary distillation column with side stream (figure 1) separating a liquid ternary mixture of BTX (Benzene, Toluene and o-Xylene, considered as the component 1, 2 and 3 respectively) [START_REF] Gmehling | Vapor-liquide equilibrium data colection, Aromatic Hydrocarbons[END_REF]. The column consists of (n-2) theoretical trays, a boiler (tray n)and a condenser (tray 1). Using the material balance, we get the DAE (1), where, j = 1, 2, 3 is the component index, i = 1, ..., n is the tray index, N i is the molar retention in the tray i, x j,i the liquid phase molar fraction of component j in the tray i, y j,i is the vapour phase molar fraction of component j in the tray i, T i is the stage temperature in the tray i, Z j,l is the feed flow molar fraction of component j in the tray l, F is the molar feed rate in the tray l, S is the molar flow rate extracted from tray s, L is the molar flow rate of reflux and V f is the molar vapour flow rate. Fractions of the third component is calculated from the equation of matter conservation (2).

                                               total condenser, i = 1 : N 1 dx j,1 dt = V f (y j,2 -x j,1 )
upper rectifying section, i = 2, ..., s : N i dx j,i dt = L(x j,i-1x j,i ) +V f (y j,i+1y j,i ) lower rectifying section, i = s + 1, ..., l -1 : N i dx j,i dt = (L -S)(x j,i-1x j,i ) +V f (y j,i+1y j,i ) feed tray, i = l : N l dx j,l dt = (L -S)(x j,l-1x j,l ) +F(Z j,lx j,l ) +V f (y j,l+1y j,l ) stripping section, i = l + 1, ..., n -1 :

N i dx j,i dt = (F + L -S)(x j,i-1 -x j,i ) +V f (y j,i+1 -y j, i) boiler, i = n : N n dx j,n dt = (F + L -S)(x j,n-1 -x j,n ) +V f (x j,n -y j,n ) thermodynamic constraints, i = 1, .., n : (∑ c j=1 y j,i ) -1 = 0 (1) 3 ∑ j=1 x j,i -1 = 0, ∀i ∈ [1, n]. (2)
From now, based on (2), we set that x 3,i is a function of x 1,i and x 2,i . From the VLE and the Antoine equation [START_REF] Gmehling | Vapor-liquide equilibrium data colection, Aromatic Hydrocarbons[END_REF], [START_REF] Gmehling | Chemical thermodynamics for process simulation[END_REF] we get:

                     rectifying section: 1 ≤ i < l y j,i = γ j,i (x 1,i ,x 2,i ,T i ) P T e l i ((x j,i -(1 -e l i )(x j,i-1 ))10 a j - b j T i +c j ) stripping section: l ≤ i ≤ n -1 y j,i = (1 -e v i ) y j,i+1 + e v i γ j,i (x 1,i ,x 2,i ,T i ) P T (x j,i 10 a j - b j T i +c j )
boiler: y j,n = 1 P T (x j,n 10

a j - b j Tn+c j γ j,n (x 1,i , x 2,i , T n )) (3)
where, γ j,i is the activity coefficient that represents the non ideality of mixture, e l i , e v i are respectively the rectifying and stripping efficiencies [START_REF] Gmehling | Chemical thermodynamics for process simulation[END_REF], and P T the internal pressure. The γ j,i 's and Antoine parameters (a j , b j and c j ) are obtained from binary mixtures experimental data of the three components BTX [START_REF] Gmehling | Vapor-liquide equilibrium data colection, Aromatic Hydrocarbons[END_REF]. This choice is the best in terms of accuracy compared to the ternary experimental data [START_REF] Baratti | A composition estimator for multicomponent distillation columns-development and experimental test on ternary mixtures[END_REF]. The difference between the classical models [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF] and this one is that the mixture here is not ideal. Consequently, γ j,i is a function of liquid molar fractions x j,i and temperature T i . We divide the process to a rectifying section and a stripping section, with exponents 1 and 2 respectively, and introduce the following notations :

                           x 1 i = x 1,i ; 1 ≤ i ≤ l -1 x 2 i = x 1,n-i+1 ; 1 ≤ i ≤ n -l -1 x1 i = x 2,i ; 1 ≤ i ≤ l -1 x2 i = x 2,n-i+1 ; 1 ≤ i ≤ n -l -1 z 1 i = T i ; 1 ≤ i ≤ l -1 z 2 i = T n-i+1 ; 1 ≤ i ≤ n -l -1 u T = (L,V f , F, S, Z 1,l , Z 2,l ) T y 1 = (x 1 1 , x 2 1 ) T = (x 1,1 , x 1,n ) T y 2 = ( x1 1 , x2 1 ) T = (x 2,1 , x 2,n ) T (4)
where y 1 contains the top and bottom measures of light component (Benzene), y 2 contains the top and bottom measures of intermediate component (Toluene). F, S, Z 1,l and Z 2,l are the unknown input disturbances, V f , L are manipulated inputs of the input vector u. x and x are the state variables. Using (4), the system (1)-( 3) belongs to the following fully interconnected systems [START_REF] Besanc ¸on | On observer design for interconnected systems[END_REF]:

(Σ)    ẋ(t) = f (x(t), x(t), z(t), u(t)) ϕ(x(t), x(t), z(t)) = 0 y 1 (t) = h(x) = [x 1 1 (t), x 2 1 (t)] T (5)
(Σ)    ẋ(t) = f (x(t), x(t), z(t), u(t)) ϕ(x(t), x(t), z(t)) = 0 y 2 (t) = h( x) = [ x1 1 (t), x2 1 (t)] T (6)
where,

x = (x 1 1 , .., x 1 l-1 , x 2 1 , .., x 2 n-l+1 ) T ; x = ( x1 1 , .., x1 l-1 , x2 1 , .., x2 n-l+1 ) T f = f 1 ( x, x, z, u) f 2 (x, x, z, u) ; f = f 1 ( x, x, z, u) f 2 (x, x, z, u) ; ϕ = ϕ 1 (x 1 , x1 , z 1 ) ϕ 2 (x 2 , x2 , z 2 )
f , f and ϕ are C 1 functions with respect to their arguments and satisfy the following triangular structure:

C1) f i =         f i 1 (x i 1 , x i 2 , xi , z i 2 , u) : f i k (x i k-1 , x i k , x i k+1 , xi , z i k , z i k+1 , u) : f i n i -1 (x i n i -2 , x i n i -1 , x i n i , xi , z i n i-1 , z i n i , u) f i n i (x, x, z, u)         C2) f i =         f i 1 ( xi 1 , xi 2 , x i , z i 2 , u) : f i k ( xi k-1 , xi k , xi k+1 , x i , z i k , z i k+1 , u) : f i n i -1 ( xi n i -2 , xi n i -1 , xi n i , xi , z i n i-1 , z i n i , u) f i n i ( x, x, z, u)         C3) ϕ i =       ϕ i 1 (x i 1 , xi 1 , z i 1 ) : ϕ i k (x i k-1 , x i k , x i k+1 , xi k-1 , xi k , xi k+1 , z i k ) : ϕ i 1 (x i n i , xi n i , z i n i )      
where n 1 = l -1, and n 2 = nl + 1.

III. EXPLICIT OBSERVER DESIGN FOR THE FULLY INTERCONNECTED IMPLICIT SYSTEMS

The system under consideration (Σ,Σ) belongs to the class of implicit systems written in the following compact format:

       ẋ(t) = f (x(t), x(t), z(t), u(t)) x ∈ R n , u ∈ U ẋ(t) = f (x(t), x(t), z(t), u(t)) x ∈ R n ϕ(x(t), x(t), z(t)) = 0 z ∈ R d y(t) = [h(x(t)) h( x(t))] T y ∈ R p , (7)
where

U is a bounded open subset of R m , (x(t), x(t), z(t)) ∈ R n × R n × R d
is the state vecteur, f , f , h and h are assumed to be sufficiently smooth with respect to their arguments. Let M denotes the constraints set,

M = (x, x, z) ∈ R n × R n × R d , s.t. ϕ(x, x, z) = 0 (8)
In [START_REF] Hammouri | An observer design for a class of implicit systems[END_REF] the authors consider a class of implicit systems of the form (5) and gave an high gain explicit observer. Based on this result, we will show that the following system ( 9) is an asymptotic observer for system [START_REF] Boukhobza | Observability of structured linear systems in descriptor form: A graph-theoretic approach[END_REF].

                     ẋ(t) = f ( x(t), x(t), ẑ(t), g(t), h( x(t)), u(t)) ẋ(t) = f ( x(t), x(t), ẑ(t), ḡ(t), h( x(t)), u(t)) ż(t) = -∂ ϕ ∂ z | ( x,ẑ, x) -1 × ∂ ϕ ∂ x | ( x,ẑ) , ∂ ϕ ∂ x | ( x,ẑ) ẋ , ẋ T + Λϕ( x, x, ẑ) ġ(t) = G(g(t), y 1 (t), u(t)) ġ(t) = G( ḡ(t), y 2 (t), u(t)), (9)
where, Λ is a d ×d symmetric positive definite matrix (SPD), g and ḡ in R n×p are the observer gain functions defined in the next sections. The following assumptions are considered.

Assumption A 1): There exists ε 0 > 0 such that ∂ ϕ ∂ z has a full rank on the tubular neighborhood Ω 0 = {(x, x, z) ∈ R n+n+d ; ϕ(x, x, z) < ε 0 }, and there exists a function σ of class K ∞ such that:

i) The following map:

Φ : Ω 0 -→ R n+n × J (x, x, z) -→ (x, x, ϕ(x, x, z)) ( 10
) is a diffeomorphism, where J is an open subset of R d . ii) ∀(x, x, z) ∈ Ω 0 , ϕ(x, x, z) ≥ σ ( z -z * ) z * , is the solution of ϕ(x, x, z) = 0.
The next hypothesis guaranties the existence of an implicit observer: Assumption A 2):

1) There exists U ⊂ L ∞ (R + ,U), such that for every u ∈ U ; for every (g(0)), ( ḡ(0)) ∈ O, the trajectory g(.) and ḡ(.) are bounded. 2) There exist positive definite functions V (t, e), V (t, ē), W (e) and W ( ē) defined on R + ×R n such that for every u ∈ U , we have:

i) ∀(x, x, z) ∈ M ; ∀( x, x, ẑ) ∈ M ; ∀t ≥ 0; we have ∂V ∂t (t, e) + ∂V ∂ e (t, e) f (e + x, x, z * , g, h( x), u) -f (x, x, z, u) ≤ -W (e) (11)
∂V ∂t (t, ē) + ∂V ∂ ē (t, ē) f ( x, ē + x, z * , ḡ, h( x), u) -f ( x, x, z, u)) ≤ -W ( ē) (12
) where e = xx and ē = xx. ii) There exist constants a ≥ 0, α > 0 and ᾱ > 0;

such that ∀e, ē; ∀t ≥ 0

∂V ∂ e (t, e) + ∂V ∂ ē (t, ē) ≤ αW (e) + ᾱW ( ē) e + ē ≥ a. iii) ∃α 1 > 0, α 2 > 0; ∀(x, x) ∈ R n × R n ; ∀(e, ē) ∈ R n × R n ; ∀t ≥ 0, ∂V ∂ e ( f (x, x + ē, z, u) -f (x, x, z, u)) ≤ α 1 W (e) W ( ē), ∂V ∂ ē ( f (x + e, x, z, u) -f (x, x, z, u)) ≤ α 2 W (e) W ( ē), 0 < α 1 + α 2 ≤ 2.
Finally, we consider the Lipschitz condition: Assumption A 3): There exists a class K ∞ functions λ 1 (.) and λ 2 (.) such that, ∀( x, x, z, z * ( x, x)) ∈ Ω 0 ; ∀g, we have:

f ( x, x, z, g, h( x), u) -f ( x, x, z * , g, h( x), u) ≤ λ 1 ( z -z * ) f ( x, x, z, ḡ, h( x), u) -f ( x, x, z * , ḡ, h( x), u) ≤ λ 2 ( z -z * )
Now consider Λ defined as above, and let ε Λ > 0 be such that

Ω Λ = {(x, x, z) ∈ R n × R n × R d ; ϕ(x, x, z) T Λϕ(x, x, z) < ε Λ } ⊂ Ω 0 . Noticing that Ω Λ = Ω 0 if Λ = I and ε Λ = ε 0 .
Using these notations, we can state our main result as follow.

Theorem 1: Under assumptions A 1), A 2) and A 3) , system (9) forms an asymptotic explicit observer for system [START_REF] Boukhobza | Observability of structured linear systems in descriptor form: A graph-theoretic approach[END_REF].

Proof: Let(x(t), x(t), z(t)) and ( x(t), x(t), ẑ(t)) be the trajectories of systems ( 7) and ( 9) respectively. From (9) we can get:

d ϕ( x(t), x(t), ẑ(t)) 2 dt = -2(ϕ( x, x, ẑ)) T Λ(ϕ( x, x, ẑ)) (13)
From assumption A 1), Φ is a diffeomorphism from Ω 0 into R n+n × J, so there exists a unique solution z * (t), such that ϕ( x(t), x(t), z * (t)) = 0. Let e(t) = x(t)x(t), ē(t) = x(t)x(t) be the estimation errors of the observer. So from [START_REF] Boukhobza | Observability of structured linear systems in descriptor form: A graph-theoretic approach[END_REF] and ( 9) we can get:

ė(t) + ė(t) = f ( x, x, z * , g, h( x), u) -f (x, x, z, u) + f ( x, x, z * , ḡ, h( x), u) -f ( x, x, z, u) + f (x, x, z, u) -f (x, x, z, u) + f ( x, x, z, u) -f (x, x, z, u) + f ( x, x, ẑ, g, h( x), u) -f ( x, x, z * , g, h( x), u) + f ( x, x, ẑ, ḡ, h( x), u) -f ( x, x, z * , ḡ, h( x), u). (14)
Considering V (t, e), V (t, ē), W (t, e) and W (t, ē) given in assumption A 2) -2) and using A 2) -2)i), we get:

dV (t, e, ē) dt = dV (t, e) dt + dV (t, ē) dt ≤ -W (e) -W ( ē) + ∂V (t, e) ∂ e ( f (x, x, z, u) -f (x, x, z, u)) + ∂V (t, ē) ∂ ē ( f ( x, x, z, u) -f (x, x, z, u)) + ∂V (t, e) ∂ e f ( x, x, ẑ, g, h( x), u) -f ( x, x, z * , g, h( x), u) + ∂V (t, ē) ∂ ē f ( x, x, ẑ, ḡ, h( x), u) -f ( x, x, z * , ḡ, h( x), u) .
(15) And from assumption A 2) -2)iii), we get:

dV (t, e, ē) dt ≤ -W (e) -W ( ē) +α 1 W (e) W (e) + α 2 W (e) W (e) + ∂V (t, e) ∂ e f ( x, x, ẑ, g, h( x), u) -f ( x, x, z * , g, h( x), u) + ∂V (t, ē) ∂ ē f ( x, x, ẑ, ḡ, h( x), u) -f ( x, x, z * , ḡ, h( x), u) .
(16) From assumption A 3), we deduce that:

dV (t, e, ē) dt ≤ - 1 2 (1 - α 1 + α 2 2 )( W (e) + W ( ē)) 2 +(( ∂V (t, e) ∂ e λ 1 ) + ( ∂V (t, ē) ∂ ē λ 2 )) ẑ -z * . (17)
And finally, from assumption A 1)ii), we get:

dV (t, e, ē) dt ≤ - 1 2 (1 - α 1 + α 2 2 )( W (e) + W ( ē)) 2 +(( ∂V 1 (t, e) ∂ e λ 1 ) + ( ∂V 2 (t, ē) ∂ ē λ 2 ))σ -1 ϕ( x, x, ẑ) . ( 18
) From ( 13), we have lim t→∞ ϕ( x(t), x(t), ẑ) = 0 and using assumption A 2)ii) in [START_REF] Hammouri | An observer design for a class of implicit systems[END_REF] we deduce that lim

t→∞ |e(t)| + | ē(t)| = 0. From assumption A 1) -ii), ( ẑ -z ) ≤ σ -1 ϕ( x, x, ẑ)
where z is the solution of: ϕ(x, x, z) = 0, using [START_REF] Becerra | Applying the extended Kalman filter to systems described by nonlinear differential algebraic equations[END_REF] we obtain lim t→∞ ẑ(t)z(t) = 0.

IV. EXPLICIT OBSERVER FOR IMPLICIT SYSTEMS WITH TRIANGULAR STRUCTURE

The first step in the design of the observer ( 9) is the synthesis of the gain functions g and ḡ. Many normal forms which characterized some subclasses of multi-output uniformly observable systems together with their associated observer syntheses have been proposed in the literature, see for instance: [START_REF] Gauthier | Observability and observers for nonlinear systems[END_REF], [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF], and [START_REF] Hammouri | Nonlinear observers for locally uniformly observable systems[END_REF]. Now, let us come back to the structures ( 5) and [START_REF] Darouach | Reduced order observer design for descriptor systems with unknown inputs[END_REF]. Referring to (C1, C2), we remark that the functions f and f have triangular structures with respect to (h(x), x) and ( h( x), x). In this case, a constant gain g = K θ (i.e. G(g,u)=0) respectively ḡ = K θ , can be designed as stated in [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF]. As can be seen, system (5) and system [START_REF] Darouach | Reduced order observer design for descriptor systems with unknown inputs[END_REF] have the same structure. Therefore, the same observer structure will be used. Now, by considering the triangular structure of (5) given by C1), C2) and C3) and assuming that the map (x, x, z) → (x, x, ϕ(x, x, z)) satisfies the assumption A 1), we deduce that for i = 1, 2 and 1 ≤ k ≤ n i ; there exist smooth functions γ i k (their analytic expressions are generally unknown) such that the solution z of ϕ(x, x, z) = 0 is given by:

   z i 1 = γ i 1 (x i 1 , xi 1 , x i 2 , xi 2 ) z i k = γ i k (x i k-1 , xi k-1 , x i k , xi k , x i k+1 , xi k+1 ) 1 < k < n i z i n i = γ i n i (x, x).
In the sequel, we use the same notation as above. We set

F i k (x, x, u) = f i k (x, x, γ i k (x, x), u), and
F i = [F i 1 , . . . , F i n i ] T , F = [F 1 T F 2 T ] T , then system (5) can be rewritten: ẋ(t) = F(x(t), x(t), u(t)) y 1 (t) = Cx(t). (19)
Moreover F(x, x, u) is a C 1 function of triangular form and satisfy the following hypothesis.

H1) For every (x, x, u), ρ 1 ≤ ∂ F i j ∂ x i j+1 (x, x, u) ≤ ρ 2 , for some constants ρ 1 > 0, ρ 2 > 0. H2) F is a global Lipschitz function w.r.t. x, ∃c > 0; ∀u ∈ U; ∀x, x , F(x, x, u) -F(x , x, u) ≤ c x -x .
In [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF] the authors gave a constant gain observer for systems of the form ( 19) with x as a known variable. Based in this result, a constant high gain observer can be constructed for system [START_REF] Hammouri | Nonlinear observers for locally uniformly observable systems[END_REF] as follows:

ẋ1 (t) = F 1 ( x(t), x, u(t)) -r 1 ∆ θ δ 1 S -1 1 C T 1 (C 1 x1 -x 1 1 ) ẋ2 (t) = F 2 ( x(t), x, u(t)) -r 2 ∆ θ δ 2 S -1 2 C T 2 (C 2 x2 -x 2 1 ), (20)
where r k ≥ 1 (k = 1, 2) is a positive constant tuning parameter of the observer,

C k = [1, . . . , 0] is a vector of dimension n k ; δ 1 > 0, δ 2 > 0 satisfy: 2l -3 2n -2l + 3 δ 1 < δ 2 < 2l -1 2n -2l + 3 δ 1 , θ > 0 and ∆ θ δ i is a (n k × n k ) diagonal matrix: ∆ θ δ k = diag(θ δ k , θ 2δ k , . . . , θ n i δ i ),
and finally S k is an n k × n k constant SPD matrix of the form:

S k =          s 11 s 12 0 0 s 12 s 22 . . . . . . 0 . . . . . . 0 . . . . . . s k-1k 0 . . . 0 s k-1k s kk .          (21)
To assure the exponential convergence of the above observer, the matrix S k is chosen to satisfy the following condition (detail of the proof see [START_REF] Hammouri | High gain observer based on a triangular structure[END_REF]):

∀t ≥ 0, A T k (t)S k + S k A k (t) -ρC T k C k ≤ -νI k , (22)
where A k (i, j) = a k i (t) for j = i+1 and A k (i, j) = 0 for j = i+ 1 such that ρ 1 ≤ a i (t) ≤ ρ 2 , and ρ 1 , ρ 2 are positive constants. Now, let Ω 0 be a tubular neighborhood of M given by the assumption A 1). In order to show that system (23) forms an observer for system (5,[START_REF] Darouach | Reduced order observer design for descriptor systems with unknown inputs[END_REF], the following assumption is used: H3): f and f is a global Lipschitz functions on Ω 0 , namely:

∃(c , c ) > 0; ∀u ∈ U; ∀(x, x, z), (x , x , z ) ∈ Ω 0 , f (x, x, z, u) -f (x , x , z , u) ≤ c (x, x, z) -(x , x , z ) . f (x, x, z, u) -f (x , x , z , u) ≤ c (x, x, z) -(x , x , z
) . Now, we can state the main result of this section:

Proposition 1: Under assumption A 2 ) and hypotheses H1), H3), the following system:

               ẋ = f ( x, x, ẑ, u) + K θ ( ŷ1 -y 1 ) ẋ = f ( x, x, ẑ, u) + K θ ( ŷ2 -y 2 ) ż(t) = -∂ ϕ ∂ z | ( x,ẑ, x) -1 × ∂ ϕ ∂ x | ( x,ẑ) , ∂ ϕ ∂ x | ( x,ẑ) ẋ , ẋ T + Λϕ( x, x, ẑ) , (23) where
K θ = diag(-r 1 ∆ θ δ 1 S -1 1 C T 1 , -r 2 ∆ θ δ 2 S -1 2 C T 2 ) K θ = diag(-r 1 ∆θ δ 1 S-1 1 CT 1 , -r2 ∆θ δ 2 S-1 2 CT 2 ),
forms an explicit observer for system [START_REF] Boukhobza | Observability of structured linear systems in descriptor form: A graph-theoretic approach[END_REF].

The proof of this proposition is accomplished by applying Theorem 1 of the above section, where f ( x, x, ẑ, h( x), g, u) = f ( x, x, ẑ, u) + K θ ( ŷ1y 1 ) and f ( x, x, ẑ, h( x), ḡ, u) = f ( x, x, ẑ, u) + K θ ( ŷ2y 2 ).

To do, it suffices to show that assumptions A 2 ) and A 3 ) hold. Setting V (t, e) = e T ∆ θ S∆ θ e and V (t, ē) = ēT ∆ θ S∆ θ ē, so V and V are positive definite quadratic functions and do not depends on t. Thus condition of assumption A 2) -2)ii) is satisfied. Let us now check the condition A 1) -1)i): By construction, we know that if (x, x, z) ∈ M , ( x, x, ẑ) ∈ M , then:

f (x, x, z, u) = F(x, x, u), f (x, x, z, u) = F(x, x, u) and f ( x, x, ẑ, g, h( x), u) = F( x, x, u) + K θ (h( x) -y 1 ), f (x, x, ẑ, ḡ, h( x), u) = F(x, x, u) + K θ ( h( x) -y 2 ).
Now since H1), H2) are satisfied for F(x, x, u), F(x, x, u), inequalities [START_REF] Gauthier | Observability and observers for nonlinear systems[END_REF] and ( 12) holds. Consequently, A 2) -2)i) is satisfied. Assumption A 3 )) is a direct consequence of hypothesis H3) and the fact that the gain of the observer is constant and the output is a linear map. In the end it remains to verify the assumption A 2) -2)iii) which is the full interconnection condition [START_REF] Besanc ¸on | On observer design for interconnected systems[END_REF], it is also a direct consequence of hypothesis H3).

V. SIMULATION RESULTS

A. Experimental validation of the model 1

In this section we will analyse the validity of the controloriented DAE model (1) developed above. It is solved using the Runge-Kutta method combined with the Newton-Raphson optimization algorithm. As first step of the model validation is the study of the model behaviour when the column is assumed to be stabilized around its steady state. To do so, we consider an experimental configuration of the steady state shown in table 1. The simulation results of this test are given in Figure 2 which shows the evolution over time of the liquid molar fractions of benzene in the condenser, the boiler and the trays 3 and 9. We note the existence of a negligible transient compared to the initial values, which prove the good accuracy of the model. 1 Authors acknowledge the company Processium for the experimental data.

Tray index: total n = 18, side stream and feed trays s = 6, l = 9

Liquid molar retention in condenser

N 1 = 20 mol
Liquid molar retention in boiler

N 18 = 20 mol
Liquid molar retention in each tray

N i = 8 mol
Liquid and vapour Murphree efficiency

e l = 1 , e v = 1
Internal pressure A second step of the model validation is done using experimental data. Figure 3 shows the simulated and experimental steady state profiles of molar fractions of the ternary mixture BTX according to the number of trays. We can note a good agreement.

B. Observer simulation

To simulate the response of the observer and compare it with the response of the implicit system, we kept the same configuration of the operating point and initial conditions given in the previous part. For the observer we initialized the estimated states to arbitrary values different from those of model. The values of the observer parameters r and θ are set to 1 and 2.8 respectively. The output measurements y 1 (t) and y 2 (t) feature a zero mean Gaussian noise of 5%. The performances of the observer are illustrated in the figures 4 and 5: in spite of the errors in the initial condition, noise measurements and the disturbances, the estimated state converge to their dynamic target values, for both the molar fraction and the temperature profiles.

VI. CONCLUSION

In this work we first developed and validated experimentally a mathematical model (a set of differential algebraic equations) that describes the dynamic of the temperature and component molar fraction profiles of a ternary mixture in a distillation column. The obtained model is an implicit systems fully interconnected based on the triangular structure of the model, an explicit high gain observer is given using bottom and top molar fractions measurement. The observer implementation is simple and requires small computational effort. Furthermore, simulation results demonstrate the good performances achieved by this observer in coping with model nonlinearities and output noise.

Fig. 1 .

 1 Fig. 1. Scheme of a distillation column with side stream.

Fig. 2 .

 2 Fig. 2. Response of the dynamic model starting from the experimental steady state, in 4 locations: liquid molar fractions of benzene.

Fig. 3 .

 3 Fig. 3. Final molar fraction profiles: simulation and experimental.

Fig. 4 .

 4 Fig. 4. Four true and estimated benzene molar fractions, in the presence of 5% noise and initial errors in the estimations, θ = 2.8.

TABLE I OPERATING

 I

POINT AND INITIAL STEADY STATE.


J Peralez
email: johan.peralez@ifpen.fr

P Tona

A Sciarretta
email: antonio.sciarretta@ifpen.fr

P Dufour
email: dufour@lagep.univ-lyon1.fr

M Nadri
email: nadri@lagep.univ-lyon1.fr

Optimal Control of a Vehicular Organic Rankine Cycle via Dynamic Programming with Adaptive Discretization Grid

Keywords: Dynamic programming, backward reachability, waste heat recovery, Rankine cycle, vehicle energy management

This paper addresses the problem of maximizing the power produced by an organic rankine cycle (ORC) waste heat recovery (WHR) system along a given vehicle mission. A simplified model of the system allows the formulation of an optimal control problem that can be solved via dynamic programming (DP). To increase the smoothness and the accuracy of the solution obtained using the state-of-the-art implementation of DP known as level-set DP, an improved version is developed, making use of adaptive grids for discretization. The adaptivegrid algorithm is first validated on a simple problem with an analytical solution, then applied to the more complex ORC case study. It is shown that the required accuracy can be achieved with a lower computational effort than the original algorithm. The analysis of the resulting optimal trajectory for the ORC control problem provides useful insight for both control design and system design.

INTRODUCTION

Over the past few decades, several problems of transient control of energy production on board vehicles in the presence of constraints have been addressed and solved via optimal control methods and in particular via dynamic programming [START_REF] Bellman | Dynamic programming[END_REF]). A typical example is the energy management (EM) problem in hybrid electric vehicles (HEVs), which consists in minimizing a cost function (most often, fuel consumption) during vehicle operation while respecting driver's demand and actuator limitations, following a prescribed battery state-of-charge trajectory [START_REF] Guzzella | Vehicle propulsion systems[END_REF]; [START_REF] Serrao | A comparative analysis of energy management strategies for hybrid electric vehicles[END_REF]). Although it is not directly implementable in real-time (as it needs perfect prediction of future driving cycle), dynamic programming proves very useful in this context, as an analysis tool to understand (optimal) HEV EM behaviour, as a benchmark to assess other EM strategies or as a method to determine the potential of a given architecture.

With new technologies continuously being introduced to further reduce vehicle fuel consumption, new control problems arise that must be efficiently solved to make the most of the innovation potential. This is for instance the case of engine waste heat recovery (WHR) systems based on the Rankine thermodynamic cycle. Applications of this technology exists, at different levels of maturity, for cars, trucks, trains and ships (see [START_REF] Sprouse | Review of organic rankine cycles for internal combustion engine exhaust waste heat recovery[END_REF] for a recent overview). For some applications, such as WHR via organic rankine cycle (ORC) from Diesel engines of long haul trucks, very interesting fuel consumption reductions have been reported (up to 6%, according to [START_REF] Stanton | Systematic development of highly efficient and clean engines to meet future commercial vehicle greenhouse gas regulations[END_REF]). As shown in [START_REF] Endo | Study on maximizing exergy in automotive engines[END_REF], and more recently in [START_REF] Peralez | Towards model-based control of a steam rankine process for engine waste heat recovery[END_REF], [START_REF] Peralez | Improving the control performance of an organic rankine cycle system for waste heat recovery from a heavy-duty diesel engine using a model-based approach[END_REF] and [START_REF] Xie | Dynamic behavior of rankine cycle system for waste heat recovery of heavy duty diesel engines under driving cycle[END_REF], control of Rankine-based WHR systems for vehicular applications is far from being trivial, especially when not addressed at initial stages of system design. On the other hand, an early assessment of the best dynamic performance achievable on a given architecture, and of the way to obtain it, can significantly facilitate control design, while providing valuable insight into the system. This paper addresses the problem of maximizing the power produced by an ORC WHR system along a given vehicle mission. A simplified model of the system allows the for-mulation of an optimal control problem that can be solved via dynamic programming. However, despite the small size of the model, numerical issues linked to discretization slow down the search for an appropriate parameterization of the dynamic programming algorithm. To increase the smoothness and the accuracy of the solution obtained using the state-of-the-art implementation of DP known as Level-Set DP [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF]), an improved version is developed and validated on a simpler problem with an analytical solution, before application to the ORC case study.

We start by recalling the principle of dynamic programming, its basic formulation and the level-set algorithm (Section 2). In Section 3, we present our improved variant of the level-set algorithm obtained via an adaptive discretization grid. The new algorithm is first validated on a simple energy management example, where the power consumption of an electrical vehicle is to be minimized (Section 4). At the expense of an acceptable increase of computational load, the adaptive-grid DP algorithm indeed provides a smoother and more accurate solution.

In Section 5, we present the optimal control application consisting in maximizing the power produced by an organic Rankine cycle waste heat recovery system on board a diesel-electric train. First, we describe the system. Then we provide an experimentally validated simplified model, and use it to formulate the optimal control problem. We show that, for this application, the adaptive-grid DP algorithm achieves the same accuracy of the reference level-set DP algorithm with fewer model calls, while using a lower level of discretization.

The analysis of the resulting optimal control trajectory provides useful insight for both control design and system design.

DYNAMIC PROGRAMMING AND LEVEL-SET METHOD

The generic optimal control problem we want to solve for a dynamic system with n state variables and m control inputs can be written:

min u(t) J (u (t)) (1) s.t. ẋ(t) = f (x(t), u(t), t) (2)
x(0) = x 0 (3) x(t) = X (t) ⊆ R n (4) u(t) = U(t) ⊆ R m (5) where J (u (t)) = g (x(t f )) + t f 0 h (x(t), u(t), t) dt, (6)
is the cost functional.

Since dynamic programming is a numerical algorithm, discretization of time, state space, and control space is required. Let the discrete-time model be x k+1 = f k (x k , u k ), k = 0, 1, ..., N -1 (7) where x k ∈ X k is the state vector and u k ∈ U k the control signal vector. Based on principle of optimality introduced by [START_REF] Bellman | Dynamic programming[END_REF], DP proceeds backwards in time to evaluate the optimal cost-to-go function J k (x i ) at every node x i of the discretized time-state space: I. Final time computation

J N (x i ) = g N (x i ), for x i ∈ X N ∞ else. ( 8
) II. Intermediate computations for k = N -1 to 0 J k (x i ) = min u k ∈U k {g k (x i , u k ) + J k+1 (f k (x i , u k ))}. (9)
The optimal control is the argument that minimizes the right-hand side of equation ( 9) for each x i at time index k, which yields the optimal control policy π = {µ 0 (x), µ 1 (x), ..., µ N -1 (x)}. This map is used to find the optimal control signal during a forward simulation of model ( 7) starting from a given initial state x 0 . Since the resulting state trajectory do not generally coincide with the nodes x i , interpolation is used.

An major issue to consider is the definition of the costto-go function for infeasible states. Grid points that are not backward-reachable should have infinite cost as in (8). However, this causes numerical problems as illustrated in Fig. 1. When running step k of backward iteration, the cost-to-go J k+1 (x) is known for all grid points x i . In order to evaluate the cost-to-go J k (x p ), all possible control candidates u ∈ U k are applied. In Fig. 1, the point f k (x p , u 1 ) is in the backward-reachable space. But since the value of J k+1 (x i+1 ) is infinite, interpolation will lead to consider that the cost-to-go of x i+2 k is infinite too.

A method to deal with this problem consists in using a large but finite value for infeasible states. The value of this penalty represents a critical parameter that is difficult to calibrate. Nevertheless, the penalty method results in a steep gradient of the cost-to-go function near the boundary and cannot completely solve the aforementioned problem.

Level-set algorithm

Another approach is based on the concept of backwardreachable space, as introduced in [START_REF] Back | Predictive powertrain control for hybrid electric vehicles[END_REF]; [START_REF] Sundström | On implementation of dynamic programming for optimal control problems with final state constraints[END_REF]. More particularly, the level-set DP algorithm, introduced by [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF], evaluates the backward-reachable space in parallel with cost-to-go function evaluation. A level-set function I is introduced that associates a real value to the state x:

I : X ⊆ R n → R, (10)
where I(x) ≤ 0 if x is in reachable space and I(x) > 0 else. The level-set function is evaluated for grid points x i at each time index k yielding a grid evaluation I k x . At the next backward iteration step k + 1, interpolation is used to determine whether a point x is backward reachable or not. Such a function can represent nonconvex regions, moreover, the number of state and control inputs is not limited. Assume that the final state constraint is given as a target set T , which is defined by a level-set function h(x)

h : X N → R, where X N ⊆ R n (11) T = {x ∈ X N |h(x) ≤ 0}. ( 12
)
Then a DP algorithm is applied where the cost-to-go function to minimize is the level-set function I:

I. Initialization I N (x i ) = h(x i ). (13) II. Backward iteration for k = N -1 to 0 I k (x i ) = min u k ∈U k {I k+1 (f k (x i , u k ))}. ( 14
)
3. ADAPTIVE-GRID DP ALGORITHM

Backward algorithm

The level-Set DP algorithm introduced in Elbert et al.

[2013] can be modified to improve the accuracy of solution.

I. Initialize k=N and the level-set and cost-to-go functions as:

I N (x i ) = h(x i ) (15) J N (x i ) = g N (x i ). ( 16
)
Considering that

x i ∈ R n let x i = x i (1) ... x i (n) .
A second discretization of state space XN is introduced:

XN = {x 1 N , x2 N , ..., xq N }, (17)
where each j th component of vector xi N is bounded below by x N (j) and above by xN (j) as follows:

x N (j) = min i∈{1,...,q}

{x i N (j) | 0 < I N (x i N )} (18) xN (j) = max i∈{1,...,q} {x i N (j) | 0 < I N (x i N )} (19)
As illustrated in Fig. 2 this second discretization may improve the accuracy of level-set and cost-to-go estimations:

ĨN (x i ) = h(x i ) (20) JN (x i ) = g N (x i ). ( 21
)
II. Reduce k by 1 and update the level-set function by

I k (x i ) = min u k ∈U k { Ĩk+1 (f k (x i , u k ))}. ( 22
)
A second discretization Xk is computed in the same way as for the first step. Then the level-set function is updated by Ĩk (x i ) = min

u k ∈U k { Ĩk+1 (f k (x i , u k ))}. ( 23
)
III. For each grid point xi , find the set of control signals for which the system trajectory ends up inside the backward-reachable space at the next time step 24) and the one control candidate that minimizes the level-set function *

U F k (x i ) = {u k ∈ U k | Ĩk+1 (f k (x i , u k )) ≤ 0}, (
u k (x i ) = arg min u k ∈U k { Ĩk+1 (f k (x i , u k ))}. ( 25
)
IV. Update the optimal cost-to-go by the following rule: if at least one valid control candidate is found, i.e., U F k (x i ) = ∅, then calculate the cost-to-go based upon the optimal candidate Jk (x i ) = min

u k ∈U F k (x i ) {g k (x i , u k ) + Jk+1 (f k (x i , u k ))}.
(26) If, however, the grid point is not backward-reachable, then calculate the cost-to-go based on the control input *

u k (x i ) Jk (x i ) = g k (x i , * u k ) + Jk+1 (f k (x i , * u k )) (27)
Fig. 2. A first estimation I k of level-set function is computed on the whole state space (left plot). A second discretization (right plot) allows to improve the accuracy of both level-set and cost-to-go functions.

Forward algorithm

To increase smoothness, the forward algorithm introduced in [START_REF] Elbert | Implementation of dynamic programming for n -dimensional optimal control problems with final state constraints[END_REF] is used with a new discretization of the control space:

Ũk = {u 1 k , ..., u r k U k , ..., u r k } (28)
I. Initialize k = 0 and x 0 at initial condition. II. Increase k by 1 and find the feasible control candidates

Ũ F k (x k ) = {u k ∈ Ũk | Ĩk+1 (f k (x k , u k )) ≤ 0}, (29) III. Find the optimal control input u 0 k (x k ) = arg min u k ∈U k {g k (x k , u k ) + Jk+1 (f k (x k , u k ))}
(30) IV. Simulate the system using the optimal control input

x k+1 = f k (x k , u o k ) (31) and repeat steps II. -IV. until k = N .
Note that the computational effort involved in the finer discretization Ũk is limited. Since steps II. -IV. of the algorithm are applied to one value of x k , the number of additional calls to the model is N • (r -r).

A SIMPLE ENERGY MANAGEMENT EXAMPLE

To illustrate the benefits of the new algorithm, a first optimal control problem with an analytic solution is studied. It concerns energy management for an electric vehicle.

System description

Neglecting the aerodynamic friction and considering the road slope null, the simplified continuous-time dynamics of the vehicle in the longitudinal direction x can be written as in [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF]:

ẋ = v, v = h 1 u -h 0 , (32)
where the control variable u is a percent torque demand, h 0 and h 1 are constant parameters depending on the rolling resistance coefficient, the motor maximum torque, vehicle mass, transmission ratio and wheel radius.

The on-board electric power consumption can be written as

P m = b 1 uv + b 2 u 2 (33)

Optimal control problem

We wish to find a control strategy that minimizes the power consumption under the constraints that the vehicle must reach a destination point at a distance D in a given time t f starting from a given point, at rest. In order to formulate this problem for DP solving, the continuous-time model is discretized using an Euler forward approximation with a time step T s . The optimal control problem can be formulated as follows:

min

u k ∈[-2, 2] N -1 k=0 b 1 u k v k + b 2 u 2 k T s ( 34
)
v k+1 = v k + (h 1 u k -h 0 ) T s ( 35
)
x k+1 = x k + v k T s ( 36
)
x 0 = 0, v 0 = 0 (37) 0 ≤ x k ≤ x, 0 ≤ v k ≤ v (38) D ≤ x N ≤ x, 0 ≤ v N ≤ v , (39)
where

N = t f
Ts . The final constraints are formulated in equations (39). In the following, a particular solution is sought for a final distance D = 200 and a final time t f = 60. The time discretization is chosen to be T s = 0.2s (then N = 300), with x = 205, v = 8, v = 0.3.

DP solution

The problem stated above is solved using the basic, the original level-set and the level-set with two-step discretization DP algorithms.

State variable spaces are first discretized respectively with N v = 81 and N x = 106, while the control space is discretized with N u = 81. Fig. 3 shows the resulting system trajectory and the corresponding control inputs for each of the three algorithms. They are compared to the analytic solution (a parabola) presented by [START_REF] Petit | Optimal drive of electric vehicles using an inversion-based trajectory generation approach[END_REF]. The results obtained using the new algorithm are closer to the analytic solution, even with a Fig. 3. Vehicle trajectory and control inputs when solving the problem with N v = 81, N x = 106 and N u = 81 relatively coarse discretization. In term of cost, the relative mean errors of the three algorithms are found to be of 1.4%, 0.78% and 0.19% respectively.

The evolution of backward-reachable space boundaries, where level-set function vanishes, helps to analyze those results. Fig. 4 shows that in the first backward steps, the optimal trajectory (x o , v o ) is close to Ĩk estimations. Here, for accurate estimation of boundaries it is important not to consider reachable points close to the optimal solution as infeasible. After those first backward steps, the optimal trajectory moves away from the boundaries. However, the reachable state space remain restricted. Thus, the second discretization step of the proposed algorithm can improve the accuracy of cost-to-go function. Two-step discretization involves a larger computational effort for a given set of parameters (N v , N x , N u ). This computation cost can be expressed in terms of number

OPTIMAL CONTROL OF AN ORC FOR WHR ON BOARD TRAINS

System description

The system under consideration is represented in Fig. 6. It is a prototype ORC system for WHR from a diesel engine generating set ("Power Pack") manufactured by Alstom Transport for installation in diesel-electric railcars.

A pump pressurizes an organic fluid (R365mfc) in liquid state and circulates it through a closed circuit. Via a bypass valve, a fraction V o,e of engine exhaust gas is fed to a heat exchanger (evaporator) which transfers heat to the working fluid. The vaporized fluid then enters a turbine, which converts some of the kinetic energy into work, driving a generator connected to the Power Pack electrical network. Lastly, the fluid is recondensed before returning to the pump at low pressure. Cooling is provided by a fan with controllable air flow * m A .

The objective here is to find the controls u = { * m A , V o,e } that maximize the net power P net produced by the WHR system, function of P exp , the power recovered by the turbine, and of P pump and P A , the powers consumed by the pump and the cooling system :

P net = P exp -P pump -P A , ( 40
)
The thermodynamic cycle of the working fluid is represented in Fig. 7. In the following, subscripts 1 to 4 will refer to the working fluid at components outlet, namely the evaporator, the expander, the condenser and the pump.

Recoverable mechanical power from the turbine is usually expressed [START_REF] Bao | A review of working fluid and expander selections for organic rankine cycle[END_REF]) in terms of a constant efficiency η exp , of the working fluid mass-flow * m and of the enthalpy gradient h 2 -h 1 :

P exp = * m(h 2 -h 1 ) η exp . ( 41
)
Similarly, the efficiency η pump can be considered constant and the power consumed by the pump can be written in terms of the pressure gradient and the fluid density:

P pump = * m p 1 -p 3 ρ 3 η pump . ( 42
)
The power consumed by the fan can be considered proportional to the delivered air flow [START_REF] Manente | An organic rankine cycle off-design model for the search of the optimal control strategy[END_REF]):

P A = * m A k A . ( 43
)
To solve this problem, we need a model to predict the evolution of the variables in the expressions of powers ( 41)-( 43), that is fluid mass-flow, pressures and enthalpies.

Model dimensions should be small so as not to be confronted with the well-known "curse of dimensionality" [START_REF] Powell | Approximate Dynamic Programming: Solving the curses of dimensionality[END_REF]). The so-called moving boundaries (MB) method allows a realistic dynamic representation of heat exchangers with a limited number of states. [START_REF] Cheng | Nonlinear observer design for two-phase flow heat exchangers of air conditioning systems[END_REF] proposes a model that considers an average wall temperature, while [START_REF] Peralez | Improving the control performance of an organic rankine cycle system for waste heat recovery from a heavy-duty diesel engine using a model-based approach[END_REF] shows experimentally that, for an ORC, the slowest dynamics are those of the wall temperatures. Thus, it is possible to neglect the fluid dynamics, thereby reducing the number of states.

In the next paragraph, we propose a new heat exchanger model with a single state variable obtained by applying the first law of thermodynamics. This model is then validated against experimental data, showing that it is able to capture the slowest dynamics of the ORC.

Fig. 6. ORC system under investigation

Two-state ORC model. Considering the thermodynamic cycle in Fig. 7, fluid superheating at evaporator outlet (point 1) is assumed to be perfectly regulated at a constant value by the pump mass-flow. The fluid at pump inlet (point 3) is assumed to be in the saturated liquid state. The slow dynamics of the system, i.e. the thermal dynamics of exchangers walls, are taken into account while the fluid is considered to be in equilibrium.

Let Tw,e and Tw,c the average wall temperatures respectively of the evaporator and the condenser. Let Tf,e and Tf,c the average temperatures of the working fluid circulating respectively in the evaporator and the condenser. The wall energy balance can be written as:

+ η A * m A cp A (T A -Tw,c ), (44)
where

η A = 1 -exp(- α A S A * m A c A
) reflects the decrease in heat exchange efficiency for large air mass flows [START_REF] Mckinley | An advanced nonlinear switched heat exchanger model for vapor compression cycles using the moving-boundary method[END_REF]).

In the following, Tf,e and Tf,c will be assumed equal to the evaporation temperatures respectively of the high-pressure and of the low-pressure sides of the ORC. Tf,e (resp. Tf,c ) is then an unimodal function of p 1 (resp. of p 3 ). Likewise, h 1 , ρ 1 depend only on p 1 , whereas h 3 depends on p 3 .

Fluid energy balances in the exchangers yield:

       0 = h 1 -S f,e ᾱf,e
Tw,e -Tf,e * m -h 3

0 = h 3 -S f,c ᾱf,c Tw,c -Tf,c * m -h 2 , (45)
where the mass flow * m is considered homogeneous throughout the circuit. The mass flow through the turbine is given by *

m = k 2 ρ 1 p 1 , (46)
where k is a constant (function of the equivalent section of the turbine nozzles).

Assuming a constant isentropic efficiency η is of the turbine

h 2 = h 1 -η is (h 1 -h is,2 ), (47)
where h is,2 corresponds to the enthalpy whose entropy is equal to that of point 1 and pressure is equal to that of point 3. h is,2 then depends on p 1 and p 3 .

Using the previous expressions for * m and h 2 , system (45) can be written in the following implicit form: 0 = Φ( Tw,e , Tw,c , p 1 , p 3 ). ( 48)

Knowing Tw,e and Tw,c , ( 48) is solved numerically by iterating over p 1 and p 3 (as illustrated in Fig. 8). Tf,e and Tf,c are then deduced from p 1 and p 3 to express the wall thermal dynamics (44).

Optimal control problem

We wish to find a control strategy that maximizes the power production P net , defined by equations ( 40)-( 43), along a given mission of duration t f . The mission profile considered herein, a portion of a typical intercity train trip, results in a set of external conditions for the ORC system, namely for the temperature T exh (t) and mass-flow * m exh (t) of engine exhaust gas. A variation of ambient temperature T A (t) is also considered as shown in Fig. 11.

The cost functional to be maximized depends on time, state variables x = { Tw,e , Tw,c } and control inputs u = {V o,e , * m A } :

J = t f 0 P exp (x) + P pump (x) + P A (u, t), dt, (49)
where the dynamic of x = { Tw,e , Tw,c } is defined by the two-state ORC model equations ( 44)-( 48). This maximization problem must be solved under a security constraint on pressure p 1 . A final constraint on wall temperature x 1 is also considered to let the system stop at the end of the mission:

p 1 (x, t) < 25 bar, ∀t ∈ [0, t f ] (50) x 1 (t f ) < 130 • C (51)
Bounds on control inputs are defined by technical considerations. The fraction of exhaust gas is bounded by [0, 1] whereas air mass flow provided by fan is comprised between 0 and 4 kg s -1 :

u(t) ∈ [0, 1] × [0, 4], ∀t ∈ [0, t f ] (52)

DP solution

The optimal control problem is solved using the original level-set DP and the new adaptive-grid DP algorithms. Different levels of state space and input variables discretization are tested with numbers of grid points set to

N = N x1 = N x2 = N u1 = N u2
. The computational efforts -in terms of number of model evaluations -required to achieve a similar level of accuracy are compared in Table 1.

The comparison shows that, like in the simple thermal management example, a coarse discretization suffices to obtain good results using the adaptive-grid algorithm.

Results analysis

Fig. 12 shows the system trajectory obtained with the discretization corresponding to N = 21. Wall temperatures of heat exchangers periodically rise in response to the increase of the exhaust gas mass flow. This yields a pressure increase (top plot) activating the safety constraint (50). A the end of the mission, the final constraint (51) leads to the complete closing of the evaporator by-pass, i.e u 1 = 0.

The influence of ambient temperature on the overall efficiency of the ORC system is also highlighted. In fact, the power recovered by the turbine is substantially reduced for high ambient temperatures, whereas the powers consumed by the pump and the cooling fan are scarcely affected.

These results can be very useful to design a control strategy. We can infer that tracking control of condenser pressure can be obtained using the air mass flow * m A to regulate condenser cooling; which would in turn ensure the pressure ratio required to run the turbine efficiently. On the other hand, the evaporator by-pass can be used to enforce safety limits in terms of pressure and temperature.

Solving the optimal control problem introduced above should also be useful at system design stage, especially considering that transient behavior is accounted for in system performance evaluation. The problem could be formulated for variants of the system under consideration, obtained with different sets of model parameters. For instance, the influence of working fluid choice or of exchangers sizing -taken into account in the two-state model described in Section 5.1 -could be studied. Thanks to the good accuracy achievable by the adaptive-grid DP algorithm at low levels of discretization, it seems feasible to study the impact of these (numerous) parameters within reasonable simulation times.

CONCLUSION

This paper proposes two main contributions.

• An original optimal control problem is formulated in order to find the maximum power produced by an organic Rankine cycle waste heat recovery system on board a diesel-electric train, along a given mission profile. An ad-hoc simplified model is developed, and experimentally validated, to make solving by dynamic programming computationally treatable. The optimal control trajectories obtained provide valuable insight into system behavior and prove useful for control design.

• An improved version of the state-of-the-art level-set DP algorithm is developed based on an adaptive grid, and validated on a simpler optimal control problem with an analytical solution, in order to achieve the required accuracy with a lower level of discretization and shorter computational time than the original algorithm.

Since the proposed improvement does not rely on a specific structure, application of the new algorithm should also be beneficial to other optimal control problems. Yet, we have not proven (nor tried to prove) that our "adaptivegrid" version consistently outperforms the original level- Nonetheless, the accuracy obtained on the ORC study case with a relatively short computational time opens up new opportunities. For instance, it becomes possible to carry out a thorough parametric study of the ORC system, in order to optimize its design by taking into account the transient behavior.

Fig. 1 .

 1 Fig. 1. Numerical issues near boundary surface 2.1 Basic DP

Fig. 4 .

 4 Fig. 4. Optimal trajectory (x o , v o ) and evolution of backward-reachable space (evaluated by adaptive-grid DP).

Fig. 5 .

 5 Fig. 5. Algorithm accuracies vs number of model calls for the simple energy management example. of model calls for evaluation. For example, for the level of discretization used above, basic and original level-set DP both require about 2.085 • 10 8 model calls whereas the new algorithm requires about 2.85 • 10 8 . For a fair comparison, algorithm accuracies are compared for different levels of discretization. Results are shown in Fig. 5. The adaptivegrid DP algorithm attains acceptable accuracies for a much lower number of calls.

M

  w,e Ṫw,e = S f,e ᾱf,e ( Tf,e -Tw,e ) + V o,e * m exh cp exh (T exh -Tw,e ) M w,c Ṫw,c = S f,c ᾱf,c ( Tf,c -Tw,c )

Fig. 7 .

 7 Fig. 7. Rankine cycle for "dry" fluid

Fig. 9 .

 9 Fig. 9. External conditions for experimental validation of one-state evaporator model.

Fig. 10 .

 10 Fig. 10. One-state evaporator model: experimental validation

Fig. 11 .

 11 Fig. 11. Input disturbances for ORC problem.

Fig. 12 .

 12 Fig. 12. Adaptive-grid DP solution: time evolution of the ORC system characteristics set algorithm in different contexts than those presented in the paper.

Table 1 .

 1 Comparison of computational efforts required for the ORC problem.

		Level-set DP Adaptive-grid DP
	Level of discretization	N = 24	N = 14
	Mean net power produced	5.485 kW	5.497 kW
	Number of model calls	8.309 • 10 7	1.24 • 10 7

The authors wish to thank ANR, the French National Research Agency, for partially funding this work within the framework of the project TRENERGY ("Train Energy Efficiency via Rankine-cycle Exhaust Gas Heat Recovery", # ANR-12-VPTT-09), and ALSTOM Transport, for providing application data.


Paulo E U Souza

Caroline P Carvalho

Frédéric Dehais

MOMDP-based target search mission taking into account the human operator's cognitive state

Keywords: MOMDP, robotics, mixed initiative, operators cognitive state estimation

This study discusses the application of sequential decision making under uncertainty and mixed observability in a mixed-initiative robotic target search application. In such a robotic mission, two agents, a ground robot and a human operator, must collaborate to reach a common goal using, each in turn, their recognized skills. The originality of the work relies in considering that the human operator is not a providential agent when the robot fails. Using the data from previous experiments, a Mixed Observability Markov Decision Process (MOMDP) model was designed, which allows to consider aleatory failure events and the partial observable human operator's state while planning for a long-term horizon. Results show that the collaborative system was in general able to successfully complete or terminate the mission, even when many simultaneous sensors, devices and operators failures happened. So, the mixed-initiative framework highlighted in this study shows the relevancy of taking into account the cognitive state of the operator, which permits to compute a policy for the sequential decision problem which prevents to re-planning when unexpected (but known) events occurs.

I. INTRODUCTION

Unmanned Vehicles (UVs) are becoming increasingly present in a wide variety of operational contexts such as military operation, border security, inspection of contaminated area for prevent human from hazard exposure. Most of scientific and technical efforts have focused on the implementation of smart sensors, complex embedded systems and autonomy to enhance the efficiency of the UVs [START_REF] Thrun | Autonomous exploration and mapping of abandoned mines[END_REF], especially when the human operator can not analyze or access visual data [START_REF] Thrun | Autonomous exploration and mapping of abandoned mines[END_REF]- [START_REF] Casper | Human-robot interactions during the robot-assisted urban search and rescue response at the world trade center[END_REF]. However these developments were generally achieved without questioning the integration of the human operators in the control loop [START_REF] Schurr | Impact of human advice on agent teams: A preliminary report[END_REF]: the human operator is considered as a providential agent that will be able to take over when sensors or automations fail [START_REF] Schurr | Improving adjustable autonomy strategies for time-critical domains[END_REF]- [START_REF] Casper | Human-robot interactions during the robot-assisted urban search and rescue response at the world trade center[END_REF]. Yet, poor user interface design, complexity of automation and high operational pressure can leave the human operator ill-equipped when mental workload exceeds human capacity [START_REF] Durantin | Using near infrared spectroscopy and heart rate variability to detect mental overload[END_REF]. For instance, careless design of authority sharing can lead to human-automation conflicts when the human operator misunderstand the automation behavior [START_REF] Dehais | Towards an anticipatory agent to help pilots[END_REF], [START_REF] Dehais | Automation surprise in aviation: Real-time solutions[END_REF]. The occurrence of such situation is critical as long as it may cause "mental confusion" (i.e. the human operator is unable to glance and process the relevant parameters) [START_REF] Dehais | Automation surprise in aviation: Real-time solutions[END_REF] or attentional tunneling (i.e. the human operator is excessively focused on a single display) [START_REF] Dehais | Mitigation of conflicts with automation use of cognitive countermeasures[END_REF] yielding to irrational behavior [START_REF] Dehais | Cognitive conflict in human-automation interactions: a psychophysiological study[END_REF]. Not surprisingly, a safety analysis report [START_REF] Williams | A summary of unmanned aircraft accident/incident data: Human factors implications[END_REF] revealed that human factors issues were involved in 80% of accidents. This trend has led Cummings and Mitchell (2008) to state: "Because of the increased number of sensors, the volume of information, and the operational demands that will naturally occur in a multiple-vehicle control environment, excessive cognitive demands will likely be placed on operators. As a result, efficiently allocating attention between a set of dynamic tasks will be critical to both human and system performance. -p. 451".

A promising avenue to deal with these issues is to consider that robot and human abilities are complementary and are likely to provide better performance when joined efficiently than when used separately. This approach, known as mixed-initiative [START_REF] Allen | Mixed-initiative interaction[END_REF], [START_REF] Ai-Chang | Mapgen: mixed-initiative planning and scheduling for the mars exploration rover mission[END_REF] defines the role of the human and artificial agents according to their recognized skills. It allows the human and the robot to set the appropriate level of autonomy of the robot [START_REF] Hardin | On using mixed-initiative control: a perspective for managing large-scale robotic teams[END_REF]. An interesting mixedinitiative proposition, presented by [START_REF] Sellner | Coordinated multiagent teams and sliding autonomy for large-scale assembly[END_REF], relies on a statistical approach to determine which entity (i.e human or UVs) is the most efficient for a given task. Interestingly enough, this approach paves the way for allocating roles and sharing authority between the human and artificial agents. In [START_REF] Finzi | A mixed-initiative approach to human-robot interaction in rescue scenarios[END_REF], a mixed-initiative planning approach is proposed to monitoring the system and to coordinate operator's interventions in a rescue scenario. In this work, the mixed-initiative planning continuously coordinates, integrates and monitors the operator's interventions and decisions. Another example can be found in [START_REF] Gianni | Awareness in mixed initiative planning[END_REF], in which a robot and an human operator collaborate for an urban search and rescue mission in order to detect and report objects of interest.

A key issue to design a mixed-initiative system is to implement a decision system. This latter defines the role and the authority of human and artificial agents, while estimating capabilities of evolved human (intention, situation awareness, sensor's failure perception) and robotic agent (sensor's status, mission task, etc). Such decision-making system can be governed by a policy resulting from the resolution of a Partially Observable Markov Decision Process (POMDP), as proposed by [START_REF] Taha | A POMDP framework for modelling human interaction with assistive robots[END_REF], which is able to adapt itself to the user's intention getting feedback from the user in terms of satisfaction. A different way to drive interaction using POMDPs is studied in [START_REF] Hoey | Assisting persons with dementia during handwashing using a partially observable markov decision process[END_REF] for assisting persons with dementia during hand-washing. Note that, the state vector of the robot can be often considered as fully observable while the operator's cognitive state is, by definition, partial observable. Such decomposition can be addressed using a Mixed Observability Markov Decision Process (MOMDP) [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF], which is a stochastic model derived from the POMDP [START_REF] Smallwood | The optimal control of partially observable Markov processes over a finite horizon[END_REF]. The MOMDP is a formal framework that considers fully and partially observable state variables under probabilistic uncertainties while decreasing the computational cost required to produce a optimal policy [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF]. In addition, these two types of agents may face unexpected random situations during the mission. It can be modeled as probabilistic effects of actions. Moreover, this kind of model allows the inclusion of the uncertainty in the observations of the agents' states (i.e the cognitive state of the human operator) and the environment. The MOMDP aim to achieve a policy that maps an optimal action for each belief state -composed by the observable state and the partially observable state estimation. Thus, it is expected that the resulting policy could help to implement a genuine adaptive interaction, because this formalism is perfectly suited to maintain a state estimation and to decide of the men-robot system dynamics based on data coming from sensors applied to the operator (e.g eye-tracker) and from sensors embedded in robots.

In this present study, we propose to test the MOMDP approach on an mission involving a human and a physical UV that cooperate to perform a target identification task. Data collected during previous experiments allowed us to set probabilities of UV failure as well as of human operators poor cognitive state. This paper is organized as follow: first we recall POMDP and MOMDP models. In the sequence we present the mission model treated. Afterwards we evaluate the results obtained for such modeling. And, finally we conclude and discuss future work.

II. BACKGROUND

A. POMDP overview

POMDPs model situations where the agent only has access to partial information about the state of the system. A POMDP is a Markov Decision Process where the agent does not have access to the state of the system: it has only a partial and imprecise observation [START_REF] Smallwood | The optimal control of partially observable Markov processes over a finite horizon[END_REF]. In this context, the agent maintains a probability distribution over states, i.e. a belief state, which is updated after each action executed and observation perceived.

A POMDP is a tuple (S, A, Ω, T, O, R, b 0 , γ) where:

• S is a bounded set of states; (b) MOMDP transition model Figure 1. Transition models of a POMDP and a MOMDP.

• A is a bounded set of actions;

• Ω is a bounded set of observations;

• T : S × A × S → [0; 1] is a transition function such that T (s t+1 , a, s t ) = p(s t+1 | a, s t ); • O : Ω × S → [0; 1] is an observation function such that O(o t , s t ) = p(o t |s t ); • R : S × A → R
is a reward function associated with a state-action pair, and; • b 0 is the initial probability distribution over states.

• γ ∈ [0, 1[ is the discount factor We note ∆ the belief state space. At each time step t, the agent updates its belief state defined as an element b t ∈ ∆ using the Bayes' rule [START_REF] Smallwood | The optimal control of partially observable Markov processes over a finite horizon[END_REF].

Solving a POMDP consists in finding a policy function π : ∆ → A that maps to each belief state an optimal action that maximizes a performance criterion. The expected discounted reward from any initial belief

V π (b) = E π [ ∞ t=0 γ t r(b t , π(b t )) | b 0 = b]
is usually optimized. The value of an optimal policy π * is defined by the optimal value function V * that satisfies the Bellman optimality equation:

V * (b) = max a∈A r(b, a) + γ o∈Ω p(o|a, b)V * (b o a ) (2)
where, γ is the discount factor. When r(b, a), can be computed as an average gain r(b, a) = s r(s, a)b(s), the optimal value of belief states is proven to be piecewise linear and convex and solution of the Bellman's equation [START_REF] Smallwood | The optimal control of partially observable Markov processes over a finite horizon[END_REF]. As it, at n th optimization stage, the value function V n can be parametrized as a set of hyperplanes over ∆ known as α-vectors. An α-vector and the associated action a(α i n ) define a region of the belief state space for which this vector maximizes V n . Thus, the value of a belief b can be defined as:

V n (b) = max α i n ∈Γn b • α i n .
The optimal policy at this step is then:

π n (b) = a(α b n ).
Recent offline solving algorithms HSVI2 [START_REF] Smith | Point-based POMDP algorithms: Improved analysis and implementation[END_REF] or SARSOP [START_REF] Kurniawati | SARSOP: Efficient point-based POMDP planning by approximating optimally reachable belief spaces[END_REF], for instance, approximate the value function with a bounded set of belief states B, where B ⊂ ∆. These algorithms implement different heuristics to explore the belief state space using probabilistic trials reaching in this way only relevant belief states, and updating the value of V for them, instead of computing the value function for all the belief state space, which is a continuous space.

B. MOMDP

The Mixed Observability Markov Decision Process (MOMDP) is an extension recently proposed for the POMDP model [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF], which explores the particular structure where certain state variables are fully observable. This factored model leads to a very significant time gain in policy computation, improving the efficiency of a point-based algorithms. According to [START_REF] Ong | POMDPs for robotic tasks with mixed observability[END_REF] the completely observable state is represented by x and the partially observable state by y. In this way, the couple (x, y) specifies the complete state with |S| = |X | × |Y|, where X represents the space with all the possible values of the variable x (resp. Y to y).

A MOMDP is a tuple (X , Y, A, Ω, T X , T Y , Ω, R, b 0 , γ), where:

• X is the bounded set of fully observable state variables;

• Y is the bounded set of partially observable state variables; • A is a bounded set of actions;

• Ω is a bounded set of observations;

• T X : X × A × X × Y → [0; 1] is a transition function such that T X (x, y, a, x ) = p(x |x, y, a); • T Y : Y ×X ×A×X ×Y → [0; 1] is a transition function such that T Y (x, y, a, x , y ) = p(y |x, y, a, x ); • O : Ω × Y → [0; 1] is an observation function such that O(o, a, x , y ) = p(o|x , y , a); • R : X × Y × A → R is a reward function associated
with a state-action pair; and:

• b 0 = (x 0 , b Y0 ) is the initial probability distribution over states. • γ ∈ [0, 1[ is the discount factor.
Note that, as the probability distribution over states concerns only the Y set, the belief state update is redefined as:

b o,a,x Y (y ) = η y ∈Y p(o|y , x , a)p(y |x, y, a, x )p(x |x, y, a)bY (y) (3)
where, η is a normalization constant. The belief state b is now noted by the couple (x, b Y ), and B Y is the belief state space y conditioned by x :

B Y (x) = {(x, b Y ), b Y ∈ B Y }. B Y (x) is a sub-space of B, such that B = x∈X B Y (x).
Solving MOMDPs consists in finding a set of policies π x : B Y → A, which maximize the criterion :

π * x ← arg max πx∈Π Eπ x ∞ t=0 γ t r((xt, b t Y ), π((xt, b t Y ))) b0 = (x0, bY 0 ) (4)
As for the POMDP, the value function at a time step n < ∞ can be also represented by a set of α-vectors:

V n (x, b Y ) = max α∈Γ n Y (x) (α • b Y ) ( 5
)
where α is the hyperplan over the space B y (x). In this way, the value function over the complete state space is parametrized by the set

Γ Y (x), i.e. Γ = {Γ Y (x), x ∈ X }.
So, given a belief state (x, b Y ) the optimal action is defined by the action associated with the α-vector that maximizes

max α∈Γ Y (x) (α • b Y ).
For more details about MOMDP algorithm resolution, please see [START_REF] Ong | Planning under uncertainty for robotic tasks with mixed observability[END_REF], [START_REF] Ong | POMDPs for robotic tasks with mixed observability[END_REF].

Next, we present previous experiments which were used as base for statistical data in order to leaning the MOMDP model for the target search mission taking into account the operator's cognitive state.

III. PREVIOUS EXPERIMENTS A. Material

The experimental set-up was developed at ISAE-SUPAERO. It was composed of a robot and a ground station. The robot was equipped with different sensors such as a GPS for autonomous navigation, an Ultrasound sensor to detect and avoid obstacle, a video camera and an Xbee transmitter to communicate with the ground station. It had a 15 minutes autonomy thanks to electrical battery. The robot could be operated in manual mode or in supervised mode. In manual mode, the robot was operated with a joystick. In supervised mode, the robot performed waypoint navigation autonomously, but any action of the operator with the joystick let her/him take over until the joystick was released. The ground station was displayed on a 24inch screen showing different kinds of information to control and supervise the robot such as a tactical map, a panoramic video scene screen; a mission synoptic ; an interactive panel sending the requests to the human operator; a status panel panel indicating the state of the GPS, the ultrasound sensor and the battery; and a guidance mode state (i.e. supervised or manual). Note that the operator could not see the robot and only gathered information through the screen. Fig. 2 shows the interface of the ground station to operate the robot during the mission.

B. Assessing probability of failures of the robot

A first experiment was conducted to assess the probability of failure of the different sensors and devices embedded in the robot. Thirty tests were run and consisted of a 10 minutes navigation task with the robot evolving among four waypoints in a field at ISAE-SUPAERO campus. The results are summarized on Table I, where FP design the failure probability.

C. Assessing human's performance

Data collected from a previous experiment [START_REF] Dehais | Cognitive conflict in human-automation interactions: a psychophysiological study[END_REF] were used to assessing the probability of the human operator to perceive robot failures. The scenario of this experiment consisted of a target localization and identification task involving a ground robot and a human operator. The target to be identified had two short messages written in white on each side (front side OK, back side KO). The mission had four main segments: S 1 -Reach the area, S 2 -Scan for target, S 3 -Identify target, and S 4 -Battery-Failure. At the beginning of the mission, the robot evolved in an autonomous supervised mode to reach the target area (S 1 ) and then start scanning the area to detect the target (S 2 ). When the robot was close to the target, the operator had to take over in manual mode and to identify the target S 3 ). While the operator was performing the identification task, a low battery event was triggered (S 4 ). In turn this event yield to a safety procedure that made the robot to go back to base autonomously. As this event occurred while the operator was excessively focused on his target identification task, it was expected that he would missed the alerts and thus persist in achieving the target detection task.

1) Assessing failure perception: 12 subjects participated to the experiment and were equipped with an electrocardiogram (ECG) and a 25 Hz Pertech head mounted eye tracker. This latter device was used to collect participants' eye gaze on the user interface. More specifically we focused our eye tracking analysis on eight areas of interest (AOI) of the user interface: 1) tactical map, 2) message panel, 3) guidance mode (supervised vs manual), 4) synoptic, 5) "back to base" warning, 6) GPS and ultrasound status, 7) battery status, 8) panoramic video. The collected ocular data were used to set the probability of the operator to perceive the sensor's status (sensors statuses are summarized in Table I). This sensor status perception probability is based on the normalized sum of the averaged fixation time (∆T ) on the related AOIs. For instance, when the GPS or the ultrasound are lost, the icons turns to red (area 6) and the robot is stopped (i.e it can be seen through the panoramic video -area 8). When the low-battery event occurs, three changes can be observed in the user interface: (i) the battery icon (area 7) turns to orange with the associated low battery message, (ii) the mode changes automatically from manual to supervised, and area 3 blinks twice and (iii) the segment status became back to base (area 5).

Thus we introduced the spGpsUs boolean state variable that can be used to model perception about a GPS or ultrasound status by the operator: With the same reasoning, for the spBatt (Battery status perception) boolean state variable, the transition probability was defined by the normalized sum of the averaged time that the participants expended looking to areas 3, 5 and 7 during the manual and autonomous operations: 2) Assessing cognitive availability: the result of the experiment revealed that 8 participants out of 12 did not understand the robot behavior, though some of them glanced at the battery failure indicator. These 8 participant persevered to achieve the no-longer-relevant identification task [START_REF] Dehais | Cognitive conflict in human-automation interactions: a psychophysiological study[END_REF]. This typical behavior is known as "attentional tunneling" and is defined as "the allocation of attention to a particular channel of information, diagnostic hypothesis or task goal, for a duration that is longer than optimal, given the expected costs of neglecting events on other channels, failing to consider other hypotheses, or failing to perform other tasks" [START_REF] Thomas | Eye-tracking and individual differences in off-normal event detection when flying with a synthetic vision system display[END_REF]. Therefore, the inference of such impaired attentional state is of great importance to design a mixed initiative system. We implemented an Adaptive Neuro-Fuzzy Inference System (ANFIS) to detect attentional tunneling that is associated with higher cardiac activity, decreased saccadic activity and long concentrated eye fixations [START_REF] Régis | Formal detection of attentional tunneling in human operator-automation interactions[END_REF]. The ANFIS classifier had a probability of 91.1% to detect Attention Tunneling (please report to [START_REF] Régis | Formal detection of attentional tunneling in human operator-automation interactions[END_REF] for more details). This detection probability was used in this study to define the observation function of the Cognitive Availability state variable as shown further. Cognitive Availability is defined here as the capability of the human operator to be available and aware of the robot's status during all mission tasks. 3) Cognitive countermeasure to assist the human operator: many studies revealed that alarms are inefficient to warn human operator during high workload situations such as performing manual control and identifying target [START_REF] Dehais | Failure to detect critical auditory alerts in the cockpit evidence for inattentional deafness[END_REF]. Rather than adding alarms during stressful situations, an optimal solution to warn the human operator consists of temporarily removing the information the human operator is focusing on, and replacing it by an explicit visual stimulus designed to change the attentional focus. The principle of this cognitive countermeasure was tested in second experiment with 11 participants facing the same scenario (i.e. target identification task and battery failure). The results revealed that these cognitive countermeasures helped 10 participants out of 11 to perceive the battery failure and to let the robot go back to base [START_REF] Dehais | Mitigation of conflicts with automation use of cognitive countermeasures[END_REF].

IV. MODELING THE COLLABORATIVE TARGET IDENTIFICATION MISSION

Using all those previous experimental data, a MOMDP model was defined in order to drive the adaptive interaction between the human operator and the UV. The choice for a Mixed Observability model comes from the nature of our problem: the robot and mission states can be considered as fully observable, while the operator's cognitive ability, here considered as Cognitive Availability is a partially observable state variable by definition.

The mission can be decomposed in six high level phases: going to zone, searching the target, handling the target, returning to base, on base and failed. The robot states can be defined by the cross product of the phases, the embedded sensors statuses, the statuses of the ground station, the on board camera, the Xbee and the joystick devices, and the Cognitive Availability as cognitive state of the operator. Next, we present the fully and partially observable state variables considered in the model.

Fully observable state variables (X ): The section III-C and tables I and III, present the fully observable state variables considered in the mission modeling. As discussed before, mission phases were classified according to the operation mode. The sensors statuses were discretized in two possibilities: OK and KO (not OK). It was also assumed that after a sensor failed, it switches to KO and remained KO until the end of the mission. The sensors' failure probabilities were shown in Table I.

Fig. 3 summarizes the transition function for the mission phase state variable. A manual mode was associated with each autonomous mode (except the on base and failed mission phases). One can argue that, in a human operator's point of view, there is only one manual mode, but for modeling propose, the manual mode was factored in four phases (see Table III) to prevent the planner from selecting a supervised mode already held when returning to the autonomous operation mode. As show in previous experiments presented before, it is relevant to define two fully observable variables that model the operator's sensor's status perception (see section III-C). Note that, operator's perception about sensor's status (GPS/ultrasound and battery) state variables are assumed as fully observable, because it is not possible to observe if the human operator perceived (i.e. his cognitive process) only with the eye-tracker device. In this case, for example, if the operator looked to the areas (6 and 8) between two decision time steps, he should detect, with a probability of 0.70, for supervised mode, or 0.86, for manual mode, if there was a GPS or a ultrasound breakdown because the related icon turns to orange (see Fig. 2 and Section III-C).

Partially observable state variable (Y): The operator's Cognitive Availability is considered in this study as the opposite of Attentional Tunneling [START_REF] Dehais | Mitigation of conflicts with automation use of cognitive countermeasures[END_REF], [START_REF] Régis | Formal detection of attentional tunneling in human operator-automation interactions[END_REF], [START_REF] Pizziol | Towards human operator state assessment[END_REF]. The measure of the allocation of attention if not a straightforward task [START_REF] Régis | Formal detection of attentional tunneling in human operator-automation interactions[END_REF]. Therefore, we consider the Cognitive Availability of the human operator as a partially observable variable. Hence: available Y models that the human operator is cognitively available (resp. available N, not cognitively available). Associated with this partially observable state variable we have two possible observations: oAvailable Y meaning that the operator is observed as cognitively available and potentially aware of the situation and oAvailable N modeling that he is observed as not cognitively available. Table II summarizes the observation probability function for this observation variable.

Actions: Discrete actions were defined as: goToZone, tgtSearch, tgtHandle, retBase, onBase, getAttention and cMeasure. Action result depends on the aleatory sensors behavior (cf. Table I). For instance, in a nominal case and based on previous works [START_REF] Dehais | Mitigation of conflicts with automation use of cognitive countermeasures[END_REF], [START_REF] Régis | Formal detection of attentional tunneling in human operator-automation interactions[END_REF], [START_REF] Pizziol | Towards human operator state assessment[END_REF], the robot is able to autonomously navigate and avoid obstacles, but if the robot chosen goToZone and the ultrasound sensor fails, the mission phase turns to manual mode (m0) (see Fig. 3) because the robot is no more able to avoid obstacles autonomously. If a low battery event arrives, the robot can return to the base automatically if it was in a supervised mode. When it was in a manual mode, it can switch to returning to base (with any action) automatically only if the human operator is aware of the failure, i.e. if he was observed as cognitively available being aware of the situation and by consequence leaving the joystick.

The getAttention is considered as a non deterministic action, since it should be used when the robot needed help and the operator's Cognitive Availability was estimated as "not available" (oAvailable N). The same occurs with the action cMeasure (countermeasure), which should be executed when a low battery event arrives during a manual operation and the operator was considered as "not available" (e.g his attention was focused on handling the robot and he would not notice the alerts on the user interface). In such case, when a cMeasure action is launched the robot wait the human operator leaves the joystick (see Fig. 3).

Rewards: The reward function (R) was designed in order to payoff suitable actions, for instance, goToZone in the phase e0 when the navigation sensors are OK, and to punish otherwise. The same occurs with the manual modes and its essential devices (cf. Table I). Note that, we have chosen to associate a increasing reward with sequential phases, i.e. reward associated with the action tgtSearch in e1 considering essential sensors are OK (R=15) is more important than the action goToZone in the e0 phase (R=10). We have considered that processing the target in manual mode is more dependable than in autonomous mode, since the interpretation done by the human operator is more reliable. In this case, the reward for the choice of tgtHandle in manual mode (R=30) is more important than tgtHandle in supervised mode (R=20). The action getAttention only has a positive payoff (R=30) if the operator's Cognitive Availability is estimated as "not available" (oAvailable N) considering that at least one of the essential devices are KO and that the human operator did not see the alert, otherwise the reward is negative (-500). Similarly, the action cMeasure has a positive payoff (R=50) only if the operator is perceived as "not available" in a manual mode when a low battery event arrives.

A mission is considered as fully accomplished if the robot had passed through the phases e2 or m2 (resp. processed the target autonomously or manually) and arrived at base e4. When the robot returns to the base (autonomously or manually) before processing the target, the mission is considered as aborted and when the robot is unable to reach the base, the mission is considered as failed.

V. SIMULATION RESULTS

The APPL 0.96win SARSOP1 [START_REF] Ong | POMDPs for robotic tasks with mixed observability[END_REF] was used as solver. The grammar of this solver has a special format that differs from the classical input POMDP file format 2 . Therefore, we have developed a script written in Python 2.7.8 to produce the MOMDP input file. For the MOMDP resolution, we have set the precision to 0.01. We recall that this precision is related to the difference between the upper and lower bound of the value function for the initial belief state, which is considered as a stop optimization criterion.

Statistical analysis was performed to process the results over 5000 policy simulations. Fig. 4(a) provides an overview of how many times the robot passed by each mission phase. Note that, the robot passed through phase e1 exactly 5000 times and never crossed the phase m0, this can be explained by the fact that the initial state, when all fully observable state variable values are known, caused a deterministic cycle in the first time stamp.

To sum up, the mission was fully accomplished 2998 times (59.96%) (cf. Fig. 4(b)), meaning that the robot has passed through the phases e2 or m2, respectively processing the target autonomously or manually, and arrived at base e4. In such cases, the target was handled in autonomous mode (e2) 1228 times, which represents 41% of successful missions, and it was handled in manual mode (m2) 1770 times (resp. 59% of successful missions). Fig. 4(b) also shows that the robot returned to the base in 68.72% of times (including aborted missions). The mission completely failed, i.e. it reached e5, 1564 times (31.28%).

The Table IV presents an example of a Fully accomplished mission where the robot changed its mode to manual (m2) for the operator process the target (bigger reward), then the GPS failed. This was not a problem at that moment because the operator did not need the GPS to handle the target. Next, the robot remained in manual (m3) but the operator seemed not to be aware, so, the robot ask for his attention (action getAttention), and the operator led the robot to the base.

Finally, a aborted mission is shown in the table V. In this aborted mission a low battery event occurred while the operator was processing the target and the robot observed him or her as "not available" (not aware) of the failure. Consequently, the robot executed a countermeasure action (cMeasure) trying to show the situation to the operator. After, the robot changed its phase to returning to base (e3) and went home. Here, is interesting to observe that the operator never looked to the areas 3, 5 or 7 (cf. Fig. 2), where the low battery event could be identified without the countermeasure.

VI. CONCLUSIONS

This study has shown the effectiveness of the MOMDP model as basis for mixed-initiative actions planning. In such cases, agents must collaborate by bringing, according to their recognized skills, the relevant elements to reach together a shared goal. In our application case, the robot counts on the operator to process a target, since the operator's interpretation is considered more reliable that the robot's. Also, the robot may needs the intervention of a human operator in cases an essential sensor for autonomous navigation breaks down. Our principal contribution in this mixed-initiative problem is that we have considered that the human operator is not a providential agent, i.e. he can be unaware of the situation. To model the problem, we have used data collected from previous experiments with an heterogeneous human-robot system. Based on it, the probability functions were assigned for the sensors failure, operator's perception about sensor's status and for the operator's cognitive availability. With the MOMDP model and a simulated environment, we checked that the collaborative system was in general able to successfully complete or terminate the mission, even when the simulated environment caused many simultaneous sensors/devices/operator failures.

Future work shall to take into account in the model more than one partially observable state variable. For the human factor community, the estimation of the operator state is obviously more complex and composed by more state variables than the one considered in this study (workload, stress, engagement, etc). In the future, we hope to take into account more cognitive states as, for instance, the operator's workload or stress level, and evaluating the policy in a real manner set-up with human operator participants.

  , a) s p(s |s, a)b(s) s p(o|s , a) s p(s |s, a)b(s)

Figure 2 .

 2 Figure 2. Operator's interface and areas of interest [24].

  p(spGpsU s = Y | auto && (gpsKO || usKO)) = = ∆T Area 6 + ∆T Area 8 ∆T all areas = 0.70 p(spGpsU s = Y | manual && (gpsKO || usKO)) = 0.86

  p(spBatt = Y | manual && battKO) = = ∆T Area 3 + ∆T Area 5 + ∆T Area 7 ∆T all areas = 0.021 p(spBatt = Y | auto && battKO) = 0.033

Figure 3 .

 3 Figure 3. Mission phases and summarized transitions. The loop transitions (red lines) indicate the transitions observed after getAttention or cMeasure actions.Table III MISSION PHASES

  Number of visits per mission phase. (b) Percentages.

http://bigbird.comp.nus.edu.sg/pmwiki/farm/appl/

http://www.pomdp.org/code/pomdp-file-spec.shtml

ACKNOWLEDGMENT

This work was supported by CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico -Brazil.