Raïssi Tarek

Efimov Denis

Based on the theory of positive systems, the goal of interval observers is to compute sets of admissible values of the state vector at each instant of time for systems subject to bounded uncertainties (noises, disturbances and parameters). The size of the estimated sets, which should be minimised, are proportional to the model uncertainties. An interval estimation can be seen as a conventional point estimation (the centre of the interval) with an estimation error given by the interval radius. The reliable uncertainties propgation performed in this context can be useful in several fields such as robust control, diagnosis and fault-tolerant control. This paper presents some recent results on interval observers for several dynamical systems classes such as continuous-time and switched systems.

Introduction

The problem of state vector estimation is very challenging and can be encountered in many applications [START_REF]Control and Observer Design for Nonlinear Finite and Infinite Dimensional Systems[END_REF][START_REF] Fossen | New Directions in Nonlinear Observer Design[END_REF][START_REF] Besançon | Nonlinear Observers and Applications[END_REF]. For linear time-invariant models there are plenty solutions, among them the most popular are Luenberger observer and Kalman filter for deterministic and stochastic settings, respectively. Among other popular solutions for estimations it is worth to mention high-gain observers [START_REF] Esfandiari | Output feedback stabilization of fully linearizable systems[END_REF] or high order sliding mode observers/differentiators [START_REF] Levant | High-order sliding modes: differentiation and output feedback control[END_REF]. In nonlinear case, observer or controller design is based on transformation of the system into a canonical form (frequently close to a linear canonical representation [START_REF]Control and Observer Design for Nonlinear Finite and Infinite Dimensional Systems[END_REF][START_REF] Fossen | New Directions in Nonlinear Observer Design[END_REF][START_REF] Besançon | Nonlinear Observers and Applications[END_REF]), since such a transformation may depend on uncertain parameters or may be unknown due to the model complexity, then application of the transformation can be an obstruction in practice. The reader can refer to a survey on application of different observers for control and estimation of nonlinear systems [START_REF] Golubev | Stabilization of nonlinear dynamic systems using the system state estimates made by the asymptotic observer[END_REF]. Il this context, the class of Linear Parameter-Varying (LPV) systems became very popular in applications: a wide class of nonlinear systems can be presented in the LPV form (in this case Raïssi Tarek, Conservatoire National des Arts et Métiers (CNAM), Cedric -laetitia 292, Rue St-Martin, 75141 Paris Efimov Denis, Inria, Univ. Lille, CNRS, UMR 9189 -CRIStAL, F-59000 Lille, France the system equations are extended). A partial linearity of LPV models allows a rich spectrum of methods developed for linear systems to be applied [START_REF] Shamma | Control of Linear Parameter Varying Systems with Applications[END_REF][START_REF] Marcos | Development of linear parameter-varying models for aircraft[END_REF][START_REF] Shamma | Gain-scheduled missile autopilot design using linear parameter-varying transformations[END_REF][START_REF] Tan | Applications of Linear Parameter-Varying Control Theory[END_REF].

In several fields such as networked control systems, electrical devices/circuits and congestion modeling [START_REF] Liberzon | Switching in systems and control[END_REF] switched systems appear as interesting modeling tools, that can be more accurate than LPV models, to take into account some complex physical behaviours. These systems are one of the most important classes of Hybrid Dynamical Systems (HDS). They consist of a family of continuous or discrete-time subsystems and a switching rule orchestrating among them. At each time only one subsystem is active. Several works treated switched systems such as [START_REF] Djemai | Hybrid Dynamical Systems[END_REF], [START_REF] Liberzon | Switching in systems and control[END_REF].

Apart of model complexity, another difficulty for an estimator design consists in the model uncertainty (unknown parameters or/and external disturbances). In the presence of uncertainty, design of a conventional estimator, converging to the ideal value of the state, cannot be realized. Application of sliding-mode tools [START_REF] Shtessel | Sliding Mode Control and Observation[END_REF] or other disturbance cancellation approaches may resolve this issue in some cases, however, in general, in the presence of uncertainty, the state estimation error is never approaching zero (it can be bounded or asymptotically bounded, and different versions of practical stability are used for analysis). In this case an interval estimation may still remain feasible: an observer can be constructed that, using input-output information, evaluates the set of admissible values (interval) for the state at each instant of time. The interval length has to be minimized by tuning the observer parameters, and it is proportional to the size of the model uncertainty. Despite such a formulation looks like a simplification of the state estimation problem (instead of the state just an interval is estimated), in fact it is an improvement since the interval mean can be used as the state pointwise estimate, while the interval width gives the admissible deviations from that value. Thus, an interval estimator provides a simultaneous accuracy evaluation for bounded uncertainty, which may not have a known statistics [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous time systems[END_REF][START_REF] Kieffer | Guaranteed nonlinear state estimator for cooperative systems[END_REF][START_REF] Olivier | Closed loop observers bundle for uncertain biotechnological models[END_REF][START_REF] Milanese | Unified set membership theory for identification, prediction and filtering of nonlinear systems[END_REF][START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF]. This idea has been initially proposed in [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF] and but recently extended to several classes of dynamical systems. This paper constitues an overview of recent results on interval observers which form a subclass of set-membership estimators and whose design is based on the monotone systems theory [START_REF] Olivier | Closed loop observers bundle for uncertain biotechnological models[END_REF][START_REF] Moisan | Near optimal interval observers bundle for uncertain bio-reactors[END_REF][START_REF] Raïssi | Interval observers design for consistency checks of nonlinear continuous-time systems[END_REF][START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF][START_REF] Efimov | Interval estimation for LPV systems applying high order sliding mode techniques[END_REF]. Only the case of continuous-times systems is considered, the reader can refer for instance to [START_REF] Efimov | On interval observer design for discrete systems[END_REF][START_REF] Efimov | Interval observers for time-varying discrete-time systems[END_REF][START_REF] Mazenc | Interval observers for discrete-time systems[END_REF][START_REF] Mazenc | Robust interval observers and stabilization design for discretetime systems with input and output[END_REF][START_REF] Briat | Interval peak-topeak observers for continuous-and discrete-time systems with persistent inputs and delays[END_REF][START_REF] Zhang | Event-triggered fault detection for a class of discrete-time linear systems using interval observers[END_REF] for the case of discrete-time systems.

In such a way the main restriction for the interval observer design consists in providing cooperativity of the interval estimation error dynamics by a proper design. Such a complexity has been overcame in [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF][START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF][START_REF] Combastel | Stable interval observers in C for linear systems with time-varying input bounds[END_REF] for Linear Time-Invariant (LTI) systems and extended to Linear Time-Varying (LTV), LPV and particular classes of nonlinear systems. In those studies, it has been shown that under some mild conditions, by applying a similarity transformation, a Hurwitz matrix could be transformed to a Hurwitz and Metzler one. In the following the main ideas of interval observer design are explained for different classes of systems.

Preliminaries

Real and integer numbers sets are denoted by R and Z respectively, R + = {𝜏 ∈ R : 𝜏 ≥ 0} and Z + = Z∩R + . Euclidean norm for a vector 𝑥 ∈ R 𝑛 will be denoted as |𝑥|, and for a measurable and locally essentially bounded input 𝑢 :

R + → R (R + = {𝜏 ∈ R : 𝜏 ≥ 0}) the symbol ||𝑢|| [𝑡0,𝑡1] de- notes its 𝐿 ∞ norm ||𝑢|| [𝑡0,𝑡1] = 𝑒𝑠𝑠 sup{|𝑢(𝑡)|, 𝑡 ∈ [𝑡 0 , 𝑡 1 ]}, if 𝑡 1 =
+∞ then we will simply write ||𝑢||. We will denote as ℒ ∞ the set of all inputs 𝑢 with the property ||𝑢|| < ∞. Denote the sequence of integers 1, ..., 𝑘 as 1, 𝑘. The symbols 𝐼 𝑛 , 𝐸 𝑛×𝑚 and 𝐸 𝑝 denote the identity matrix with dimension 𝑛 × 𝑛, the matrix with all elements equal 1 with dimensions 𝑛 × 𝑚 and 𝑝 × 1, respectively. For a matrix 𝐴 ∈ R 𝑛×𝑛 the vector of its eigenvalues is denoted as 𝜆(𝐴), ||𝐴|| 𝑚𝑎𝑥 = max 𝑖=1,𝑛,𝑗=1,𝑛 |𝐴 𝑖,𝑗 | (the elementwise maximum norm, it is not sub-multiplicative) and ||𝐴|| 2 = √ max 𝑖=1,𝑛 𝜆 𝑖 (𝐴 T 𝐴) (the induced 𝐿 2 matrix norm).

Interval relations

𝐴 + 𝑥 -𝐴 -𝑥 ≤ 𝐴𝑥 ≤ 𝐴 + 𝑥 -𝐴 -𝑥. ( 1
)
(2) If 𝐴 ∈ R 𝑚×𝑛 is a matrix variable and 𝐴 ≤ 𝐴 ≤ 𝐴 for some 𝐴, 𝐴 ∈ R 𝑚×𝑛 , then

𝐴 + 𝑥 + -𝐴 + 𝑥 --𝐴 -𝑥 + + 𝐴 -𝑥 -≤ 𝐴𝑥 (2) ≤ 𝐴 + 𝑥 + -𝐴 + 𝑥 --𝐴 -𝑥 + + 𝐴 -𝑥 -.
Furthermore, if -𝐴 = 𝐴 ≤ 0 ≤ 𝐴, then the inequality (2) can be simplified: -𝐴(𝑥 + + 𝑥 -) ≤ 𝐴𝑥 ≤ 𝐴(𝑥 + + 𝑥 -). [START_REF] Briat | Robust stability analysis of uncertain linear positive systems via integral linear constraints: 𝑙 1and 𝑙 ∞ -gain characterizations[END_REF][START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF], for this kind of systems these gains are interrelated.

Nonnegative continuous-time linear systems

A matrix 𝐴 ∈ R 𝑛×𝑛 is
Lemma 2. [START_REF] Briat | Robust stability analysis of uncertain linear positive systems via integral linear constraints: 𝑙 1and 𝑙 ∞ -gain characterizations[END_REF][START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF] Let the system (3) be nonnegative (i.e. 𝐴 is Metzler, 𝐵 ≥ 0, 𝐶 ≥ 0 and 𝐷 ≥ 0), then it is asymptotically stable if and only if there exist 𝜆 ∈ R 𝑛 + ∖{0} and a scalar 𝛾 > 0 such that the following LP problem is feasible:

[ 𝐴 T 𝜆 + 𝐶 T 𝐸 𝑝 𝐵 T 𝜆 -𝛾𝐸 𝑞 + 𝐷 T 𝐸 𝑝 ] < 0.
Moreover, in this case 𝐿 1 gain of the operator 𝜔 → 𝑦 is lower than 𝛾.

Lemma 3. [START_REF] Briat | Robust stability analysis of uncertain linear positive systems via integral linear constraints: 𝑙 1and 𝑙 ∞ -gain characterizations[END_REF][START_REF] Ebihara | L1 gain analysis of linear positive systems and its application[END_REF] Let the system (3) be cooperative(i.e. 𝐴 is Metzler, 𝐵 ≥ 0, 𝐶 ≥ 0 and 𝐷 ≥ 0), then it is asymptotically stable if and only if there exist 𝜆 ∈ R 𝑛 + ∖{0} and a scalar 𝛾 > 0 such that the following LP problem is feasible:

[ 𝐴𝜆 + 𝐵𝐸 𝑞 𝐶𝜆 -𝛾𝐸 𝑝 + 𝐷𝐸 𝑞 ] < 0.
Moreover, in this case 𝐿 ∞ gain of the transfer 𝜔 → 𝑦 is lower than 𝛾.

The conventional results and definitions on 𝐿 2 /𝐿 ∞ stability for linear systems can be found in [START_REF] Khalil | Nonlinear Systems[END_REF].

Nonnegative systems with delays

Consider a linear system with constant delays

ẋ(𝑡) = 𝐴 0 𝑥(𝑡) + 𝑁 ∑ 𝑖=1 𝐴 𝑖 𝑥(𝑡 -𝜏 𝑖 ) + 𝜔(𝑡), (4)
where

𝑥(𝑡) ∈ R 𝑛 is the state, 𝑥 𝑡 ∈ 𝒞 𝑛 𝜏 for 𝜏 = max 1≤𝑖≤𝑁 𝜏 𝑖 where 𝜏 𝑖 ∈ R + are the delays (𝒞 𝜏 = 𝐶([-𝜏, 0], R) is the set of continuous maps from [-𝜏, 0] into R, 𝒞 𝜏 + = {𝑦 ∈ 𝒞 𝜏 : 𝑦(𝑠) ∈ R + , 𝑠 ∈ [-𝜏, 0]}); a piecewise continuous function 𝜔 ∈ ℒ 𝑛
∞ is the input; the constant matrices 𝐴 𝑖 , 𝑖 = 0, 𝑁 have appropriate dimensions. The system (4) is called cooperative or nonnegative [START_REF] Haddad | Stability theory for nonnegative and compartmental dynamical systems with time delay[END_REF] if it admits 𝑥(𝑡) ∈ R 𝑛 + for all 𝑡 ≥ 𝑡 0 provided that 𝑥 𝑡0 ∈ 𝒞 𝑛 𝜏 + and 𝜔 : R → R 𝑛 + .

Lemma 4. [START_REF] Dambrine | Stability analysis of time-delay systems[END_REF][START_REF] Haddad | Stability theory for nonnegative and compartmental dynamical systems with time delay[END_REF] The system ( 4) is nonnegative for all 𝜏 ∈ R + iff 𝐴 0 is Metzler and 𝐴 𝑖 , 𝑖 = 1, 𝑁 are nonnegative matrices.

Interval observers for LTI systems

Let us start the line of design of interval observers with the simplest cases of time-invariant linear models. Consider the following system

ẋ(𝑡) = 𝐴𝑥(𝑡) + 𝑑(𝑡), 𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝑣(𝑡), 𝑡 ∈ R + (5)
where

𝑥(𝑡) ∈ R 𝑛 is the state, 𝑦(𝑡) ∈ R 𝑝 is the output; 𝑑(𝑡) ∈ R 𝑛 is the disturbance, 𝑑 ∈ ℒ 𝑛 ∞ ; 𝑣(𝑡) ∈ R 𝑝 is the measurement noise, 𝑣 ∈ ℒ 𝑝
∞ ; the matrices 𝐴, 𝐶 have appropriate dimensions. This model has three sources of uncertainty: initial conditions for 𝑥(0), instant values of 𝑑 and 𝑣. It is assumed that all these uncertain factors belong to known intervals.

Assumption 1. Let 𝑥(0) ∈ [𝑥 0 , 𝑥 0 ] for some known 𝑥 0 , 𝑥 0 ∈ R 𝑛 , let also two functions 𝑑, 𝑑 ∈ ℒ 𝑛
∞ and a constant 𝑉 > 0 be given such that

𝑑(𝑡) ≤ 𝑑(𝑡) ≤ 𝑑(𝑡), |𝑣(𝑡)| ≤ 𝑉 ∀𝑡 ≥ 0.
Using the available information, the goal is to calculate two bounds 𝑥, 𝑥 ∈ ℒ 𝑛 ∞ verifying

𝑥(𝑡) ≤ 𝑥(𝑡) ≤ 𝑥(𝑡) ∀𝑡 ≥ 0. ( 6
)
An interval observer, composed of two conventional ones, is a solution to this problem:

⎧ ⎨ ⎩ ẋ(𝑡) = 𝐴𝑥(𝑡) + 𝐿[𝑦(𝑡) -𝐶𝑥(𝑡)] -|𝐿|𝐸 𝑝 𝑉 + 𝑑(𝑡), ẋ(𝑡) = 𝐴𝑥(𝑡) + 𝐿[𝑦(𝑡) -𝐶𝑥(𝑡)] + |𝐿|𝐸 𝑝 𝑉 + 𝑑(𝑡), 𝑥(0) = 𝑥 0 , 𝑥(0) = 𝑥 0 , (7)
where 𝐿 ∈ R 𝑛×𝑝 is the observer gain to be designed. The conditions to satisfy for 𝐿 are given below.

Theorem 1. [19]

Let Assumption 1 hold, the matrix 𝐴 -𝐿𝐶 be Metzler and 𝑥 ∈ ℒ 𝑛 ∞ , then the solutions of (5) and (7) satisfy the relations [START_REF] Golubev | Stabilization of nonlinear dynamic systems using the system state estimates made by the asymptotic observer[END_REF]. In addition, 𝑥, 𝑥 ∈ ℒ 𝑛 ∞ if 𝐴 -𝐿𝐶 is Hurwitz.

In the following only the proof of this theorem is given as an illustration.

Proof. Define two estimation errors

𝑒(𝑡) = 𝑥(𝑡) -𝑥(𝑡), 𝑒(𝑡) = 𝑥(𝑡) -𝑥(𝑡), which yield differential equations: ė(𝑡) = [𝐴 -𝐿𝐶]𝑒(𝑡) + 𝐿𝑣(𝑡) + |𝐿|𝐸 𝑝 𝑉 + 𝑑(𝑡) -𝑑(𝑡), ė(𝑡) = [𝐴 -𝐿𝐶]𝑒(𝑡) -𝐿𝑣(𝑡) + |𝐿|𝐸 𝑝 𝑉 + 𝑑(𝑡) -𝑑(𝑡).
By Assumption 1, for all 𝑡 ≥ 0

|𝐿|𝐸 𝑝 𝑉 ± 𝐿𝑣(𝑡) ≥ 0, 𝑑(𝑡) -𝑑(𝑡) ≥ 0, 𝑑(𝑡) -𝑑(𝑡) ≥ 0.
If 𝐴 -𝐿𝐶 is a Metzler matrix, since all inputs of 𝑒(𝑡), 𝑒(𝑡) are positive and 𝑒(0) ≥ 0, 𝑒(𝑡) ≥ 0, then 𝑒(𝑡) ≥ 0, 𝑒(𝑡) ≥ 0 for all 𝑡 ≥ 0 [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF][START_REF] Smith | Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems[END_REF]. The property [START_REF] Golubev | Stabilization of nonlinear dynamic systems using the system state estimates made by the asymptotic observer[END_REF] follows from these relations. If 𝐴-𝐿𝐶 is Hurwitz, since all inputs of 𝑒(𝑡), 𝑒(𝑡) are bounded, then 𝑒, 𝑒 ∈ ℒ 𝑛 ∞ and boundedness of 𝑥, 𝑥 is implied by boundedness of 𝑥.

From this proof we can conclude that the idea to design an interval observer is to guarantee nonnegativity of the estimation error dynamics. The observer gain 𝐿 has to be designed such that the matrix 𝐴 -𝐿𝐶 is Metzler and Hurwitz. In addition, in order to optimize the width of estimated interval [𝑥(𝑡), 𝑥(𝑡)] the problem of 𝐿 1 or 𝐿 ∞ optimization of that gain value can be posed. Using lemmas 2 and 3 this problem can be formulated as a LP computational procedure. The reader can also refer to [START_REF] Mccarthy | Cardinality constrained robust optimization applied to a class of interval observers[END_REF] for another LP problem formulation and solution for a more generic LPV system, and to [START_REF] Ellero | Unknown input interval observer with ℎ∞ and dstability performance[END_REF] for a 𝐻 ∞ oprimization procedure. In the case of a 𝐿 1 optimization, it is necessary to find

𝜆 ∈ R 𝑛 , 𝑤 ∈ R 𝑝 and a diagonal matrix 𝑀 ∈ R 𝑛×𝑛 such that [ 𝐴 T 𝜆 -𝐶 T 𝑤 + 𝐸 𝑛 𝜆 -𝛾𝐸 𝑛 ] < 0, 𝜆 > 0, 𝑀 ≥ 0, ( 8
)
𝐴 T 𝜆 -𝐶 T 𝑤 + 𝑀 𝜆 ≥ 0,
then 𝑤 = 𝐿 T 𝜆 [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF]. Formulation [START_REF] Marcos | Development of linear parameter-varying models for aircraft[END_REF] provides an effective computational tool to design interval observers, but it is only sufficient condition and in some cases this LP problem may have no solution, but it does not imply that it is not possible to design an interval observer. Roughly speaking in this case it is not possible to find 𝐿 such that 𝐴 -𝐿𝐶 is simultaneously Metzler and Hurwitz. However, it is well known that Hurwitz property of a matrix is preserved under a similarity transformation of coordinates, then to overcome the issue it is possible to design the gain 𝐿 such that the matrix 𝐴 -𝐿𝐶 is Hurwitz and next to find a nonsingular matrix 𝑆 ∈ R 𝑛×𝑛 such that in the new coordinates 𝑧 = 𝑆𝑥 the state matrix 𝐷 = 𝑆(𝐴 -𝐿𝐶)𝑆 -1 is Metzler (it is Hurwitz by construction). The conditions of existence of such a real transformation matrix 𝑆 are given in the following lemma. Lemma 5. [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF] Given the matrices 𝐴 ∈ R 𝑛×𝑛 , 𝐷 ∈ R 𝑛×𝑛 and 𝐶 ∈ R 𝑝×𝑛 . If there is a matrix 𝐿 ∈ R 𝑛×𝑝 such that the matrices 𝐴 -𝐿𝐶 and 𝐷 have the same eigenvalues, then there is a matrix 𝑆 ∈ R 𝑛×𝑛 such that 𝐷 = 𝑆(𝐴 -𝐿𝐶)𝑆 -1 provided that the pairs (𝐴 -𝐿𝐶, 𝜒 1 ) and (𝐷, 𝜒 2 ) are observable for some

𝜒 1 ∈ R 1×𝑛 , 𝜒 2 ∈ R 1×𝑛 .
Note that (5) can be rewritten as follows:

ẋ(𝑡) = (𝐴 -𝐿𝐶)𝑥(𝑡) + 𝐿𝑦(𝑡) -𝐿𝑣(𝑡) + 𝑑(𝑡).
Under conditions of this lemma and in the new coordinates 𝑧 = 𝑆𝑥 the system (5) takes the form:

ż(𝑡) = 𝐷𝑧(𝑡) + 𝑆𝐿𝑦(𝑡) + 𝛿(𝑡), 𝛿(𝑡) = 𝑆[𝑑(𝑡) -𝐿𝑣(𝑡)]. (9)
And using Lemma 1 we obtain that 𝛿(𝑡) ≤ 𝛿(𝑡) ≤ 𝛿(𝑡), where 𝛿(𝑡) = 𝑆 + 𝑑(𝑡) -𝑆 -𝑑(𝑡) -|𝑆𝐿|𝐸 𝑝 𝑉 and 𝛿(𝑡) = 𝑆 + 𝑑(𝑡) -𝑆 -𝑑(𝑡) + |𝑆𝐿|𝐸 𝑝 𝑉 . For the system (9) all conditions of Theorem 1 are satisfied and an interval observer similar to (7) can be designed:

ż(𝑡) = 𝐷𝑧(𝑡) + 𝑆𝐿𝑦(𝑡) + 𝛿(𝑡), ż(𝑡) = 𝐷𝑧(𝑡) + 𝑆𝐿𝑦(𝑡) + 𝛿(𝑡), 𝑧(0) = 𝑆 + 𝑥 0 -𝑆 -𝑥 0 , 𝑧(0) = 𝑆 + 𝑥 0 -𝑆 -𝑥 0 , 𝑥(𝑡) = (𝑆 -1 ) + 𝑧(𝑡) -(𝑆 -1 ) -𝑧(𝑡), ( 10
)
𝑥(𝑡) = (𝑆 -1 ) + 𝑧(𝑡) -(𝑆 -1 ) -𝑧(𝑡),
where the relations (1) are used to calculate the initial conditions for 𝑧, 𝑧 and the estimates 𝑥, 𝑥. It is easy to show that in [START_REF] Tan | Applications of Linear Parameter-Varying Control Theory[END_REF] the inclusion ( 6) is satisfied and 𝑥, 𝑥 ∈ ℒ 𝑛 ∞ . If conditions of Lemma 5 are not satisfied, then it is possible also to find a time-varying change of coordinates [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF][START_REF] Combastel | Stable interval observers in C for linear systems with time-varying input bounds[END_REF], for instance, based on the Jordan canonical form, as in the next lemma. Lemma 6. [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF] Let 𝐴 -𝐿𝐶 be Hurwitz, then there exists an invertible matrix function

𝑃 : R → R 𝑛×𝑛 , of class 𝐶 ∞ elementwise, ||𝑃 (𝑡)|| 2 < +∞ for all 𝑡 ∈ R, such that for all 𝑡 ∈ R Ṗ (𝑡) = 𝐷𝑃 (𝑡) -𝑃 (𝑡)(𝐴 -𝐿𝐶),
where 𝐷 ∈ R 𝑛×𝑛 is a Hurwitz and Metzler matrix.

Interval observers for time-delay systems

Consider the system (4) with an output 𝑦 ∈ R 𝑝 subject to a bounded noise 𝑣 ∈ ℒ 𝑝 ∞ :

𝑦 = 𝐶𝑥, 𝜓 = 𝑦 + 𝑣(𝑡), ( 11
)
where 𝐶 ∈ R 𝑝×𝑛 . Below the relation 𝑎 ≤ 𝑏 for 𝑎, 𝑏 ∈ 𝒞 𝑛 𝜏 is understood in the sense that 𝑎(𝜃) ≤ 𝑏(𝜃) for all 𝜃 ∈ [-𝜏, 0].

Assumption 2. Let 𝑥 ∈ ℒ 𝑛

∞ with 𝑥 0 ≤ 𝑥 𝑡0 ≤ 𝑥 0 for some 𝑥 0 , 𝑥 0 ∈ 𝒞 𝑛 𝜏 ; ||𝑣|| ≤ 𝑉 for a given 𝑉 > 0; and 𝜔(𝑡) ≤ 𝜔(𝑡) ≤ 𝜔(𝑡) for all 𝑡 ≥ 𝑡 0 for some known 𝜔, 𝜔 ∈ ℒ 𝑛 ∞ .

In this assumption it is supposed that the state of the system (4) is bounded with an unknown upper bound, but with a specified admissible set for initial conditions [𝑥 0 , 𝑥 0 ]. Uncertainty of the system is collected in the external input 𝜔 with known bounds on the incertitude 𝜔, 𝜔, the measurement noise 𝑣 and the interval of initial conditions [𝑥 0 , 𝑥 0 ].

As we see above, interval observers have an enlarged dimension (the examples given above have 2𝑛 variables to estimate 𝑛 states). Thus, design of reduced order interval observers is of a big importance for applications. A reduced order interval observer for time-delay system (4), [START_REF] Liberzon | Switching in systems and control[END_REF] has been proposed in [START_REF] Efimov | Interval estimation for uncertain systems with timevarying delays[END_REF][START_REF] Efimov | On reduced-order interval observers for time-delay systems[END_REF], those ideas are explained below.

For the system (4), [START_REF] Liberzon | Switching in systems and control[END_REF] there exists a nonsingular matrix 𝑆 ∈ R 𝑛×𝑛 such that 𝑥 = 𝑆 [𝑦 T 𝑧 T ] T for an auxiliary variable 𝑧 ∈ R 𝑛-𝑝 (define

𝑆 -1 = [𝐶 T 𝑍 T ] T for a matrix 𝑍 ∈ R (𝑛-𝑝)×𝑛 ), then ẏ(𝑡) = 𝑅 1 𝑦(𝑡) + 𝑅 2 𝑧(𝑡) + 𝑁 ∑ 𝑖=1 [𝐷 1𝑖 𝑦(𝑡 -𝜏 𝑖 ) + 𝐷 2𝑖 𝑧(𝑡 -𝜏 𝑖 )] + 𝐶𝜔(𝑡), (12) ż(𝑡) = 𝑅 3 𝑦(𝑡) + 𝑅 4 𝑧(𝑡) + 𝑁 ∑ 𝑖=1 [𝐷 3𝑖 𝑦(𝑡 -𝜏 𝑖 ) + 𝐷 4𝑖 𝑧(𝑡 -𝜏 𝑖 )] + 𝑍𝜔(𝑡),
for some matrices 𝑅 𝑘 , 𝐷 𝑘𝑖 , 𝑘 = 1, 4, 𝑖 = 1, 𝑁 of appropriate dimensions. Introducing a new variable 𝑤 = 𝑧 -𝐾𝑦 = 𝑈 𝑥 for a matrix 𝐾 ∈ R (𝑛-𝑝)×𝑝 with 𝑈 = 𝑍 -𝐾𝐶, from (12) the following equation is obtained

ẇ(𝑡) = 𝐺 0 𝜓(𝑡) + 𝑀 0 𝑤(𝑡) + 𝑁 ∑ 𝑖=1 [𝐺 𝑖 𝜓(𝑡 -𝜏 𝑖 ) + 𝑀 𝑖 𝑤(𝑡 -𝜏 𝑖 )] + 𝛽(𝑡), 𝛽(𝑡) = 𝑈 𝜔(𝑡) -𝐺 0 𝑣(𝑡) - 𝑁 ∑ 𝑖=1
𝐺 𝑖 𝑣(𝑡 -𝜏 𝑖 ), [START_REF] Shtessel | Sliding Mode Control and Observation[END_REF] where 𝜓(𝑡) is defined in [START_REF] Liberzon | Switching in systems and control[END_REF], 𝐺 0 = 𝑅 3 -𝐾𝑅 1 + (𝑅 4 -𝐾𝑅 2 )𝐾, 𝑀 0 = 𝑅 4 -𝐾𝑅 2 , and 𝐺 𝑖 = 𝐷 3𝑖 -𝐾𝐷 1𝑖 + {𝐷 4𝑖 -𝐾𝐷 2𝑖 }𝐾, 𝑀 𝑖 = 𝐷 4𝑖 -𝐾𝐷 2𝑖 for 𝑖 = 1, 𝑁 . Under Assumption 2 and using the relations (1) the following inequalities follow:

𝛽(𝑡) ≤ 𝛽(𝑡) ≤ 𝛽(𝑡), 𝛽(𝑡) = 𝑈 + 𝜔(𝑡) -𝑈 -𝜔(𝑡) - 𝑁 ∑ 𝑖=0 |𝐺 𝑖 |𝐸 𝑝 𝑉, 𝛽(𝑡) = 𝑈 + 𝜔(𝑡) -𝑈 -𝜔(𝑡) + 𝑁 ∑ 𝑖=0 |𝐺 𝑖 |𝐸 𝑝 𝑉.
Then the next interval reduced-order observer can be proposed for (4):

ẇ(𝑡) = 𝐺 0 𝜓(𝑡) + 𝑀 0 𝑤(𝑡) + 𝑁 ∑ 𝑖=1 [𝐺 𝑖 𝜓(𝑡 -𝜏 𝑖 ) + 𝑀 𝑖 𝑤(𝑡 -𝜏 𝑖 )] + 𝛽(𝑡), (14) ẇ(𝑡) = 𝐺 0 𝜓(𝑡) + 𝑀 0 𝑤(𝑡) + 𝑁 ∑ 𝑖=1 [𝐺 𝑖 𝜓(𝑡 -𝜏 𝑖 ) + 𝑀 𝑖 𝑤(𝑡 -𝜏 𝑖 )] + 𝛽(𝑡), 𝑤 0 = 𝑈 + 𝑥 0 -𝑈 -𝑥 0 , 𝑤 0 = 𝑈 + 𝑥 0 -𝑈 -𝑥 0 .
The applicability conditions for ( 14) are given below. Theorem 2. [START_REF] Efimov | Interval estimation for uncertain systems with timevarying delays[END_REF] Let Assumption 2 be satisfied and the matrices 𝑀 0 , 𝑀 𝑖 , 𝑖 = 1, 𝑁 form an asymptotically stable cooperative system. Then 𝑥, 𝑥 ∈ ℒ 𝑛 ∞ and

𝑥(𝑡) ≤ 𝑥(𝑡) ≤ 𝑥(𝑡) ∀𝑡 ≥ 0,
where

𝑥(𝑡) = 𝑆 + [𝑦(𝑡) T 𝑧(𝑡) T ] T -𝑆 -[𝑦(𝑡) T 𝑧(𝑡) T ] T , ( 15
)
𝑥(𝑡) = 𝑆 + [𝑦(𝑡) T 𝑧(𝑡) T ] T -𝑆 -[𝑦(𝑡) T 𝑧(𝑡) T ] T , 𝑦(𝑡) = 𝜓(𝑡) -𝑉, 𝑦(𝑡) = 𝜓(𝑡) + 𝑉, ( 16
)
𝑧(𝑡) = 𝑤(𝑡) + 𝐾 + 𝑦 -𝐾 -𝑦, 𝑧(𝑡) = 𝑤(𝑡) + 𝐾 + 𝑦 -𝐾 -𝑦.
The main condition of Theorem 2 is rather straightforward: the matrices 𝑀 0 , 𝑀 𝑖 , 𝑖 = 1, 𝑁 have to form a stable cooperative system. It is a standard LMI problem to find a matrix 𝐾 such that the system composed by 𝑀 0 , 𝑀 𝑖 , 𝑖 = 1, 𝑁 is stable, but to find a matrix 𝐾 making the system stable and cooperative simultaneously could be more complicated. However, the advantage of Theorem 2 is that its main condition can be reformulated using LMIs following the idea of [START_REF] Rami | Positive observation problem for linear time-delay positive systems[END_REF].

Proposition 1. [41] Let there exist 𝜍 ∈ R + , 𝑝 ∈ R 𝑛-𝑝 + , 𝑞 ∈ R 𝑛-𝑝 +
and 𝐵 ∈ R (𝑛-𝑝)×𝑝 such that the following LMIs are satisfied: [START_REF] Jaulin | Nonlinear bounded-error state estimation of continuous time systems[END_REF].

𝑝 T Π 0 -𝐸 T 𝑛-𝑝 𝐵Π 1 + 𝑞 T ≤ 0, 𝑝 > 0, 𝑞 > 0; diag[𝑝]𝑅 4 -𝐵𝑅 2 + 𝜍𝐼 𝑛-𝑝 ≥ 0, 𝜍 > 0; diag[𝑝]𝐷 4𝑖 -𝐵𝐷 2𝑖 ≥ 0, 𝑖 = 1, 𝑁 ; Π 0 = 𝑅 4 + 𝑁 ∑ 𝑖=1 𝐷 4𝑖 , Π 1 = 𝑅 2 + 𝑁 ∑ 𝑖=1 𝐷 2𝑖 , then 𝐾 = diag[𝑝] -1 𝐵 and the matrices 𝑀 0 = 𝑅 4 -𝐾𝑅 2 , 𝑀 𝑖 = 𝐷 4𝑖 -𝐾𝐷 2𝑖 , 𝑖 = 1, 𝑁 represent a stable cooperative system in
If these LMIs are not satisfied, the assumption that the matrix 𝑀 0 is Metzler and the matrices 𝑀 𝑖 , 𝑖 = 1, 𝑁 are nonnegative can be relaxed using Lemma 5[41]. This approach can also be applied to the system (4) with timevarying and uncertain delays [START_REF] Efimov | Interval estimation for uncertain systems with timevarying delays[END_REF]. The cases of delayed measurements are studied in [START_REF] Mazenc | Exponentially stable interval observers for linear systems with delay[END_REF][START_REF] Efimov | On reduced-order interval observers for time-delay systems[END_REF]. An extension to descriptor delay systems is presented in [START_REF] Efimov | Interval observer design for estimation and control of timedelay descriptor systems[END_REF].

The case of LPV models

Consider a LPV system:

ẋ = [𝐴 0 + Δ𝐴(𝜌(𝑡))]𝑥 + 𝑑(𝑡), 𝑦 = 𝐶𝑥 + 𝑣(𝑡), ( 17
)
where 𝑥 ∈ R 𝑛 is the state, 𝑦 ∈ R 𝑝 is the output available for measurements, 𝜌(𝑡) ∈ Π ⊂ R 𝑟 is the vector of scheduling parameters, 𝜌 ∈ ℒ 𝑟 ∞ . The values of the scheduling vector 𝜌 are not available for measurements, and only the set of admissible values Π is known. The matrices 𝐴 0 ∈ R 𝑛×𝑛 and 𝐶 ∈ R 𝑝×𝑛 are known, the matrix function Δ𝐴 : Π → R 𝑛×𝑛 is piecewise continuous and also known for a given value of 𝜌. The signals 𝑑 : R + → R 𝑛 and 𝑣 : R + → R 𝑝 are the external input and measurement noise respectively.

Assumption 3. 𝑥 ∈ ℒ 𝑛

∞ and 𝑥(0) ∈ [𝑥 0 , 𝑥 0 ] for some known 𝑥 0 , 𝑥 0 ∈ R 𝑛 ; 𝑑(𝑡) ≤ 𝑑(𝑡) ≤ 𝑑(𝑡) and |𝑣(𝑡)| ≤ 𝑉 for all 𝑡 ≥ 0 and some known 𝑑, 𝑑 ∈ ℒ 𝑛 ∞ and 𝑉 > 0; Δ𝐴 ≤ Δ𝐴(𝜌) ≤ Δ𝐴 for all 𝜌 ∈ Π and some known Δ𝐴, Δ𝐴 ∈ R 𝑛×𝑛 .

In the system [START_REF] Milanese | Unified set membership theory for identification, prediction and filtering of nonlinear systems[END_REF], the uncertainties are due to the initial state, the noises, the disturbances and also to the fact that the matrix Δ𝐴(𝜌) belongs into the interval [Δ𝐴, Δ𝐴] for all 𝑡 ≥ 0. The interval [Δ𝐴, Δ𝐴] is easy to compute for a given set Π and known Δ𝐴(𝜌).

The objective is to design an interval observer for the system [START_REF] Milanese | Unified set membership theory for identification, prediction and filtering of nonlinear systems[END_REF]. Before introduction of interval observer equations note that for a matrix 𝐿 ∈ R 𝑛×𝑝 the system (17) can be rewritten as follows:

ẋ = [𝐴 0 -𝐿𝐶]𝑥 + Δ𝐴(𝜌(𝑡))𝑥 + 𝐿[𝑦 -𝑣(𝑡)] + 𝑑(𝑡),
and according to Lemma 1 and Assumption 3 we have for all 𝜌 ∈ Π:

Δ𝐴 + 𝑥 + -Δ𝐴 + 𝑥 --Δ𝐴 -𝑥 + + Δ𝐴 -𝑥 -≤ Δ𝐴(𝜌)𝑥 (18)
≤ Δ𝐴 + 𝑥 + -Δ𝐴 + 𝑥 --Δ𝐴 -𝑥 + + Δ𝐴 -𝑥 - provided that 𝑥 ≤ 𝑥 ≤ 𝑥 for some 𝑥, 𝑥 ∈ R 𝑛 .
The following interval observer structure is proposed [START_REF] Rami | Tight robust interval observers: an LP approach[END_REF][START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF] for the LPV system [START_REF] Milanese | Unified set membership theory for identification, prediction and filtering of nonlinear systems[END_REF]:

ẋ = [𝐴 0 -𝐿𝐶]𝑥 + [Δ𝐴 + 𝑥 + -Δ𝐴 + 𝑥 - -Δ𝐴 -𝑥 + + Δ𝐴 -𝑥 -] + 𝐿𝑦 -|𝐿|𝑉 𝐸 𝑝 + 𝑑(𝑡), ẋ = [𝐴 0 -𝐿𝐶]𝑥 + [Δ𝐴 + 𝑥 + -Δ𝐴 + 𝑥 - (19)
-Δ𝐴

-𝑥 + + Δ𝐴 -𝑥 -] + 𝐿𝑦 + |𝐿|𝑉 𝐸 𝑝 + 𝑑(𝑡), 𝑥(0) = 𝑥 0 , 𝑥(0) = 𝑥 0 .
Note that due to the presence of 𝑥 + , 𝑥 -, 𝑥 + and 𝑥 -, the interval observer ( 19) is a globally Lipschitz nonlinear system. In addition, in [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF] the dynamics of 𝑥 and 𝑥 are coupled.

Theorem 3. [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF] Let assumption 3 be satisfied and the matrices 𝐴 0 -𝐿𝐶, 𝐴 0 -𝐿𝐶 be Metzler. Then the relations 𝑥(𝑡) ≤ 𝑥(𝑡) ≤ 𝑥(𝑡) hold for all 𝑡 ≥ 𝑡 0 . Furthermore, if there exist 𝑃 ∈ R 2𝑛×2𝑛 , 𝑃 = 𝑃 T ≻ 0 and 𝛾 > 0 such that the following Riccati matrix inequality is verified

𝐺 T 𝑃 + 𝑃 𝐺 + 2𝛾 -2 𝑃 2 + 𝛾 2 𝜂 2 𝐼 2𝑛 + 𝑍 T 𝑍 ≺ 0, where 𝜂 = 2𝑛||Δ𝐴 -Δ𝐴|| 𝑚𝑎𝑥 , 𝑍 ∈ R 𝑠×2𝑛 , 0 < 𝑠 ≤ 2𝑛 and 𝐺 = [ 𝐴 0 -𝐿𝐶 + Δ𝐴 + -Δ𝐴 - -Δ𝐴 - 𝐴 0 -𝐿𝐶 + Δ𝐴 + ] , then 𝑥, 𝑥 ∈ ℒ 𝑛 ∞ .
In addition, the operator

[ 𝑏 𝑏 ] → 𝑍 [ 𝑥 𝑥
] in [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF] has an 𝐿 2 gain less than 𝛾.

Note that if there are no gains such that 𝐴 0 -𝐿𝐶 and 𝐴 0 -𝐿𝐶 are Metzler, then a transformation of coordinates 𝑇 can be used for [START_REF] Gouzé | Interval observers for uncertain biological systems[END_REF] in order to relax such requirement. The Riccati matrix inequality from Theorem 3 can be reformulated in terms of LMIs with respect to 𝐿, 𝐿 and 𝑃 [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF].

Example 1. Consider a nonlinear system [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF]:

ẋ = ⎡ ⎣ 𝜖 cos 𝑡 1 + 𝜖 sin 𝑥 3 𝜀 sin 𝑥 2 𝜀 sin 𝑥 3 -0.5 + 𝜖 sin 𝑡 1 + 𝜖 cos 2𝑡 𝜀 sin 𝑥 2 0.3 + 𝜀 cos 2𝑡 -1 + 𝜖 sin 𝑡 ⎤ ⎦ 𝑥 + ⎡ ⎣ 6 cos 𝑥 1 sin 𝑡 + 0.1 sin 𝑥 3 -cos 3𝑡 + 0.1 sin 2𝑥 2 ⎤ ⎦ , 𝑦 = 𝑥 1 + 𝑣(𝑡),
where 𝜖 = 0.01 and 𝜀 = 0.001. We assume that 𝑉 = 0.1, and for simulation we selected 𝑣(𝑡) = 𝑉 (sin 5𝑡 + cos 3𝑡)/2. For initial conditions |𝑥 𝑖 (0)| ≤ 5 the system has bounded The results of simulations for an LPV system [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF] solutions. This system can be presented in the form of [START_REF] Milanese | Unified set membership theory for identification, prediction and filtering of nonlinear systems[END_REF] for

𝐴 0 = ⎡ ⎣ 0 1 0 0 -0.5 1 0 0.3 -1 ⎤ ⎦ , Δ𝐴 = ⎡ ⎣ 𝜖 𝜖 𝜀 𝜀 𝜖 𝜖 𝜀 𝜀 𝜖 ⎤ ⎦ = -Δ𝐴, 𝑏(𝑡, 𝑦) = ⎡ ⎣ 6𝑓 (𝑦) sin 𝑡 -0.1 -cos 3𝑡 -0.1 ⎤ ⎦ , 𝑏(𝑡, 𝑦) = ⎡ ⎣ 6𝑓 (𝑦) sin 𝑡 + 0.1 -cos 3𝑡 + 0.1 ⎤ ⎦ , 𝑓 (𝑦) = { cos 𝑦 cos 𝑉 if cos 𝑦 ≥ 0 cos 𝑦 if cos 𝑦 < 0 -| sin 𝑦| sin 𝑉, 𝑓 (𝑦) = { cos 𝑦 if cos 𝑦 ≥ 0 cos 𝑦 cos 𝑉 if cos 𝑦 < 0 + | sin 𝑦| sin 𝑉
and a properly selected 𝜌, clearly assumption 3 is satisfied. Solution of LMIs formulated for Theorem 3 gives the 𝐿 2 optimal solution [START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF]:

𝐿 = ⎡ ⎣ 82.923 -3𝑒 -4 -4𝑒 -4 ⎤ ⎦ , 𝐿 = ⎡ ⎣ 97.16 -2𝑒 -5 -1𝑒 -5 ⎤ ⎦
for 𝛾 = 31.4. The results of interval simulations for the variables 𝑥 2 and 𝑥 3 are given in Fig. 1 (the variable 𝑥 1 is omitted since it is available for measurements).

Case of LTV systems

LTV models can be considered as a special case of LPV ones when the vector of scheduling parameters is available for measurements. Interval observers have been proposed for LTV systems in [START_REF] Efimov | Interval state observer for nonlinear time varying systems[END_REF][START_REF] Thabet | An effective method to interval observer design for time-varying systems[END_REF], where different variants of transformations of time-varying systems to nonnegative forms are presented. Consider a LTV system described by:

⎧ ⎨ ⎩ ẋ(𝑡) = 𝐴(𝑡)𝑥(𝑡) + 𝑑(𝑡), 𝑦(𝑡) = 𝐶(𝑡)𝑥(𝑡) + 𝑣(𝑡), 𝑥(𝑡 0 ) ∈ R 𝑛 , 𝑡 ≥ 𝑡 0 ≥ 0, (20)
where 𝑥(𝑡) ∈ R 𝑛 , 𝑑(𝑡) ∈ R 𝑛 , 𝑦(𝑡) ∈ R 𝑝 and 𝑣(𝑡) ∈ R 𝑝 are respectively the state vector, an unknown but bounded input, the output vector and a bounded noise. Let 𝑥(𝑡 0 ) ∈ [𝑥 0 , 𝑥 0 ] for some known 𝑥 0 , 𝑥 0 ∈ R 𝑛 ; 𝑑(𝑡) ∈ [𝑑(𝑡), 𝑑(𝑡)] for all 𝑡 ≥ 𝑡 0 , where 𝑑, 𝑑 ∈ ℒ 𝑛 ∞ ; 𝑣(𝑡) ∈ [𝑣(𝑡), 𝑣(𝑡)] for all 𝑡 ≥ 𝑡 0 , where 𝑣, 𝑣 ∈ ℒ 𝑝 ∞ . In [START_REF] Thabet | An effective method to interval observer design for time-varying systems[END_REF], an effective technique has been proposed to build an interval observer for systems described by [START_REF] Raïssi | Interval observers design for consistency checks of nonlinear continuous-time systems[END_REF]. Assumption 4 is a conventional requirement for LTV systems [START_REF] Amato | Necessary and sufficient conditions for quadratic stability and stabilizability of uncertain linear time-varying systems[END_REF]. Under this assumption, the observer gain 𝐿(𝑡) and the matrix function 𝑀 (𝑡) are such that the stability of the LTV system ẋ(𝑡) = 𝐷(𝑡)𝑥(𝑡) can be proven by taking 𝑉 (𝑡) = 𝑥(𝑡) 𝑇 𝑀 (𝑡)𝑥(𝑡) as a Lyapunov function. It determines the output stabilization conditions of the system (20) which can be rewritten as:

{

ẋ(𝑡) = 𝐷(𝑡)𝑥(𝑡) + d(𝑡), 𝑦(𝑡) = 𝐶(𝑡)𝑥(𝑡) + 𝑣(𝑡) [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF] with d(𝑡) = 𝑑(𝑡) -𝐿(𝑡)𝑣(𝑡) + 𝐿(𝑡)𝑦(𝑡). When the gain 𝐿(𝑡) is computed such that the closed loop matrix 𝐷(𝑡) is stable and Metzler, the observer design is similar to the ones for LTI and LPV systems. In [START_REF] Thabet | An effective method to interval observer design for time-varying systems[END_REF], such a restrictive condition has been avoided. The methodology is based on the D-similarities approach developed in [START_REF] Zhu | Unified canonical forms for linear time-varying dynamical systems under D-similarity transformations. Part I[END_REF][START_REF] Zhu | Unified canonical forms for linear time-varying dynamical systems under D-similarity transformations. Part II[END_REF][START_REF] Zhu | Unified canonical forms for matrices over a differential ring[END_REF], where any time-varying matrix can be transformed under a Metzler form.

Proposition 2. There exists a time-varying transformation 𝑧 = 𝑇 (𝑡)𝑥 transforming 𝐷(𝑡) into a Metzler matrix Γ(𝑡):

Γ(𝑡) = 𝑇 (𝑡) ( 𝐷(𝑡)𝑇 -1 (𝑡) -𝑑(𝑇 -1 (𝑡))/𝑑𝑡 ) ( 22
)
where 𝑇 (𝑡) = 𝑃 (𝑡)𝐿 -1 (𝑡) with

Γ(𝑡) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 𝜆 1 (𝑡) 1 • • • 0 0 𝜆 2 (𝑡) . . . . . . . . . . . . . . . 1 0 • • • 0 𝜆 𝑛 (𝑡) ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 23
)
The elements 𝜆 𝑖 (𝑡) are called Essential D-eigenvalues of 𝐷(𝑡) (or ED-eigenvalues, see Definition 3.1 in [START_REF] Zhu | Unified canonical forms for linear time-varying dynamical systems under D-similarity transformations. Part I[END_REF]).

Theorem 4. [START_REF] Thabet | An effective method to interval observer design for time-varying systems[END_REF] Let assumption 4 be satisfied for [START_REF] Raïssi | Interval observers design for consistency checks of nonlinear continuous-time systems[END_REF] and ‖𝑦(𝑡)‖ ≤ 𝑌 for all 𝑡 ≥ 𝑡 0 for a known constant 𝑌 > 0.

Given the matrix 𝑇 defined in Proposition 2, assume that

∃ 𝑀 1 ∈ R + such that ‖𝑇 (𝑡)‖ + ‖𝑇 -1 (𝑡)‖ ≤ 𝑀 1 for all 𝑡 ≥ 𝑡 0 . Then, ż(𝑡) = Γ(𝑡)𝑧(𝑡) + 𝑑 𝑜𝑏𝑠 (𝑡) + Ψ 𝑜𝑏𝑠 (𝑡) + 𝑇 𝑜𝑏𝑠 (𝑡)𝑦(𝑡) ż(𝑡) = Γ(𝑡)𝑧(𝑡) + 𝑑 𝑜𝑏𝑠 (𝑡) + Ψ 𝑜𝑏𝑠 (𝑡) + 𝑇 𝑜𝑏𝑠 (𝑡)𝑦(𝑡)
is an interval observer for [START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF] and

𝑧(𝑡) ≤ 𝑇 (𝑡)𝑥(𝑡) ≤ 𝑧(𝑡), ∀𝑡 ≥ 𝑡 0 , where 𝑑 𝑜𝑏𝑠 (𝑡) = 𝑇 + (𝑡)𝑑(𝑡)-𝑇 -(𝑡)𝑑(𝑡), 𝑇 𝑜𝑏𝑠 (𝑡) = 𝑇 (𝑡)𝐿(𝑡), Ψ 𝑜𝑏𝑠 (𝑡) = 𝑇 - 𝑜𝑏𝑠 (𝑡)𝑣(𝑡) -𝑇 + 𝑜𝑏𝑠 (𝑡)𝑣(𝑡), 𝑑 𝑜𝑏𝑠 (𝑡) = 𝑇 + (𝑡)𝑑(𝑡) - 𝑇 -(𝑡)𝑑(𝑡), Ψ 𝑜𝑏𝑠 (𝑡) = 𝑇 - 𝑜𝑏𝑠 (𝑡)𝑣(𝑡) -𝑇 + 𝑜𝑏𝑠 (𝑡)𝑣(𝑡).

Case of switched systems

Consider a continuous-time switched system described by:

{ ẋ (𝑡) = 𝐴 𝜎(𝑡) 𝑥 (𝑡) + 𝐵 𝜎(𝑡) 𝑢 (𝑡) + 𝑤 (𝑡) 𝑦 𝑚 (𝑡) = 𝐶 𝜎(𝑡) 𝑥 (𝑡) + 𝑣 (𝑡) , 𝜎 (𝑡) ∈ ℐ (24)
where 𝑥 ∈ R 𝑛 , 𝑢 ∈ R 𝑚 , 𝑦 𝑚 ∈ R 𝑝 , 𝑤 ∈ R 𝑛 , 𝑣 ∈ R 𝑝 are the state vector, the input, the output, the disturbance and the measurement noise respectively. 𝑞 = 𝜎 (𝑡) is the index of the active subsystem and assumed to be known. 𝐴 𝑞 , 𝐵 𝑞 and 𝐶 𝑞 are constant matrices of appropriate dimensions.

The aim is to derive two signals 𝑥 (𝑡) and 𝑥 (𝑡) such that 𝑥 (𝑡) ≤ 𝑥 (𝑡) ≤ 𝑥 (𝑡) , ∀𝑡 ≥ 0 holds despite the disturbances and uncertainties provided that 𝑥 0 ≤ 𝑥 0 ≤ 𝑥 0 is satisfied. As in the previous sections, assume that the measurement noise and the state disturbances are assumed to be unknown but bounded with a priori known bounds such that -𝑤 (𝑡) ≤ 𝑤 (𝑡) ≤ 𝑤 (𝑡) , and |𝑣 (𝑡)| ≤ 𝑉 𝐸 𝑝 , ∀𝑡 ≥ 0, where 𝑤 : R + → R 𝑛 , 𝑤 (𝑡) ∈ ℒ 𝑛 ∞ and 𝑉 is a positive scalar.

As in [START_REF] He | Interval state observer for nonlinear switched systems with average dwell time[END_REF][START_REF] He | Control of non-linear switched systems with average dwell time: interval observerbased framework[END_REF], assume that there exist gains 𝐿 𝑞 such that the matrices (𝐴 𝑞 -𝐿 𝑞 𝐶 𝑞 ) are Metzler for all 𝑞 ∈ ℐ, where the matrices 𝐿 𝑞 (𝑞 ∈ ℐ) denote the observer gains associated with each subsystem 𝑞.

A candidate interval observer structure to compute 𝑥 and 𝑥 is described by:

{ ẋ = (𝐴 𝑞 -𝐿 𝑞 𝐶 𝑞 ) 𝑥 + 𝐵 𝑞 𝑢 + 𝑤 + 𝐿 𝑞 𝑦 𝑚 + |𝐿 𝑞 | 𝑉 𝐸 𝑝 ẋ = (𝐴 𝑞 -𝐿 𝑞 𝐶 𝑞 ) 𝑥 + 𝐵 𝑞 𝑢 -𝑤 + 𝐿 𝑞 𝑦 𝑚 -|𝐿 𝑞 | 𝑉 𝐸 𝑝 (25)
The following theorem gives the conditions for achieving the desired design goal. Theorem 5. [START_REF] Ethabet | Interval observers design for continuous-time linear switched systems[END_REF] [START_REF] Mazenc | Robust interval observers and stabilization design for discretetime systems with input and output[END_REF] then, the solutions of ( 25) and ( 24) satisfy 𝑥, 𝑥 ∈ ℒ 𝑛 ∞ and the inclusion 𝑥(𝑡) ≤ 𝑥(𝑡) ≤ 𝑥(𝑡)

provided that 𝑥 0 ≤ 𝑥 0 ≤ 𝑥 0 , where 𝐿 𝑞 = 𝑆 -1 𝑊 𝑞 .

Remark 1. The system ( 25) is initialized with the initial conditions 𝑥 0 and 𝑥 0 for the first active subsystem. At the switching time instant 𝑡 𝑖+1 , the output of the previous active subsystem (𝑞 = 𝜎(𝑡 𝑖 )) is used to initialize [START_REF] Mazenc | Interval observers for discrete-time systems[END_REF] with the subsystem (𝑞 = 𝜎(𝑡 𝑖+1 )).

In the LMI [START_REF] Mazenc | Robust interval observers and stabilization design for discretetime systems with input and output[END_REF] only the stability of the interval observer (25) is ensured. As in the case of non-switched systems, this LMI can be reformulated to find observer gains 𝐿 𝑞 minimizing the total interval estimation error 𝑥 -𝑥 as in [START_REF] Chebotarev | Interval observers for continuous-time lpv systems with l1/l2 performance[END_REF][START_REF] Briat | Interval peak-topeak observers for continuous-and discrete-time systems with persistent inputs and delays[END_REF][START_REF] Ellero | Unknown input interval observer with ℎ∞ and dstability performance[END_REF]. The main limitation of the approach detailed in this section is the existence of gains 𝐿 𝑞 such that the matrices (𝐴 𝑞 -𝐿 𝑞 𝐶 𝑞 ) are Metzler for all 𝑞 ∈ ℐ. To overcome this problem, some techniques which consists in finding a change of coordinates that transforms the observation errors into cooperative forms are proposed in [START_REF] Ethabet | Interval observers design for continuous-time linear switched systems[END_REF][START_REF] Rabehi | Interval estimation for linear switched system[END_REF][START_REF] Ethabet | [END_REF]. The changes of coordinates proposed for instance in [START_REF] Mazenc | Interval observers for linear time-invariant systems with disturbances[END_REF][START_REF] Raïssi | Interval state estimation for a class of nonlinear systems[END_REF] for continuous systems can be used to transform the matrices (𝐴 𝑞 -𝐿 𝑞 𝐶 𝑞 ) into a Metzler form.

Conclusion

Some recent results on interval observers for continuoustime systems have been presented in this overview. Interval observers can be considered as conventional pointwise observers with an estimate and an uncertainty quantification respectively given by the midpoint and the radius of each interval. This overview can be completed in a further paper by some results on discrete-time systems [START_REF] Efimov | On interval observer design for discrete systems[END_REF][START_REF] Efimov | Interval observers for time-varying discrete-time systems[END_REF][START_REF] Mazenc | Interval observers for discrete-time systems[END_REF][START_REF] Mazenc | Robust interval observers and stabilization design for discretetime systems with input and output[END_REF][START_REF] Briat | Interval peak-topeak observers for continuous-and discrete-time systems with persistent inputs and delays[END_REF][START_REF] Zhang | Event-triggered fault detection for a class of discrete-time linear systems using interval observers[END_REF], joint unknown and input estimation [START_REF] Gucik-Derigny | A note on interval observer design for unknown input estimation[END_REF][START_REF] Robinson | Interval observer design for unknown input estimation of linear time-invariant discrete-time systems[END_REF] and applications of interval observers in the field of robust control [START_REF] Efimov | Control of nonlinear and lpv systems: interval observer-based framework[END_REF][START_REF] Mazenc | Robust interval observers and stabilization design for discretetime systems with input and output[END_REF][START_REF] Cai | Stabilisation for a class of non-linear uncertain systems based on interval observers[END_REF][START_REF] Rauh | Intervalbased implementation of robust variable-structure and backstepping controllers of single-input singleoutput systems[END_REF][START_REF] Lamouchi | Interval observer framework for fault-tolerant control of linear parameter-varying systems[END_REF], diagnosis [START_REF] Puig | Passive robust fault detection using interval observers: Application to the DAMA-DICS benchmark problem[END_REF][START_REF] Blesa | Fault detection using interval LPV models in an open-flow canal[END_REF][START_REF] Blesa | FDI and FTC of wind turbines using the interval observer approach and virtual actuators/sensors[END_REF][START_REF] Thabet | Set-membership fault detection under noisy environment with application to the detection of abnormal aircraft control surface positions[END_REF], prognosis [START_REF] Yousfi | Setmembership methodology for model-based prognosis[END_REF] and automotive [START_REF] Efimov | Application of interval observers to estimation and control of air-fuel ratio in a direct injection engine[END_REF][START_REF] Goffaux | Continuous-discrete confidence interval observerapplication to vehicle positioning[END_REF].

In addition, sliding mode control and other estimation techniques are well known for their compensation of matched disturbances and finite-time convergence. Combination of sliding mode differentiators and interval observers is investigated in [START_REF] Efimov | Interval estimation for LPV systems applying high order sliding mode techniques[END_REF]. More generally, the combination of standard and interval observers can be considered as an appealing research issue for future works. Furthermore, a lot of open problems are still open in the field of interval observers design, namely for time-delay, distributed and hybrid systems. The optimization of the width of interval estimation remains also an interesting issue in the field of linear parameter varying systems.

  Fig. 1:The results of simulations for an LPV system[START_REF] Chebotarev | Interval observers for continuous-time LPV systems with 𝑙 1 /𝑙 2 performance[END_REF]

Assumption 4 .

 4 There exist bounded matrix functions𝐿 : R → R 𝑛×𝑝 , 𝑀 : R + → R 𝑛×𝑛 , 𝑀 (•) = 𝑀 (•) T ≻ 0 such that for all 𝑡 ≥ 𝑡 0 , Ṁ (𝑡) + 𝐷(𝑡) 𝑇 𝑀 (𝑡) + 𝑀 (𝑡)𝐷(𝑡) ≺ 0, 𝐷(𝑡) = 𝐴(𝑡) -𝐿(𝑡)𝐶(𝑡).

  For two vectors 𝑥 1 , 𝑥 2 ∈ R 𝑛 or matrices 𝐴 1 , 𝐴 2 ∈ R 𝑛×𝑛 , the relations 𝑥 1 ≤ 𝑥 2 and 𝐴 1 ≤ 𝐴 2 are understood elementwise. The relation 𝑃 ≺ 0 (𝑃 ≻0) means that the matrix 𝑃 ∈ R 𝑛×𝑛 is negative (positive) definite. Given a matrix 𝐴 ∈ R 𝑚×𝑛 , define 𝐴 + = max{0, 𝐴}, 𝐴 -= 𝐴 + -𝐴 (similarly for vectors) and denote the matrix of absolute values of all elements by |𝐴| = 𝐴 + + 𝐴 -. Let 𝑥 ∈ R 𝑛 be a vector variable, 𝑥 ≤ 𝑥 ≤ 𝑥 for some 𝑥, 𝑥 ∈ R 𝑛 .(1) If 𝐴 ∈ R 𝑚×𝑛 is a constant matrix, then

	Lemma 1. [22]

  Given 𝑥 0 , 𝑥 0 ∈ R 𝑛 . Let Assumption 1 be satisfied and 𝑥 ∈ ℒ 𝑛 ∞ . If there exist 𝜆, 𝜂 ∈ R 𝑛 , 𝜆 > 0, 𝜂 > 0, matrices 𝑊 𝑞 ∈ R 𝑛×𝑛 and positive scalars 𝛼 𝑞 , ∀𝑞 ∈ ℐ such that 𝑆𝐴 𝑞 -𝐶 𝑇 𝑞 𝑊 𝑇 𝑞 -𝑊 𝑞 𝐶 𝑞 + 𝛼 𝑞 𝑆 ≺ 0 𝐴 𝑇 𝑞 𝑆 -𝐶 𝑇 𝑞 𝑊 𝑇 𝑞 + 𝑑𝑖𝑎𝑔 (𝜂) ≥ 0 𝑆 = 𝑑𝑖𝑎𝑔 (𝜆)

	⎧ ⎨	𝐴 𝑇 𝑞 𝑆 +
	⎩


Pauline Puteaux
email: pauline.puteaux@lirmm.fr

Dave Trinel
email: dave.trinel@lirmm.fr

William Puech
email: william.puech@lirmm.fr

High-Capacity Data Hiding in Encrypted Images using MSB Prediction

Keywords: Image pre-processing, image encryption, image recovery, data hiding, MSB prediction

In the last few years, visual privacy has become a major problem. Because of this, encrypted image processing has received a lot of attention within the scientific and business communities. Data hiding in encrypted images (DHEI) is an effective technique to embed data in the encrypted domain. The owner of an image encrypts it with a secret key and it is still possible to embed additional data without knowing the original content nor the secret key. This secret message can be extracted and the initial image can be recovered in the decoding phase. Recently, DHEI has become an investigative field, but the proposed methods do not allow a large amount of embedding capacity. In this paper, we present a new method based on the MSB (most significant bit) prediction. We suggest to hide one bit per pixel by pre-processing the image to avoid prediction errors and, thereby, to improve the quality of the reconstructed image. We have applied our method to various images and, in every cases, the obtained image is very similar to the original one in terms of PSNR or SSIM.

I. INTRODUCTION

Reversible data hiding (RDH) consists of embedding secret data in a signal (e.g. an image). After its extraction, it is fundamental to reconstruct the original image with a minimum of errors or preferably, none at all. Indeed, in some strict areas, like in the military or medical world, distortion is unacceptable: each bit of information is important.

Fridrich et al. were the first to describe a RDH method [START_REF] Fridrich | Reliable detection of LSB steganography in color and grayscale images[END_REF]. They suggested compressing the least significant binary (LSB) planes of an image and to insert the to-be-embedded message in the vacating room. Although easy to implement and with good quality for the reconstructed image, this technique is not robust to steganalysis. Thereafter, a lot of schemes were proposed. Zhang et al. suggested exploiting the set of modification direction for a pixel (EMD) [START_REF] Zhang | Efficient steganographic embedding by exploiting modification direction[END_REF]. Tian et al. presented their method of difference expansion (DE) [START_REF] Tian | Reversible data embedding using a difference expansion[END_REF], [START_REF] Tian | Reversible watermarking by difference expansion[END_REF], using the Haar transform. It consists of calculating the differences between one pixel and its adjacent neighbors values and to select some of them to define the DE, where all the information is dissimulated. Methods based on histogram modification have also been described. Some proposed to build and to exploit the histogram according to the grayscale values [START_REF] Ni | Reversible data hiding[END_REF], [START_REF] Tai | Reversible data hiding based on histogram modification of pixel difference[END_REF] and others by using statistical data [START_REF] Gao | Lossless data embedding using GSQ histogram[END_REF], [START_REF] Tsai | Reversible image hiding scheme using predictive coding and histogram shifting[END_REF].

In RDH, the challenge lies in finding the best trade-off between the embedding capacity (in bpp) and the reconstructed image quality (in terms of PSNR or SSIM). Recently, for the purpose of answering to this problem, new methods based on prediction error analysis (PE) and their expansion (PEE) have emerged. The main idea is to exploit the correlation between a pixel and its adjacent neighbors [START_REF] Li | Efficient reversible watermarking based on adaptive prediction-error expansion and pixel selection[END_REF], [START_REF] Li | Reversible data hiding based on multiple histograms modification[END_REF], [START_REF] Ou | Pairwise prediction error expansion for efficient reversible data hiding[END_REF], [START_REF] Peng | An adaptive PEE-based reversible data hiding scheme exploiting referential prediction-errors[END_REF], [START_REF] Sachnev | Reversible watermarking algorithm using sorting and prediction[END_REF], [START_REF] Thodi | Expansion embedding techniques for reversible watermarking[END_REF], [START_REF] Wan | Reversible data hiding scheme based on prediction error sorting and double prediction[END_REF].

Otherwise, for data privacy, it is sometimes necessary to make an image unreadable. For this reason, a lot of encryption methods exist. These can be divided into two groups, depending if a block cipher or a stream cipher is used. Sometimes, in the second group, a new kind of algorithms, based on chaos, has been designed [START_REF] Baranovsky | Design of one-dimensional chaotic maps with prescribed statistical properties[END_REF], [START_REF] Kocarev | Chaos-based cryptography: a brief overview[END_REF], [START_REF] Wang | Design of pseudo-random bit generator based on chaotic maps[END_REF]. Moreover, in image processing, encryption can be full or selective.

In some cases, it may be interesting to combine image encryption and RDH (RDHEI). For example, in the medical community, a radiologist, bound by professional secrecy, has to encrypt a patient's X-rays. It must be possible to insert some information about the patient in these images without knowing their original content. Methods were also proposed to overcome this problem. Some suggested vacating room to embed data after the encryption phase (VRAE), others reserving room before image encryption (RRBE). In addition, encryption and data hiding can be joint, when data extraction and image reconstruction occur at the same time, or separate.

Puech et al. proposed to encrypt an image before embedding the secret message by using AES [START_REF] Puech | A reversible data hiding method for encrypted images[END_REF]. After that, they embedded a bit of the hidden message at a randomly selected position in each block. The original image was reconstructed by comparing values of the standard deviation in each block, calculated by assuming that the value of the hidden bit is 0 or 1. Zhang et al. method [START_REF] Zhang | Reversible data hiding in encrypted image[END_REF] consisted of image encryption, with a XOR operation, and data insertion by changing a small part. The image was divided into blocks with the same size and partitioned it into two groups. The three LSB of each pixel were modified in one of these groups, according to the to-be-inserted bit. Hong et al. [START_REF] Hong | An improved reversible data hiding in encrypted images using side match[END_REF] improved this technique by considering the local complexity for each pixel during the encryption step. Zhang [START_REF] Zhang | Separable reversible data hiding in encrypted images[END_REF] proposed a separate method: a part of the encrypted image was compressed and the vacated space served to embed secret data. Data extraction could be done before or after image decryption. Ma et al. were the first to introduce a RRBE technique [START_REF] Ma | Reversible data hiding in encrypted images by reserving room before encryption[END_REF]: a part of the original image was released before encryption, by using a RDH method of histogram shifting and some LSB were substituted in the encrypted image to embed data. Zhang et al. analyzed the prediction errors (PE) of some pixels and made space to hide data by PE-histogram shifting before image encryption [START_REF] Zhang | Reversibility improved data hiding in encrypted images[END_REF]. Wu and Sun also proposed a joint and a separate methods [START_REF] Wu | High-capacity reversible data hiding in encrypted images by prediction error[END_REF]. In the first one, they encrypted the original image in the same way as Zhang in [START_REF] Zhang | Separable reversible data hiding in encrypted images[END_REF]. After that, they partitioned the encrypted image to select pixels to embed data, according to a data hiding key. Then, they used a histogram shifting method. In the separative method, the to-be-inserted bits were hidden by MSB substitution. During the decoding phase, these bits could be extracted with the data hiding key and the original image was reconstructed by using a median filter on the watermarked image.

None of these methods succeeds in combining high payload embedding and high visual quality. Indeed, some of them are considered as reversible though PSNR is not equal to +∞. In [START_REF] Ma | Reversible data hiding in encrypted images by reserving room before encryption[END_REF], the payload can be high (0.5 bpp), but the reconstructed image is altered when compared with the original one (PSNR ≈ 40 dB). Moreover, other methods, such as Wu and Sun, propose a "high" embedding capacity, but it is only possible to embed approximately 0.1 bit per pixel at most [START_REF] Wu | High-capacity reversible data hiding in encrypted images by prediction error[END_REF].

In this paper, we introduce a new data hiding method for encrypted images based on MSB prediction with a very high capacity. During a pre-processing step, some pixel values are modified to avoid prediction errors, without significantly altering image quality. Thanks to this, it is possible to hide one bit per pixel by MSB substitution.

The rest of this paper is organized as follows. The proposed method is described in Section 2. Experiment results are provided in Section 3. And finally, the conclusion is drawn and the future work is discussed in Section 4.

II. PROPOSED METHOD

In this section, we introduce our separate data hiding method in encrypted images. The encoding phase consists of three steps: image pre-processing, encryption and data hiding by MSB substitution, as shown in Fig. 1. The original image owner pre-processes and encrypts it by using the encryption key K e and another person embeds an additional message by using the data hiding key K w . For the decoding phase, there are three possible schemes. If the recipient has only the encryption key, they can only obtain the pre-processed image but not the embedded message. Conversely, if they only have the watermarking key, they can just extract the message. Obviously, when both the encryption and the watermarking keys are obtained, the recipient can extract the original message and reconstruct the pre-processed image. The overview of this decoding method is presented in Fig. 2.

A. Pre-processing

Our method proposes to hide the secret message by replacing each MSB value of the image pixels by one bit of the hidden message. Thereby, the original MSB values are lost during the data hiding phase and, to reconstruct the original image, it is important to be able to predict them without any errors.

p(i, j) p(i, j -1)

p(i -1, j) Fig. 3: Context of a pixel p(i, j).
Then, we propose to process the original image I in order to obtain an image I without any error prediction. For the prediction, for each current pixel p(i, j), we propose to use its two adjacent pixels p(i, j -1) and p(i -1, j), as illustrated Fig. 3. First, we calculate the average of the two values. Then, we calculate the difference between this average and p(i, j) but also between this average and the inverse value of p(i, j), which is (p(i, j) + 128) mod 256. If the difference with the inverse value is smaller than with p(i, j), there will be a prediction error in the decoding phase. To avoid this, we must change the pixel value p(i, j) to get p (i, j). The detailed method is presented in Algorithm 1.

Algorithm 1: Pre-processing algorithm.

Data: Original m × n image I Result: Pre-processed image I for i ← 0 to m do for j ← 0 to n do inv(i, j) ← (p(i, j) + 128)%256; if i = 0 and j = 0 then pred(i, j) ← p(i, j);

else if i = 0 then pred(i, j) ← p(i, j -1)
;

else if j = 0 then pred(i, j) ← p(i -1, j); else pred(i, j) ← p(i-1,j)+p(i,j-1) 2 ; if | pred(i, j) -p(i, j) | ≥ | pred(i, j) -inv(i, j) | then if p(i, j) < 128 then p (i, j) = pred(i, j) -63;
else p (i, j) = pred(i, j) + 63; else p (i, j) = p(i, j)

B. Image encryption

During this phase, the pre-processed image I is encrypted to obtain the encrypted image I e . Piecewise Linear Chaotic Map (PWLCM) is one of the simplest chaotic systems [START_REF] Baranovsky | Design of one-dimensional chaotic maps with prescribed statistical properties[END_REF], [START_REF] Wang | Design of pseudo-random bit generator based on chaotic maps[END_REF]. Only simple operations are needed for each iteration:

xi = F (xi-1) =        xi-1 × 1 p if 0 ≤ xi-1 < p, (xi-1 -p) × 1 0.5-p if p ≤ xi-1 < 0.5, F (1 -xi-1) else, (1)
where p ∈ [0, 0.5] and

x i ∈ [0, 1].
Use elements of the secret key K e = (p, x 0 ) as parameters of this chaotic generator. As shown in Fig. 4, a sequence of pseudo-random bits b(i, j) k is obtained and used to encrypt the original image, pixel by pixel:

p e (i, j) k = b(i, j) k ⊕ p (i, j) k , (2)
where 0 ≤ k < 8 and refers to the number of the bit in a pixel (from MSB to LSB) and ⊕ represents the XOR operation.

C. Data embedding

In the data embedding phase, we can embed data in the encrypted image even if we do not have the encryption key and, so, we cannot access the original image content. By using the data hiding key, the to-be-inserted message is encrypted. In this way, it is not possible to detect its presence after the embedding in the marked encrypted image. After that, we scan pixels of the encrypted image from left to right, then from top to bottom (S-order) and substitute the MSB of each pixel by one bit b l , with 0 ≤ l < m × n, of the secret message, except for the first pixel:

p ew (i, j) = b l × 128 + (p e (i, j) mod 128). (3)

D. Data extraction and image recovery

In this phase, three cases are considered: (1) the recipient has only the data hiding key, (2) the recipient has only the encryption key and (3) the recipient has both the encryption and the watermarking keys.

In the first case, the recipient scans the pixels in the Sorder from the marked encrypted image I ew , and extracts the MSB of each pixel, except for the first pixel. After that, they just need to use the data hiding key to obtain the clear text. However, they cannot retrieve the pre-processed image I .

In the second case, the recipient can reconstruct I by adopting the following approach:

1) Use the encryption key to generate the pseudo-random chaotic sequence. 2) Scan the pixels of the marked-encrypted image I ew in the S-order and for each pixel, retrieve the seven least significant bits (LSB) of p (i, j) by XORing the marked encrypted pixel value p ew (i, j) with the associated binary sequence in the pseudo-random chaotic stream.

Only the MSB value could be wrong.

3) Predict the MSB value:

• If the considered pixel is not located on the first line or column, average p (i, j -1) and p (i -1, j) values. On the other hand, just consider p (i, j -1), if it is on the first line, or p (i -1, j), if it is on the first column. Record this value as a predictor pred(i, j). • Consider p (i, j) MSB=0 and p (i, j) MSB=1 as the pixel value with MSB = 0 and MSB = 1, respectively. Note that there is a difference equal to 128 between these two values. • Calculate the absolute difference between pred(i, j)

and each of these two values; finally, the smaller value gives the original pixel value:

If | p (i, j) MSB=0 -pred(i, j) | < | p (i, j) MSB=1 -pred(i, j) | then p (i, j) = p (i, j) MSB=0 , Else p (i, j) = p (i, j) MSB=1 .
In the third case, if the receiver has the data hiding and encryption keys, he can extract the secret message and reconstruct the pre-processed image in the same way as explained previously.

III. EXPERIMENTAL RESULTS

For data hiding methods in encrypted images, we have to measure different performances: embedding rate, number of incorrect extracted bits and recovered image quality after data extraction. It is necessary to find a trade-off between all of these parameters.

To evaluate image quality, we used two indicators: peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]. We show the results of the proposed method applied on the test image Lena (512 × 512 pixels), Fig. 5.a and on another image with a lower number of errors (512 × 512 pixels), Fig. 6.a. For these two images, we used the secret key (p, x o ) = (0.123456789, 0.567894123). Fig. 5.b and Fig. 6.b show the location of the problematic pixels, i.e. pixels of the original image whose MSB would be badly predicted if their value is not changed. Histograms in Fig. 5.c and Fig. 6.c display the required pixels modifications to avoid these prediction errors. After this pre-processing step, the adapted images are obtained (Fig. 5.d and Fig. 6.d). In Fig. 5.e and Fig. 6.e, one can notice that the initial information is not visible anymore after the encryption and watermarking steps. Finally, after data extraction, reconstructed images are exactly the same as pre-processed images and very similar to original images, even when some pixel values were changed during the pre-processing step (Fig. 5.f and Fig. 6.f). In both cases, PSNR is high and SSIM is close to 1: PSNR = 48.67 dB, SSIM = 0.9998 and PSNR = 61.11 dB, SSIM = 0.9999 for Lena and the second image respectively.

We applied our method on 500 different 512 × 512 gray level images 1 . In all cases, the embedding rate is 1 bpp and all bits of the hidden message are correctly extracted. In 9.2% of cases, when there is no prediction error (i.e. all the original pixel values are below or up to 128), our method is totally reversible. In this case, the original image is recovered without any errors, as indicated by PSNR of +∞ and SSIM of 1. In other cases, the reconstructed image quality is high, as presented in Table 1. In order to generate Fig. 7, we randomly selected 100 images among the 500 tested and applied three methods: LSB substitution, naive MSB method (without pre-processing) and our proposed method. All of these methods have the same payload (1 bpp). However, ours allows to reconstruct images globally with a better visual quality. Indeed, by using the LSB substitution, PSNR is close to 51 dB in all cases. Note that it was easy to speculate this using with a simple probability calculation. With a naive MSB method, results are generally worse because there is a phenomenon of error propagation during the pixel values prediction, excepting when there is no prediction error (9.2% of the cases). Ignoring these 9.2% of cases where PSNR is equal to +∞, PSNR is often below 20 dB. Sometimes, it is a little better -when there are only a few errors -but there are artifacts in the reconstructed image for all the pixels with the bad MSB value (difference of 128 with the correct value). Finally, with our method, PSNR is high, as illustrated in Table 1 with an average of 54.84 dB for 500 images: we avoided all the prediction errors thanks to the pre-processing step. We also made some comparisons between our proposed method and two existing ones: Zhang's method [START_REF] Zhang | Separable reversible data hiding in encrypted images[END_REF] and Wu and Sun's method [START_REF] Wu | High-capacity reversible data hiding in encrypted images by prediction error[END_REF], according to which is explained in the paper [START_REF] Wu | High-capacity reversible data hiding in encrypted images by prediction error[END_REF]. We used very well known images such Lena, Baboon, Airplane and Lake. Note in Table 2 that our method has a very high payload (1 bpp), contrary to that obtained by Zhang and Wu and Sun (0.1563 bpp). Regarding the reconstructed image quality, none of these three methods is totally reversible in every cases. Only the Lena image is exactly the same as the original one by using Wu and Sun's method. The results of our method are better than those of Zhang. They are not quite as good as those of Wu and Sun, but very similar. In order to compare the recovered image quality exclusively, we chose to embed regularly one pixel every six in such a way as to have approximately the same embedding capacity (0.1667 bpp) as Zhang and Wu and Sun's methods. Our results are better than those of Zhang and Wu and Sun: PSNR is higher, except when the original image is perfectly reconstructed. Finally, note that the results of our method and Wu and Sun's one are similar in terms of image quality when we watermark regularly one pixel every two (payload = 0.5 bpp).

In conclusion, our method allows a very good trade-off between the embedding rate and the recovered image quality after data extraction. Table 2: Performance comparisons between Zhang's method [START_REF] Zhang | Separable reversible data hiding in encrypted images[END_REF], Wu and Sun's method [START_REF] Wu | High-capacity reversible data hiding in encrypted images by prediction error[END_REF] and our proposed method.

IV. CONCLUSION

In this work, we proposed a new data hiding method in encrypted images based on prediction with a very high capacity (1 bpp). Indeed, by replacing all the MSB in the image, it is possible to hide one bit per pixel. In addition to this excellent embedding capacity, the reconstructed image quality is high (SSIM close to 1, PSNR ≈ 55 dB).

Future work on this method is to improve the image recovery. Indeed, during the pre-processing step, we modify only the current pixel to avoid all prediction errors. In most cases, it is a small modification, but sometimes -when the inverse pixel value is very close to the predictor valueit is necessarily greater. To obtain the best possible image quality, it is better to slightly change some pixels rather than making a big modification of one value, according to the mean squared error calculation. So, it may be interesting to use a new predictor, similar to the MED predictor. Finally, our ultimate aim would be to find a fully reversible method with the same high payload.

Fig. 2 :

 2 Fig. 2: Overview of the decoding method.

Fig. 4 :

 4 Fig. 4: Encryption step.

Fig. 5 :

 5 Fig. 5: Experiment using our method, embedding rate = 1 bpp. a) Original image, b) Errors' location, number of errors = 448 (0.2%), c) Histogram of the estimated prediction errors, d) Pre-processed image, PSNR = 48.67 dB, e) Marked encrypted image, f) Reconstructed image, PSNR = 48.67 dB, SSIM = 0.9998.

Fig. 6 :

 6 Fig. 6: Experiment using our method, embedding rate = 1 bpp. a) Original image, b) Errors' location, number of errors = 174 (0.1%), c) Histogram of the estimated prediction errors, d) Pre-processed image, PSNR = 61.11 dB, e) Marked encrypted image, f) Reconstructed image, PSNR = 61.11 dB, SSIM = 0.9999.

Fig. 7 :

 7 Fig. 7: Image quality comparison between our method and other methods with the same payload (1 bpp).

Table 1 :

 1 Images quality measurements on a database of 500 images.

		Best case Worst case	Average
		(9.2%)
	Number of MSB prediction	0%	3.2%	0.3%
	errors in the original image
	PSNR (dB)	+∞	36.06	54.84
	SSIM	1	0.9966	0.9998

By using the image data base of BOWS-2: http://bows2. ec-lille.fr/


Pierre Cordesse
email: pierre.cordesse@polytechnique.edu

Angelo Murrone
email: angelo.murrone@onera.fr

Marc Massot
email: marc.massot@polytechnique.edu

Marc Massot Coupling

Coupling a hierarchy of diffuse interface models with kinetic-based moment methods for spray atomization simulations in cryogenic rocket engines

Coupling a hierarchy of diffuse interface models with kinetic-based moment methods for spray atomization simulations in cryogenic rocket engines

Introduction

In a cryogenic combustion chamber, the multi-scale and various physical phenomena are very complex and their interaction is a current research area. In particular, the primary atomization plays a crucial part in the way the engine works, thus must be thoroughly studied to understand its impact on high-frequency instabilities. The latter have been encountered in the past and can lead to critical damages of the rocket. Even though experimentations must be conducted to enable simulation validation and to understand the observed physical phenomena, predictive numerical simulations are mandatory, at least as a complementary tool to understand the physics and even more to conceive new combustion chambers and predict instabilities they may generate in a given configuration. In subcritical conditions, downside a coaxial injector, the liquid oxygen atomization by co-current gaseous hydrogen engenders three two-phase flow topologies : at the exit of the injector, the two phases are separated by an interface. Downstream a polydisperse spray of droplets is carried by the gaseous phase. In between, shear stress caused by strong velocity gradients tears the liquid core apart and ligaments are formed. This process is called primary atomization. The ligaments get thinner and thinner until they break into droplets. In this mixed region, the subscale physics and the topology of the flow are very complex. Direct numerical simulations of these two-phase flows in a real engine configuration are still out of reach, CPU needs being too high and defining the smallest scale of a two-phase mixture is still unclear. Therefore predictive numerical tools must be developed using reduced-order models. However great care must be taken on the choices of these models in order to both have solid mathematics properties leading to robust and accurate numerical simulations after a validation process. Two main strategies are encountered in the literature to build such reduced-order models :

-Coupled models: a first model is used in the separated-phase can be derived either by statistical averaging of the instantaneous Navier-Stokes equations for each phase [7] or by front tracking methods including some level of space filtering [14]. In the dispersed-phase zone, the polydispersity in size, velocity and temperature of the droplets is taken into account with different methods: tracking the particles in a Lagrangian way [START_REF] Zamansky | [END_REF], or using an Eulerian approach where the droplets distribution is rebuilt thanks to the method of moments [23]. Usually, the methods applied to the separated phase zone are extended to the mixed region, but it implies either a high level of the phase disequilibrium description or an extremely refined mesh. On a top that, a numerical strategy must be provided to communicate informations between the different models.

-Unified models: In this strategy, a unique model takes care of all the flow topologies without the need of model coupling methods. A recent proposal in this spirit can be found in [8], where a unified model accounting for micro-inertia and micro-viscosity associated to bubble pulsation is proposed. The objective is here to model the sub-scale using a unified description, which degenerate into a predictive spray model in the disperse flow area [11] [10].

The present work aims at using the first approach by coupling a diffuse interface method to a member of the kinetic based moment methods (KBMM). The two models transport different variables. Therefore they need to be coupled in some way, transporting the information from one zone of the domain to the other. A reference implementation of the coupled approach can be found in [16], but the prediction of the polydispersity of the generated spray has to rely on some parameters impacting the resulting structure of the spray. The two-fluid model used, called a fourequation model is not capable of describing the phase disequilibria in terms of temperature, velocity and pressure.

Even if simulations including a temperature disequilibrium have already been conducted on real configurations, it seems it has never been done with also a velocity disequilibrium [21]. Figure 1 illustrates the numerical strategy applied in this work.

Figure 1: Topologies of a cryogenic jet atomization in sub-critical conditions.

Dispersed phase

Spray of droplets.

Evaporation and combustion.

Gas liquid interface.

Multi-species compressible �low.

Mixed zone

Primary atomization, ligaments.

Secondary fragmentation.

Separated phases

Diffuse interface method. Kinetic-Based-Moment Methods.

VL=VG

VL≠VG

The paper is organized as follows. First, a mathematical overview of the various models used is adressed. The numerical methods applied to these models are described in the following section. Later on details of the numerical simulations conducted are given and a comparison between the different models is drawn to finally conclude on the work.

Mathematical modeling

Two models are used: a two-phase flow model based on diffuse interface methods and a multi-fluid model derived by Kinetic Based Moment Methods. They are two-way coupled.

Hierarchy of diffuse interface methods

Diffuse interface methods rely on the assumption that at any time t, at any place x in a given domain Ω, two phases can coexist. Hence, three configurations occur : time-space point (t, x) with only one or the other phase or with both of them. The latter forms a mixture. To obtain a continuum-mechanical approach for describing two phase flows, a statistical averaging process is applied on the microscopic Navier-Stokes equations for each phase [7]. A system of conservative equations is obtained allowing full disequilibrium between the two phases. Two assumptions are made to obtain these conservation equations. First, each phase behaves independently to one another, the exchanges between the phases are modeled by additional mass, momentum and energy exchange terms. Second, the mixture follows the same conservation equations obtained by summing those of each phase. Then, the entropy inequality derived from the second law of thermodynamics imposes restrictions on the constitutive equations, resulting in models for each exchange term. However the system of six equations, namely a mass, momentum and energy equations for each phase, is not closed. Baer and Nunziato introduced for the first time a supplementary evolutionary equation for the volume fraction of one phase to ensure the model closure, resulting in a so-called seven-equation model describing the deflagration-to-detonation (DDT) in gas-permeable, reactive granular materials [2]. This model has been generalized for two-phase flow by introducing two interfacial quantities, u I and p I , the averaged interfacial velocity and pressure respectively [19], and writes in a quasi-linear form :

∂ t U + A 1 ∂ x U = R(U ) + S u→q -S drag -S h (1)
with the conservative variables

U = (α 2 , U 2 , U 1 ) t , U k = (α k ρ k , α k ρ k u k , α k ρ k E k )
, the conservative fluxes F (U ) and the non-conservative fluxes G(U ) defining the matrix A 1 as:

A 1 =   u I 0 0 ∂ α2 F (U 2 ) + G(U 2 ) ∂ U2 F (U 2 ) 0 ∂ α2 F (U 1 ) + G(U 1 ) 0 ∂ U1 F (U 1 )   (2)
where

F (U k ) = α k ρ k u k , α k ρ k u 2 k + α k p k , α k (ρ k E k + p k )u k t and G(U 2 ) = -G(U 1 ) = (0, p I , p I u I )
t . α k is the volume fraction of phase k = 1, 2, ρ k the partial density, u k the phase velocity, p k the phase pressure,

E k = k + 1/2u 2
k the total energy per unit of mass, k the internal energy. A two-parameter equation of state associated with a Gibbs law will be used hereafter. Mass transfer between the two phases is neglected. S u→q models the atomization of the liquid phase into droplets and the pseudo coalescence of the droplets into the liquid phase, S drag the drag force of the gas acting upon the droplets, S h the conducto-convective heat transfer at the surface of the droplet.

Because a mixture can not stay in full disequilibrium, R is an application describing the mechanical, hydrodynamic and thermal relaxations between the two phases and is the characteristic time for each of these processes. The condition R(U ) = 0 R 7 imposes three constraints (pressure, velocity and temperature equilibria) and defines the constrained manifold

E = U ∈ R 7 , R(U ) = 0 R 7 .
Usually, the constraint R is defined by physical processes respecting the entropy inequality. When only pressure and velocity relaxations are accounted for, it decomposes classically onto

R = R u u + R p p , with R u u = 0, R u 2 u , R u 1 u t and R p p = p 2 -p 1 p , R p 2 p , R p 1 p t (3)
where

R u 2 = -R u 1 = (0, u 2 -u 1 , u I (u 2 -u 1 )), R p 2 = -R p 1 = (0, 0, p I (p 2 -p 1 )
). The interfacial terms u I and p I modeling varies and an expression for the relaxation parameters u and p have been for example derived using the DEM technic in [20]. Relaxing the pressure and the velocity can be interpreted as projecting the full disequilibrium state U ∈ R 7 on the constraint manifold E where the relaxation R(U ) = 0 R 7 occurs. One can thus look for a reduced model by projecting U ∈ R 7 on R p , p ≤ 7. It can be shown that there exists an admissible map M from R p to E ⊂ R 7 which leads to semi-stable systems. In particular, one can obtain the five-equation model where the phases have now the same pressure and the same velocity. In [18], the authors interpret the reduced model as the asymptotic limit of the seven-equation model when → ∞. Following their lines, without the coupling source terms, the five-equation model writes :

∂ t U + A 1 ∂ x U = R(U ) (4)
with the conservative variables

U = (α 2 , U 2 , U 1 ) t , U k = (α k ρ k , α k ρ k u, α k ρ k E k )
t , the conservative fluxes F (U k ) and the non-conservative fluxes G(U k ) defining the matrix A 1 as in Equation ( 2) but with u I = u 1 = u 2 = u, p I = p 1 = p 2 . Finally taking the Taylor series of the relaxing source terms with respect to gives

u I = u + o( ), p I = p + o( ) and u 2 -u 1 u = α 1 ρ 1 α 2 ρ 2 ρ ( 1 ρ 1 - 1 ρ 2 )∂ x p + o( ) p 2 -p 1 p = α 1 α 2 ρ 1 a 2 1 -ρ 2 a 2 2 k=1,2 α k ρ k a 2 k ∂ x u + o( ) (5)
where ρ is the mixture density [2]. The zero order terms are injected into Equation (3) to define R in Equation (4). In practice, for the five-equation model, the two partial momentum and energy equations are summed and replaced by an equation on the total momentum ρu and an equation on the total energy ρE respectively. Only one relaxing term remains in the system of equation: that of the volume fraction equation. Thus the only difference between the five-equation model and the instantaneously relaxed seven-equation model lies in the volume fraction equation.

ρ = α 1 ρ 1 + α 2 ρ 2 , a k is the phase sound of speed, a 2 k = ∂p k /∂ρ k | s k , and k = k + 1
When relaxing the temperatures, one obtains the compressible Navier-Stokes equations, called also four-equation model. This last model, associated with the two previous one, defines the hierarchy of diffuse interface models at disposal.

To accurately describe the mixed region, one needs to use the out-of-equilibrium seven-equation model to feed more accurately the KBMM for the velocity, temperature and pressure of the droplets.

Mathematical properties of the seven-equation model

The mathematical properties of the seven-equation model have been studied by [START_REF] Coquel | [END_REF] among others. It is hyperbolic and admits 7 eigenvalues Sp (A 1 ) = u I , {u k , u k ± a k } k=2,1 . As noticed by [9], the system becomes non strictly hyperbolic when two or more eigenvalues coincide. Assuming the interface velocity u I is defined as a linear average of the two phase velocities u I = βu 1 + (1 -β)u 2 , β ∈ [0, 1] then the eigenvalues λ 1 = u I , λ 2 = u 2 and λ 5 = u 1 remain distinct. Thus the only condition leading to hyperbolic degeneracy, designated as the non resonance condition in [5], is

(u k -u I ) 2 = a 2 k , k = 1, 2.
A KBMM element: sampling method

It would have been possible to choose one of the best KBMM element such as a multi-fluid modelling using a continuous discretization of the droplet size through sections, which has been validated on evaporating polydisperse sprays [23]. Nevertheless since our primary concern is to increase the disequilibrium in the interface diffuse model, for sake of simplicity, a simple KBMM element has been chosen: a multi-fluid modelling with sampling methods [15] and monokinetic and monotemperature assumptions. Evaporation and coalescence of the droplets are neglected, thus only one sample of droplet is needed to attest the success of the coupling strategy of the two models. The system of equation is weekly hyperbolic [6] and writes for one class:

∂ t Q + ∂ x F q (Q) = S drag + S h + S q→u (6)
where the variables are

Q = (n, α d ρ d , α d ρ d u d , α d ρ d d ) t the conservative fluxes F q (Q) = (nu d , α d ρ d u d , α d ρ d u 2 d , α d ρ d d u d ) t
, and the coupling source term S q→u . n is number of droplets, ρ d is the liquid density depending only on the temperature, u d the droplet velocity, d the droplet internal energy. To be conservative, the coupling source terms of Equation ( 1) and Equation ( 6) must balance out each other, thus S q→u = -S u→q .

Numerical methods

The models used hereafter to simulate the jet atomization have been implemented in the CEDRE software which is a multi-physics platform working on general unstructured meshes and organized as a set of solver [13]. Two solvers are used, SEQUOIA for the diffuse interface model and SPIREE for the KBMM. Both are two-way coupled by the source terms S q→u and S u→q exchanged at every time step. A Lie splitting technic is used resulting in the following system of equations:

U n+1 = S u→q -S h -S drag R u u ∆t R p p H u (U n ) Q n+1 = [S q→u S h S drag ] ∆t H q (U n ) (7)
The hyperbolic operators H u and H q corresponding to convection of the two-phase model and the dispersed model respectively, are calculated using an approximate Riemann solver, HLLC, and a Presureless Gas Dynamics exact Riemann solver respectively [START_REF] Boileau | ICMF Proceedings[END_REF]. Then relaxation operators and source terms are applied in the order showed in Equation ( 7) to define the new states U n+1 and Q n+1 of the conservative variables.

HLLC scheme for the two-phase flow model

Ideally, an exact Riemann solver should be able to solve the hyperbolic operator H u [22], however the computational cost is way too expensive. Thanks to hyperbolicity of the model, a VFRoe solver has been tested [1] using primitive variables U = (α 1 , ρ 1 , u 1 , p 1 , ρ 2 , u 2 , p 2 ) t . Nevertheless, it was not robust enough to handle the strong gradients occurring in the chamber. Therefore, choice has been made to use a HLLC solver which shows usually good robustness [24]. Other difficulties for solving Equation (1) are owed to the non-conservative terms and modelling the interfacial terms. The approach proposed in [12] tackles all the issues by assuming (1) the interfacial terms p I and u I to be local constants in the Riemann problem, (2) the volume fraction to vary only across the interfacial contact discontinuity u I . As a result, the non conservative terms in Equation (1) vanish, u I and p I are determined locally by Discrete Equation Method (DEM) [20] at each time step and stay constant during the update. Thus, phases are now decoupled and Equation (1) splits into two conservative sub-systems to which a classic HLLC solver is applied. As proposed by [12], a second order reconstruction is developed. It uses the MUSCL-Hancock strategy but the spatial reconstruction is done with a multi-slope MUSCL scheme developed by [17]. Finally, two slope limiters have been tested: the Van-Leer limiter and a hybrid limiter developped in [16], which combines the advantages of a CFL-Superbee limiter for high gradient regions and a third order CFL dependent limiter for the regular regions.

Relaxation procedure

In many applications, the pressures of a two-phase flow are assumed to relax instantaneously. Thus, after the hyperbolic update, the following ODE is solved: ∂ t U = R p (U )/ p , with p → ∞ which infers u k remain constant. Manipulating the equations, an equilibrium pressure is obtained by solving a second order equation with an iterative procedure such as a Newton method. Detailed of the equation can be found in [12]. As for the velocities, within the context of assisted jet atomization in sub-critical conditions, the relaxation time is finite. The following ODE is solved: ∂ t U = R u (U ) / u which implies α k , ρ k are conserved during the relaxation. Subtraction of the momentum equations gives the following ODE on the split velocity

u d = u 2 -u 1 : ∂ t u d -A o u d / u = 0
where superscript o denotes the state before relaxation. A first numerical approach is to fix a remaining slip velocity ratio target at each computational time step ∆t. It defines the characteristic relaxing time:

u A o = ln(X)∆t with X = u d (∆t) u o d and A o = α o 1 ρ o 1 + α o 2 ρ o 2 α o 1 ρ o 1 α o 2 ρ o 2 (8)
An instantaneous velocity relaxation is in pratice also possible and manipulating the ODE leads to a unique relaxed velocity, which is the mass weighted average of the two velocities before relaxing.

Coupling source terms

The physical processes accounted for in the coupling source terms are the atomization of the liquid phase, S atom , and the pseudo-coalescence of the liquid droplets, S coal , defined respectively as S atom = α 2 ρ 2 f atom λ atom and S coal = α d ρ d f coal λ coal where f atom is the atomization frequency, λ atom describes the efficiency of the atomization, f coal the pseudo-coalescence frequency and λ coal the pseudo-coalescence efficiency, defined all in [16].

Thanks to the increase of disequilibrium by using the seven-equation model with finite velocity relaxation, S atom has been revisited making full use of the existence of two velocities to track regions with high shear stress and thus regions where atomization should occur.

Results and Discussion

In the following, numerical simulations are performed on a subcritical cryogenic single-injector configuration using the hierarchy of diffuse interface models introduced previously, namely the five equation model (5eq), the instantaneously relaxed seven equation model (IR7eq) and the non-instantaneously relaxed seven equation model (NIR7eq). On the way to validation of the numerical results with experimental data, the first interest is to attest the robustness of the numerical method by increasing difficulties one step at a time: first increase the disequilibrium of the diffuse interface model then activates the coupling via atomization and pseudo-coalescence source terms. Doing so, we can study properly the impact of the models and the coupling on quantities of interest such as the liquid core length, its dynamics, the sharpness of the interface and the velocities at the interface. However, we choose operating conditions as close as possible to a real configuration, the results may thus also be predictive.

Description of the configuration

The configuration choice meets several criteria. First, the geometry mimics the experimental test-bench MAS-COTTE [START_REF] Vingert | EAC Proceedings[END_REF], a cryogenic rocket engine combustion chamber. It adopts a unique co-axial injector of liquid oxygen O 2 (l) circumscribed by gaseous hydrogen H 2 (g) . Second, it must offer three dimensions in space to capture the dynamics of the jet. However, the computational time should not be too heavy to conduct numerical tests and validations. Hence, only a portion θ ∈ [0, π/3] of the cylindrical chamber is meshed making use of the symmetric axis of the cylinder preventing the liquid jet to flap. As seen in Figure 2, the mesh is refined inside the injector and at its exit where the liquid core flows in order to capture the interface dynamics. The mesh size is around a million tetra cells. The injector lip of length L lip is meshed by four cells at the minimum mesh size ∆x min . Walls are set to adiabatic slip boundaries, the variables (ρ k , u k , T k , α k ) define the inlets, the outlet is pressure defined at p ∞ . Table 1 summarizes key values of the configuration. The numerical simulation has been conducted as followed: first, the oxygen and the hydrogen have been injected with a ramp-up to reach at τ = τ 0 the operating point using Figure 2: Geometry and mesh of the configuration Table 1: Physical parameters of the configuration

J = ρ 1 u 2 1 ρ 2 u 2 2 T 1 /T 2 Gas Liquid L lip /∆x min p ∞ α ∼ 3 ∼ 3 H 2 (g) O 2 (l) 4 10 bar 10 -6
the 5eq with no coupling. Then, the simulation has run approximately ten times the characteristic convective time of the liquid core τ conv . At this point, designated τ 1 , start the comparisons of the models.

Case 1 : comparison of the five-equation model and the instantaneously relaxed seven-equation model without atomization

The first case aims at comparing the 5eq with the IR7eq. The unique difference emphasized by the mathematical study lies on the order with respect to of the source term of the volume fraction equation. Therefore Figure 3 compares qualitatively the liquid core obtained after two τ conv simulation time starting at τ = τ 1 . Looking at To better describe the flow in the region close to the interface, it is necessary to increase the disequilibrium of the phases. In the present configuration, a strong velocity gradient occurs at the interface, and the ratio of kinetic energy, J = (ρ 1 u 2 1 )/(ρ 2 u 2 2 ), is about 3. Physically, the liquid is expected to be accelerated by the gas, but not instantaneously. Therefore it is physically wrong to assume an instantaneous relaxation time for the velocities of the phase. The same comment can be done on the temperatures. So far in the literature [21], it has been achieved to allow the phases to have their own temperature but always have a common velocity using for example the 5eq. Even when the 7eq. model has been used, it was always with an instantaneous relaxation time of the velocities mainly due to numerical difficulties [12]. In such case, the 7eq. model lacks of interest in practice compared to the 5eq.

Here for the first time, we have successfully conducted a simulation of jet atomization with a non-instantaneously relaxed seven-equation model. Each phase retains its own velocity in the region of the interface. The simulation has run for 2.5τ conv attesting the success of the implementation. in the interface region, on the so-to-say "gaseous side", where atomization occurs. It permits to give to the atomized liquid droplets the speed of the liquid phase. It is a major gain of accuracy as long as the characteristic time of velocity relaxation is physically well-defined. At the present time, the characteristic time is finite and constant, but it will be revisited in future works to match physical reality. Finally, a comparison between the IR7eq and no coupling and the NIR7eq with atomization is proposed hereafter. The reason why the former was not coupled with atomization in due to the fact that the atomization source term S u→q depends on the existence of a slip velocity u d . Nevertheless comparing these two models helps us to verify that the activation of the atomization does not destroy the liquid core. Figure 6 compares qualitatively the liquid interface after 1.5τ conv simulation time. Interestingly, the appearance is similar for the two models meaning the atomization process is not interfering with the liquid core, which is what was expected based on the choice of the efficiencies of S atom and S coal . Furthermore the length of the liquid core L lc is also very similar up to the point that the two models show the same time averaged length (Figure 4b) which confirms that the atomization procedure has been carefully implemented.

Conclusion and future prospects

An original fully Eulerian modelling strategy of jet atomization in sub-critical combustion chambers has been successfully implemented. First, it relies on the coupling of a hierarchy of diffuse interface models from a fiveequation model to a non-instantaneously relaxed seven-equation model with an eulerian model derived by kineticbased moment method (KBMM). Temperature and velocity disequilibrium levels for the diffuse interface model have been included to describes the separated and mixed zone. The numerical strategy has copped with the strong discontinuities encountered in jet atomization thanks to a robust and accurate numerical method using multi-slope MUSCL technique. All the models have captured the pulsating movement of the liquid core and the distortion of the interface due to high-shear stress. The seven-equation model when non-instantaneously relaxed, has permitted to feed the KBMM element with a liquid velocity.

In future works, we will first aim at defining a physical dynamic velocity relaxation time. Then, we will couple the diffuse interface model with a more accurate KBMM element such as multi-sectional models including coalescence and evaporation source terms to obtain a predictive spray of droplets and compare it to experimental data obtained on the test-bench MASCOTTE [START_REF] Vingert | EAC Proceedings[END_REF].

Figure 3 :

 3 Figure 3: Comparison of the liquid volume fraction α 2 : slice at θ=π/6, at τ =τ 1 (left) and τ =τ 1 +2τ conv (right), α 2 =1.0 α 2 =0.01

Case 2 :

 2 comparison of the instantaneously relaxed seven-equation model and the non-instantaneously relaxed seven-equation model with atomization

Figure 4 : 2 :

 42 Figure 4: Model influence on the liquid core length L lc over time and its time average L lc at isovalue α 2 =0.95±0.04. Figure 4a compares 5eq and IR7eq. Figure 4b compares NIR7eq and IR7eq

Figure 5 Figure 5 :

 55 Figure 5: NIR7eq with atomization (Legend: blue isovolume of α 2 = 0.99, volume fraction of droplets in the SPIREE solver α d , slip velocity norm |u slip |)

Figure 6 :

 6 Figure 6: Comparison of the liquid volume fraction α 2 : slice at θ=π/6, at τ =τ 1 (left) and τ =τ 1 +1.5τ conv (right), α 2 =1.0 α 2 =0.01 atomized NIR7eq

Acknowlegdments

This work has been cofunded by the French Aerospace Lab (ONERA) and the French Government Space Agency (CNES). Simulations have been successfully conducted using the CEDRE computational fluid dynamic software.