Cristobal Bertoglio

David Barber

Nicholas Gaddum

Israel Valverde

Marcel Rutten

Philipp Beerbaum

Philippe Moireau

Rodney Hose

Jean-Frédéric Gerbeau

J-F Gerbeau

Identification of artery wall stiffness: in vitro validation and in vivo results of a data assimilation procedure applied to a 3D fluid-structure interaction model

Keywords: fluid-structure interaction, data assimilation, stiffness estimation, medical imaging, experimental validation, coarctation of the aorta

teaching and research institutions in France or abroad, or from public or private research centers.

Identification of artery wall stiffness: in vitro validation

and in vivo results of a data assimilation procedure applied to a 3D fluid-structure interaction model

Introduction

The in vivo estimation of arterial stiffness can provide valuable information about the cardiovascular condition of a patient [START_REF] Laurent | Expert consensus document on arterial stiffness: methodological issues and clinical applications[END_REF]. In clinical practice, the Pulse Wave Velocity [START_REF] Nichols | McDonald's blood flow in arteries: theoretical, experimental and clinical principles[END_REF] gives an estimation of the average stiffness of a portion of the arterial network by solving the Moens-Korteweg equation, which assumes a linear solid in an infinite cylindrical domain. Other methods include Vascular Elastography, either based on Ultrasound [START_REF] De Korte | Vascular ultrasound for atherosclerosis imaging[END_REF][START_REF] Balocco | Noninvasive Young's modulus evaluation of tissues surrounding pulsatile vessels using ultrasound doppler measurement[END_REF][START_REF] Balocco | Estimation of the viscoelastic properties of vessel walls using a computational model and Doppler ultrasound[END_REF] or Magnetic Resonance [START_REF] Woodrum | Vascular wall elasticity measurement by magnetic resonance imaging[END_REF][START_REF] Xu | Assessment of stiffness changes in the ex vivo porcine aortic wall using magnetic resonance elastography[END_REF], where the models used for estimating arterial stiffness are again based on simplified models.

We can go further and obtain the local mechanical properties of a vessel if we adopt a data assimilation procedure, based on a three-dimensional fluid-structure interaction model and the distributed kinematic data obtained from dynamic 3D acquisitions (e.g., MRI or CT). In this way, stiffness values can be obtained under physiological conditions without requiring invasive pressure measurements.

Recently, [START_REF] Bertoglio | Sequential parameter estimation in fluidstructure problems. Application to hemodynamics[END_REF] proposed an efficient methodology to estimate uncertain parameters (for example Young's modulus or fluid boundary conditions) from measurements of the displacement of the wall in an idealized 3D fluid-structure interaction (FSI) problem. The procedure is based on the Reduced Order Unscented Kalman Filter (ROUKF) [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF], which consists of a tractable filtering that allows to solve the estimation problem with a computational effort comparable to one forward simulation. Note that variational (in space and time) approaches have been adopted by other authors but for simplified models in [START_REF] Lagrée | An inverse technique to deduce the elasticity of a large artery[END_REF][START_REF] Martin | Parameter identification for a one-dimensional blood flow model[END_REF][START_REF] Stalhand | Determination of human arterial wall parameters from clinical data[END_REF], and for threedimensional problems in [START_REF] D'elia | Applications of variational data assimilation in computational hemodynamics[END_REF][START_REF] Perego | A variational approach for estimating the compliance of the cardiovascular tissue: An inverse fluid-structure interaction problem[END_REF] but simplifying the minimization problem in order to avoid the resolution of the adjoint equations in time.

Whereas in [START_REF] Bertoglio | Sequential parameter estimation in fluidstructure problems. Application to hemodynamics[END_REF] validations are only based on synthetic data, we consider in the present work the experimental validation of this method using real data in a fluid-structure interaction context. The data come from a laboratory experiment performed on a silicone tube mimicking an aorta. The results show an excellent agreement between the values obtained from data assimilation and those obtained by independent mechanical tests.

We also apply the stiffness estimation methodology to a real clinical case, namely an aorta with a mild coarctation. Hence, this work complements the work of [START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF] since here we estimate the stiffness parameters and not the boundary conditions.

The rest of the paper is structured as follows. In Section 2 we present the FSI-equations, in Section 3 we summarize the ROUKF algorithm and in Section 4 we overview the image segmentation algorithm. In Section 5 we present the results for stiffness estimation with experimental data and we conclude in Section 6 with the aorta study.

The fluid-structure interaction equations

Mathematical model We deal with the numerical resolution and the data assimilation of the mechanical interaction between an incompressible fluid and an elastic structure. The fluid is described by the Navier-Stokes equations (NSE), in a moving domain Ω f (t) ⊂ R d , d = 2, 3, in an Arbitrary Lagrangian Eulerian (ALE) formulation, and the structure by the elastodynamic equations in Ω s (t) ⊂ R d . The fluid-structure interface is denoted by Σ = ∂Ω s ∩ ∂Ω f , and

∂Ω f = Γ in ∪ Γ 1 ∪ • • • ∪ Γ n 0D ∪ Σ and ∂Ω s = Γ n ∪ Γ d are
given partitions of the fluid and solid boundaries, respectively.

The coupled FSI problem reads as follows: for t > 0, find the fluid velocity u f (t) : Ω f (t) → R d , the fluid pressure p(t) : Ω f (t) → R, the structure displacement y s (t) : Ω s (t) → R d and structure velocity u s (t) : Ω s (t) → R d such that

• Fluid equations:              ρ f ∂u f ∂t ξ + ρ f (u f -w) • ∇u f -∇ • σ f (u f , p) = 0, in Ω f , ∇ • u f = 0, in Ω f , u f = u in , on Γ in , σ f (u f , p)n f = -P i n f , on Γ i , i = 1, . . . , n 0D , (1a) with σ f (u f , p) = -pI + 2µ f ε(u f ),
the Cauchy stress tensor of the fluid and ε(u f ) = (∇u f + (∇u f ) )/2 denotes the deformation rate tensor, µ f the dynamic viscosity, and ∂ ∂t ξ the ALE derivative (see e.g. [START_REF] Fernández | Algorithms for fluid-structure interaction problems[END_REF]). In the hemodynamics problems considered in this work, the outlet pressure P i is obtained by solving the differential-algebraic equation, known as the Windkessel model [START_REF] Peiró | Reduced models of the cardiovascular system[END_REF]   

P i = π i + R p,i Q i , C i dπ i dt + π i R d,i = Q i , Q i = Γ i u f • n f . (1b)
The distal resistance R d,i , the proximal resistance R p,i and the capacitance C i are assumed to be given.

• Structure equations:

         ∂ t y s = u s , in Ω s , ρ s ∂ t u s -η s ∇ • σ s (u s ) -∇ • σ s (y s , θ) = 0, in Ω s , y s = y in , on Γ d , η s σ s (u s )n s + σ s (y s )n s = -c Γ u s -k Γ y s , on Γ n , (1c)
where σ s is the Cauchy stress tensor of the solid and the vector θ ∈ R κ denotes the set of solid constitutive parameters that will be estimated later. In this work, both quantities are related by:

-Considering a linear constitutive relation, namely σ s (y s , θ) = C(θ) (y s ) with (y s ) = (∇y s + (∇y s ) )/2 and C the classical elasticity tensor, where the Young's modulus is given by E = E 0 • 2 θ , E 0 > 0 , while the Poisson's ratio is fixed. -Considering a nonlinear, hyperelastic, neo-hookean material model, such that the elastic energy density function is given by W = c 0 • 2 θ • (I 1 -3), c 0 > 0 and I 1 the first invariant of the right Cauchy-Green strain tensor.

Moreover, the parameters c Γ and k Γ model in a simple way the external tissue effect on the vessel of interest [START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF][START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF] and η s is a viscoelastic coefficient.

• Coupling conditions:

     y f = Ext f Σ(0) y s| Σ , w = ∂ t y f , Ω f (t) = I Ω f + y f (t) (Ω f ), u f = u s , on Σ, η s σ s (u s )n s + σ s (y s )n s + σ f (u f , p)n f = 0, on Σ, (1d)
where Ext f Σ is an extension operator from Σ to Ω f .

FSI numerical algorithm

The spatial discretization is performed with a first order stabilized finite element. For the time marching, we use a semi-implicit partitioned FSI-algorithm [START_REF] Fernández | A projection algorithm for fluidstructure interaction problems with strong added-mass effect[END_REF][START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF] with a first order Chorin-Temam projection method in the fluid, and a Newmark scheme in the solid. The coupling between the solid and the fluid is explicit in the ALEadvection-diffusion step and implicit for the pressure projection step and the Windkessel [START_REF] Bertoglio | Fractional-step schemes for the coupling of distributed and lumped models in hemodynamics[END_REF].

3 The data assimilation procedure

In this section we summarize the data assimilation procedure based on [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] and already studied in [START_REF] Bertoglio | Sequential parameter estimation in fluidstructure problems. Application to hemodynamics[END_REF] for idealized (but realistic) three-dimensional arterial FSI problems.

Inverse problem statement Assume that we can write the fully discrete FSI-problem as

X n+1 = A n (X n , θ), n ≥ 0, (2)
with X n = (u n f , y n s , u n s ) ∈ R q the discrete dynamical state, θ ∈ R κ the physical parameters, and X 0 the given inititial condition. Suppose that θ and X 0 are uncertain, and that we have (noised) measurements of real physical system Z n ∈ R z related to the numerical state by the observation error Γ n (Z n , X n ) = ζ n , with ζ n ∈ R m the noise. Data assimilation consists in reducing the uncertainties of a model by minimizing a cost function like (see, e.g., [START_REF] Banks | Estimation techniques for distributed parameter systems, volume 1 of Systems & Control: Foundations & Applications[END_REF]):

J(X 0 , θ 0 ) = N n=0 Γ n (Z n , X n ) 2 W -1 n + θ 0 -θ+ 0 2 (P θ 0 ) -1 + X 0 -X+ 0 2 (P X 0 ) -1 ,
with X n satisfying (2). In this expression, X0 and θ0 are given a priori values and

• W -1 ,
• (P X 0 ) -1 and • (P θ 0 ) -1 denote some matrix norms used to give a different weight to the different terms.

This minimization problem can be addressed by many methods that are classically divided in two groups: the variational and the sequential approaches. Variational approaches minimize this cost function by an optimization algorithm -usually a gradient-based computed from the solution of an adjoint model -and require numerous computations of the forward -and possibly the adjoint -problem, see e.g. [START_REF] Martin | Parameter identification for a one-dimensional blood flow model[END_REF]. Here, we follow a sequential, or filtering, approach which modifies the forward dynamics (2) with a correction term proportional to the observation error. More precisely, we adopt the Reduced-Order Unscented Kalman Filter (ROUKF) [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF], inspired from [START_REF] Julier | A new approach for filtering nonlinear systems[END_REF][START_REF] Julier | A new method for the nonlinear transformation of means and covariances in filters and estimators[END_REF]. It does not require any tangent operator and allows one to run the estimation with a computational cost of the same order of the forward problem since it is highly parallelizable. Moreover, it requires only superficial modifications to the existing solvers. Its main assumption consists of neglecting the uncertainty in X 0 . Hence P X 0 → ∞ and a reduced-order form of the Kalman filter can be obtained. Note that it is also possible to take into account the uncertainties in the initial condition in the same framework with an additive nudging based correction on the dynamics as introduced and analyzed in [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspective in medical imaging[END_REF][START_REF] Bertoglio | State observers of a vascular fluid-structure interaction model through measurements in the solid[END_REF] and illustrated in [START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF].

In fact, dealing with uncertainties in the initial condition is a difficult task, which is discussed in detail in [START_REF] Bertoglio | State observers of a vascular fluid-structure interaction model through measurements in the solid[END_REF] for the FSI case.

The ROUKF algorithm Consider the matrix of simplex sigma-points computed recursively as

[I ( * ) 1 ] = - 1 √ 2α , 1 √ 2α , α = 1 κ + 1 , and for 2 ≤ d ≤ κ [I ( * ) d ] =          0 [I ( * ) d-1 ] . . . 0 1 αd(d + 1) • • • 1 αd(d + 1) -d αd(d + 1)         
. [START_REF] Pham | Filtres de Kalman singuliers évolutifs pour l'assimilation de données en océanographie[END_REF][START_REF] Hoteit | A simplified reduced order Kalman filtering and application to altimetric data assimilation in Tropical Pacific[END_REF]. Then, for given values of X+ 0 , θ+ 0 and P θ 0 , perform • Initialization:

L θ 0 = I, L X 0 = 0, U 0 = (P θ 0 ) -1 , (3a)
with L θ 0 and L X 0 the initial sensitivities for the parameters and initial condition, respectively. Note that this assumes that the uncertainty at the initial time step is present only in the parameters and not in the initial condition, see [START_REF] Moireau | Reduced-order Unscented Kalman Filtering with application to parameter identification in large-dimensional systems[END_REF] for details. Then for n = 0, . . . , N -1 perform: • Propagation:

• Sampling:      C n = U -1 n (Cholesky factorization), X(i)+ n = X+ n + L X n C n I (i) κ , 1 ≤ i ≤ κ + 1, θ(i)+ n = θ+ n + L θ n C n I (i) κ , 1 ≤ i ≤ κ + 1, (3b) where
( X(i)- n+1 , θ(i)- n+1 ) = A n+1 ( X(i)+ n , θ(i)+ n ) , i = 1, . . . , κ + 1, Γ(i) n+1 = Γ(Z n+1 , X(i)- n+1 ) , i = 1, . . . , κ + 1, (3c)
which can be performed independently for each i, leading to the intrinsical parallel feature of this algorithm.

• Prediction:

                     X- n+1 = α κ+1 i=1 X(i)- n+1 , θ- n+1 = α κ+1 i=1 θ(i)- n+1 , Γ- n+1 = α κ+1 i=1 Γ(i)- n+1 . (3d)
• Correction:

                   L X n+1 = α[ X( * )- n+1 ][I ( * ) p ] , L θ n+1 = α[ θ( * )- n+1 ][I ( * ) p ] , L Γ n+1 = α[Γ ( * ) n+1 ][I ( * ) p ] , U n+1 = P α + (L Γ n+1 ) W -1 n+1 L Γ n+1 , X+ n+1 = X- n+1 -L X n+1 U -1 n+1 L Γ n+1 W -1 n+1 Γ- n+1 , θ+ n+1 = θ- n+1 -L θ n+1 U -1 n+1 L Γ n+1 W -1 n+1 Γ- n+1 . (3e)
Here we choose W -1 n+1 = γM Σ , with M Σ ∈ R m×m the L 2 -mass matrix on the fluidstructure interface and γ > 0.

Observation operator As a meaningful observation we rely here on the available segmentations of medical images usually given as a set of surfaces, here the fluid-structure interface at each acquisition time. Hence, we quantify an observation error by computing the (signed) distance field between the segmented surfaces and the deformed fluid-structure interface [START_REF] Moireau | Filtering for distributed mechanical systems using position measurements: Perspective in medical imaging[END_REF].

Image registration based segmentation

The data assimilation process requires estimates of the displacements of points on the vessel wall at intervals over the cardiac cycle, typically measured from segmentations of dynamic image data. In this work, the segmentation is performed by registration of one reference image to each frame (target image) of the dynamic images using the Sheffield Image Registration Toolkit (ShIRT). There follows a brief summary of the algorithm and implementation: for more details see [START_REF] Barber | Automatic segmentation of medical images using image registration: diagnostic and simulation applications[END_REF][START_REF] Barber | Efficient computational fluid dynamics mesh generation by image registration[END_REF].

Image registration algorithm

The registration process seeks to determine the displacement field that brings the reference image into alignment with the target image, assuming that the intensities of corresponding points on the two images are related by

f (x + u) = m(x) (4)
with f and m the intensities of the fixed and moving images, respectively, and x = [x, y, z] the cartesian spatial coordinates and u = [u, v, w] the displacement field.

If the objects in the images have the same topology and, crucially, the same ranges of intensity then it is possible in principle to find a displacement field u that satisfies Equation (4) exactly. Otherwise it is satisfied with minimum error described by a cost function. In fact small variation in image intensity ranges, including gradients across the image that are commonly encountered in some modalities, can be accommodated by adding an additional image dimension, s, that represents the intensity of the image. Thus f (x) and m(x) are then transformed to f ( x) and m( x), respectively, with x = [x, y, z, s], and f , m taken binary values, with the height of each column in the s direction equal to the intensity of the image at that point. The displacement u is redefined now including a term in the s direction, i.e., f

( x + û) = m( x) , û = [u, v, w, u s ]. (5)
For small values of the displacements, the relationship between the intensities of the two images at the common coordinate x is given by m

( x) = f ( x) + û • ∇ x f | x . (6)
Note that the registration problem can also be formulated in terms of the same displacements applied to the image m, i.e.

f ( x) = m( x) -û • (∇ x m) | x . (7)
Combining both Equations ( 6) and ( 7) we can define the following residual for the registration

R( x) def = m( x) -f ( x) - û 2 • ∇ x( m + f ) | x (8)
which in our experience improves the performance of the registration compared with considering the residuals only based on Equations ( 6) or [START_REF] Bertoglio | Forward and inverse problems in fluid-stucture interaction[END_REF].

The registration algorithm solves then a sequence of regularized linear least-squares problems based on a cost function that evaluates the discrete residual in Equation ( 8) after registration, namely, for p > 1 solve

argmin U p 1 2 (M -F ) -[T ]U p 2 + λ 2 [∆] p-1 q=1 U q 2
with U p the degrees-of-freedom of the trilinear interpolation of û, [∆] corresponds to the discrete Laplacian operator included as the regularisation term (with λ > 0 the corresponding weighting factor), and the matrix [T ] arises from the discretisation of the gradients in Equation ( 8).

The process is assumed to be converged when the change of U p over an iteration is less than a pre-determined value. At this stage the grid resolution is increased and the operation starts again, producing a converged solution at the next level of grid refinement. In principle this can continue until the grid on which the displacement is defined is equal to the voxel resolution of the images, but in practice there is little benefit in reducing the grid below 2x2x2 voxels.

Registration-based segmentation of the MR-data For the elastic tube phantom reported in Section 5, the reference image is determined by registering each frame in the time series back to the first frame and then averaging the result. The mesh points that will be passed to the data assimilation step are determined by registering an image of a hollow cylinder on which a mesh is defined to the reference image and then propagating the interpolated displacements from the registration of the time series images. The initial segmentation of the human aorta reported in Section 6 is performed by operation on a 3D contrast-enhanced MR image, using an extension of ShIRT that performs a series of registrations to grow the segmentation from a regional seed. This is preferred to direct segmentation of one frame of the dynamic image because the image quality is better. A triangular surface mesh is defined on the surface of the resulting lumen, and this mesh is mapped onto the dynamic image at one (diastolic) frame, again by registration, with trilinear interpolation of the registration field. This reference frame is registered to each frame of the dynamic image in turn, and the displacements of the nodes of the mesh are again determined by trilinear interpolation. Note that despite the non-physical assumptions in the registration procedure leading to potential errors in the displacements computed at this stage, the surfaces generated by the segmentation are accurate and therefore meaningful as an input in the data assimilation procedure presented in Section 3 (see also [START_REF] Moireau | Sequential identification of boundary support parameters in a fluid-structure vascular model using patient image data[END_REF]).

In vitro validation: a silicon rubber aortic phantom

Experimental setting The experimental setup is described in Figure 1, and consists of a MR-compatible emulator of the systemic circulation. The fluid pump corresponds to a linear servo actuator (ETB-32, Parker Hannifin) and piston assembly, located outside the 5 Gauss line coupled to the ventricle piston via a 2 m rigid boom. The piston ejects water through a tri-leaflet polyurethane valve, (Hemolab, Eindhoven, The Netherlands), and into a 450 mm length, 18.5 mm diameter and 1 mm thickness silicone tube, (TU Eindhoven, The Netherlands), with an elastic strap distally to make it stiffer there. Afterload is maintained by a mechanical Windkessel model at the distal end of the tube comprising fluid resistance, then compliance, then another resistance. After that the working fluid returns to the pump through a venous reservoir and an additional valve. At the centre of the aorta length, a rubber collar was prepared and applied to simulate a local coarctation. The collar was held in place by a spot of glue at its proximal end. Measured data Pressure measurements were taken at 13 locations between the venous and arterial compliance chambers with a measurement spacing of 20 mm using a pressure wire, (Certus, St. Jude Medical Systems, Uppsala, Sweden) at a temporal resolution of 1000 Hz.

The diastolic values for the pressure vary along the tube between 26 and 29 mmHg, and the systolic between 105 and 110 mmHg 0.3 s later.

Then, MR-acquisitions were performed on a 3T scanner, (Philips Achieva, Philips Healthcare, Best, The Netherlands). Geometry imaging was captured in a balanced turbo field echo (B-TFE) cine scan with a temporal resolution of 0.01 s and 13 slices, the same of the pressure measurements, were taken with an excitation slice thickness of 8 mm, and a pixel size of 1.8 × 2.5 mm. Image segmentation was then carried out using the methodology described in Section 4.

A non-destructive mechanical test was then performed consisting in placing the tube (without the strap) vertically, filling it with water, and then injecting an additional controlled water volume while measuring the pressure at the bottom, see Figure 2. A pressure dependent Young's modulus was then constructued using basic linear solid mechanics theory as

E(p) ≈ (2 -ν)V 3 2 δp √ L π δV h (9)
with ν the Poisson's ratio, h the tube thickness, L the tube length and δp and δV the measured volume and pressure increments, respectively. The FSI model The geometry used for the mesh construction (and also initial condition) was chosen as the segmented surfaces at time 0.19 s, where the pressure varies between 44 and 50 mmHg along the tube. The surface was then remeshed using 3-matic (Materialise, Leuven, Belgium), and the volume was closed with additional surfaces at the outlet and the inlet. Then, using Gmsh [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF], the surface was extruded outwards the cylinder in two tetrahedral layers with total thickness of 1 mm. Finally both solid's and fluid's volume meshes were constructed with 135840 and 255023 tetrahedra, respectively, resulting from setting the mesh size in the mesh generation processs between 1 mm and 1.5 mm. The boundary conditions for the fluid were directly taken from the pressure measurements at inlet and outlet. In the experiment, the real tube was radially and axially constrained at the physical boundaries. However, since the limits of the computational domain are away from the constrained boundaries, in the FSI-model we only fixed both inlet and outlet axially. Hence, in practice we neglected the radial stiffness of the parts of the tube not considered in the analysis. The physical parameters of Model (1) were chosen as: c Γ = 10 4 , k Γ = 0, µ = 0.035, ρ f = 1, ρ s = 1.2, ν = 0.45, η s = 0. Since the initial pressure field is not zero, the loads received by the solid are corrected by the initial one, so that the system starts at rest.

Estimation results

We subdivide the solid volume in 10 regions and estimate simultaneously the Young's modulus for each region. As it can be seen in Figure 3, the estimation algorithm succesfully detects the stiffer regions (first and second, which correspond to the actual placement of the strap), while in the regions where the strap is not presented the stiffness values are close to the one derived from the mechanical test. Specifically, the estimated stiffness in regions 3 and 4 match perfectly the experimental value, while the rest gradually increases when approaching the end of the tube. This could be explained by the fact that we did not include any radial stiffness in order to represent the neglected part of the tube in the computational model.

In vivo results: an aorta with mild coarctation

Clinical data The available clinical data corresponds to a 19 year old male subject with a mild (repaired) coarctation of the descending aorta, see Figure 4. A static 3D MR-image was obtained with a 10-seconds untriggered acquisition and a gadolinium contrast agent with voxel size is 1.36 × 1.36 × 1.8 mm and then segmented using with a registration based segment growing algorithm [START_REF] Barber | Registration based segment growing for vascular segmentation[END_REF] (see also Section 4). A 4D-SSFP sequence was acquired (25 time frames) along the cardiac cycle of 1.324 s were acquired with a voxel size of 1.29×1.29× 2.4 mm. Consequently, the 25 surfaces were segmented with the template mesh obtained from the 3D-Gd, see Figure 5. Four fixed-in-space slices of in-plane velocity phase contrast MRimages Also invasive pressures were measured at the ascending, descending, and abdominal aorta. The FSI model For the mesh construction, we proceeed concordingly as for the phantom, starting from the first segmented surface from the 4D-SSFP secquence, obtaining both solid and fluid volume meshes with 103406 and 138648 tetrahedra, respectively, resulting from setting the mesh size in the mesh generator between 0.5 mm and 1.5 mm, where an arterial wall thickness of 2 mm was assumed. The physical parameters were chosen as: µ = 0.035, ρ f = 1, ρ s = 1.2 and η s = 10 5 . The solid volume was divided into 5 regions (see Figure 6-left) where a Mooney-Rivlin constitutive model with parameters c j 1 = 3 • 10 5 • 2 θ j , c j 2 = 10 4 and κ j = 10 8 for each region j = 1, . . . , 5. The parameters θ j will be estimated later. The factor 3 • 10 5 was chosen so that the overall vessel deformation is of the order of the one observed in the images with θ j = 0. The solid boundary conditions parameters where chosen as: k Γ = 0 for the supra-aortic branches, and k Γ = 10 5 for the abdominal outlet, while c Γ = 10 4 everywhere on the outer surface. Inspired by [START_REF] Moireau | External tissue support and fluid-structure simulation in blood flows[END_REF], we extract a displacement field from the segmented surfaces and impose it at the model's inlet for taking the heart motion into account, see [START_REF] Bertoglio | Forward and inverse problems in fluid-stucture interaction[END_REF] for details.

A plane velocity profile in the inlet at the level of the ascending aorta was imposed, with amplitude obtained by the flow from the IP-velocity data (20 lit/min at 0.15 s). The Windkessel models were calibrated manually in order to approach the measured fluxes and pressures, see [START_REF] Bertoglio | Forward and inverse problems in fluid-stucture interaction[END_REF] for details. The load sent from fluid to structure was again corrected by the initial one.

Estimation results

The estimation results are presented in Figure 6. We see that the parameters adapt in time in order for the simulation to pursue the MR images. In fact, we present in Figure 7 some time frames of the MR-images, segmentations and models with the initial and calibrated parameters, showing that the model with estimated parameters is closer to the segmented surfaces. However, the results do not converge as in the in vitro validation. This is most probably due to modelling errors, since for example external tissue supports are not taken into account in this case. Another likely cause lies in the low signal-to-noise ratio (SNR) (see Figure 5), in particular in the coarctation region. Hence, the estimated stiffness parameters evolve and oscillate in time in order to compensate for these modelling and data errors. Despite this lack of convergence, the parameter estimation process shows reasonable robustness, as verified with multiple estimation cycles, see in Figure 8 where we relaunch the estimation starting from the previously estimated parameters -but keeping the same initial unitary covariance matrix. We remark, indeed, that the initial and final values get closer after each estimation run and the oscillations stabilize. Hence, we prove here that our method is of practical use for complex clinical data and simulation objectives, and represents a powerful tool to identify patient-specific parameters and to assess necessary improvements for the model or the data processing.

Conclusions

We performed an in vitro validation of the parameter estimation algorithm for coupled fluid-structure systems presented in [START_REF] Bertoglio | Sequential parameter estimation in fluidstructure problems. Application to hemodynamics[END_REF]. For this purpose, we considered a silicon rubber tube emulating the aorta, connected with a mechanical model of the main components of the cardiovascular system. We constructed a FSI computational model of the cylinder based on the data segmented from dynamic MR-images and we used pressure measurements for the fluid boundary conditions. We estimated the Young's modulus for a linear solid. We verified

  particles, and (•)(i) denotes the ith-column of the matrix [(•) ( * ) ].

Figure 1 :

 1 Figure 1: MRI compatible mock ventricle and aorta pulsatile phantom, where; 1, reciprocating ventricle piston; 2, ventricle and working fluid; 3, arterial valve; 4, suspension fluid; 5, outer cylindrical case; 6, venous compliance chamber; 7, silicone aorta; 8, venous channel; 9, venous resistance screw; 10, arterial resistance screw and arterial catheter access; 11, arterial compliance chamber; 12, venous valve.

Figure 2 :

 2 Figure 2: Pressure P i and volume V i , i = 1, . . . , 27, acquired from the mechanical test. The raw data is plot with + and the filtered data (4th order interpolation of Volume-Pressure data) with continues lines.

Figure 3 :

 3 Figure 3: Estimation results for the 10 regions, with γ = 10 and reference value obtained from the mechanical test and the formula (9) with p = 47 mmHg.

Figure 4 :

 4 Figure 4: Measurements available for the study case 10 of the euHeart database. Static segmentation (pink geometry); dynamic segmentation (green geometry); invasive pressure measurements at three locations: ascending (red), descending (green) and abdominal (blue) aorta; and In-Plane velocity measurements.

  (a) t = 0 s (b) t = 0.110 s (c) t = 0.221 s (d) t = 0.386 s

Figure 5 :

 5 Figure 5: Four sample time frames and the segmented surfaces (in magenta) for a representative slice of the 4D-SSFP sequence. Notice the changes of contrast of the images due to the changes in blood velocity.

Figure 6 :

 6 Figure 6: Estimation results for the constitutive parameters θ j , j = 1, . . . , 5, using the segmented surfaces. Left: geometry of the model's aorta and subdivisions of the solid domain. Right: mean estimated values with γ = 10 3 and P θ 0 = I.

  (a) t = 0.21 s (b) t = 0.53 s

Figure 7 :

 7 Figure 7: Comparison for some time instants between image, segmented surface (magenta), and FSI-model with initial (blue) and optimal (green) parameters. The view corresponds to a foot-head cut of the ascending aorta at the level of the pulmonary trunk.

Acknowledgements

The research leading to these results has received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement 224495 (euHeart project). We also thank sincerely to Cristina Staicu (U. Sheffield) for assistance in the data processing, and to Miguel A. Fernández and Dominique Chapelle (INRIA) for valuable discussions about fluid-structure interaction and data assimilation.

one of the estimated parameters from data coming from an independent non-destructive mechanical test. In addition, we show that we can straightforwardly use the algorithm in a real clinical case, namely for the estimation of the Mooney-Rivlin constitutive parameters in an aorta with repaired coarctation from the 4D-SSFP sequence, obtaining meaningful results for the estimated parameter values, in the sense that the model with estimated parameters follows closer the segmented surfaces. However, in order to fully exploit the potential of the parameter estimation process, we should also improve both modelling and data acquisition procedures, including the sliding-contact of the aorta with surrounding tissues, more realistic constitutive models and prestress, together with improvements in the quality of the acquired geometry data and fluid field measurements for a better pressure field estimation.

In conclusion, this work illustrates that the combination of automatic image segmentation tools, detailed fluid-structure models, efficient numerical algorithms for their simulation and robust data assimilation procedures provide a powerful toolchain for the in vivo analysis of arterial tissue properties.

Conflict of interests statement

The authors of this manuscript have no conflict of interest with the presented work.


Anne Benoit
email: [anne.benoit@ens-lyon.fr

Louis-Claude Canon
email: louis-claude.canon@univ-fcomte.fr

Loris Marchal
email: loris.marchal]@ens-lyon.fr

Non-clairvoyant reduction algorithms for heterogeneous platforms

come

Introduction

Reduction is one of the most common collective operations, together with the broadcast operation. Contrarily to a broadcast, it consists in gathering and summarizing information scattered at different locations. A classical example is when one wants to compute the sum of (integer) values distributed over a network: each node owns a single value and can communicate with other nodes and perform additions to compute partial sums. The goal is to compute the sum of all values. Reductions have been used in distributed programs for years, and standards such as MPI usually include a "reduce" function together with other collective communications (see [START_REF] Pjesivac-Grbovic | Performance analysis of MPI collective operations[END_REF][START_REF] Thakur | Optimization of Collective communication operations in MPICH[END_REF] for experimental comparisons). Many algorithms have been introduced to optimize this operation on various platforms, with homogeneous [START_REF] Rabenseifner | Optimization of Collective Reduction Operations[END_REF] or heterogeneous communication costs [START_REF] Liu | An Approximation Algorithm and Dynamic Programming for Reduction in Heterogeneous Environments[END_REF][START_REF] Legrand | Optimizing the steady-state throughput of scatter and reduce operations on heterogeneous platforms[END_REF]. Recently, this operation has received more attention due to the success of the MapReduce framework [START_REF] Dean | Mapreduce: Simplified data processing on large clusters[END_REF][START_REF] Zaharia | Improving mapreduce performance in heterogeneous environments[END_REF], which has been popularized by Google. The idea of MapReduce is to break large workloads into small tasks that run in parallel on multiple machines, and this framework scales easily to very large clusters of inexpensive commodity computers. We review similar algorithms and use cases in the related work section in the companion research report [START_REF] Benoit | Non-clairvoyant reduction algorithms for heterogeneous platforms[END_REF].

Our objective in this paper is to compare the performance of various algorithms for the reduce operations in a non-clairvoyant setting, i.e., when the algorithms are oblivious to the communication and computation costs (the time required to communicate or compute). This models well the fact that communication times cannot usually be perfectly predicted, and may vary significantly over time. We would like to assess how classical static algorithms perform in such settings, and to quantify the advantage of dynamic algorithms (if any). We use various techniques and models, ranging from worst-case analysis to probabilistic methods such as Markov chains.

The design and the analysis of algorithms in a dynamic context has already received some attention. The closest related work is probably [START_REF] Canon | Comparative Evaluation of the Robustness of DAG Scheduling Heuristics[END_REF], in which the authors study the robustness of several task-graph scheduling heuristics for building static schedules. The schedules are built with deterministic costs and the performance is measured using random costs. [START_REF] Benoit | Computing the throughput of probabilistic and replicated streaming applications[END_REF] studies the problem of computing the average performance of a given class of applications (streaming applications) in a probabilistic environment. With dynamic environments comes the need for robustness to guarantee that a given schedule will behave well in a disturbed environment. Among others, [START_REF] Canon | Evaluation and optimization of the robustness of dag schedules in heterogeneous environments[END_REF] studies and compares different robustness metrics for makespan/reliability optimization on task-graph scheduling.

The rest of the paper is organized as follows. Section 2 describes four algorithms and shows that they are approximation algorithms with identical computation costs. In Section 3, we provide more involved probabilistic analysis of their expected performance. Section 4 presents simulated executions of the previous algorithms and compares their respective performance. Finally, we conclude and discuss future research directions in Section 5.

Model and Algorithms

We consider a set of n processors (or nodes) P 0 , . . . , P n-1 . Each processor P i owns a value v i . We consider an associative operation ⊕. Our goal is to compute

the value v = v 0 ⊕ v 1 ⊕ • • • ⊕ v n-1
as fast as possible, i.e., to minimize the total execution time to compute the reduction. We do not enforce a particular location for the result: at the end of the reduction, it may be present on any node.

There are two versions of the problem, depending on whether the ⊕ operation is commutative or not. For example, when dealing with numbers, the reduction operation (sum, product, etc.) is usually commutative while some operations on matrices (such as the product) are not. The algorithms proposed and studied below deal with both versions of the problem.

We denote by d i,j the time needed to send one value from processor P i to processor P j . A value may be an initial value or a partial result. When a processor P i receives a value from another processor, it immediately computes the reduction with its current value. We assume that each processor can receive at most one result at a time. The communication costs are heterogeneous, that is we may well have different communication costs depending on the receiver (d i,j = d i,j ′ ), on the sender (d i,j = d i ′ ,j ) and non-symmetric costs (d i,j = d j,i ).

Even though these costs are fixed, we consider non-clairvoyant algorithms that make decisions without any knowledge of these costs.

The computation time of the atomic reduction on processor P i is denoted by c i . In the case of a non-commutative operation, we ensure that a processor sends its value only to a processor that is able to perform a reduction with its own value. Formally, assume that at a given time, a processor owns a value that is the reduction of v i ⊕ • • • ⊕ v j , which we denote by [v i , v j ]; the processor may only send this value to a processor owning a value

[v k , v i-1 ] or [v j+1 , v k ],
which is called a neighbor value in the following.

Reduction Algorithms. During a reduction operation, a processor sends its value at most once, but may receive several values. It computes a partial reduction each time it receives a value. Thus, the communication graph of a reduction is a tree (see Figure 1): the vertices of the tree are the processors and its edges are the communications of values (initial or partially reduced values). In the example, P 0 receives the initial value from P 1 , and then a partially reduced value from P 2 . In the following, we sometimes identify a reduction algorithm with the tree it produces.

We now present the four algorithms that are studied in this paper. The first two algorithms are static algorithms, i.e., the tree is built before the actual reduction. Thus, they may be applied for commutative or non-commutative reductions. The last two algorithms are dynamic: the tree is built at run-time and depends on the durations of the operations.

The first algorithm, called Binomial-stat, is organized with ⌈log 2 n⌉ rounds. Each round consists in reducing a pair of processors that own a temporary or initial data using a communication and a computation. During round k = 1, . . . , ⌈log 2 n⌉, each processor i2 k +2 k-1 (i = 0, . . . , 2 ⌈log 2 n⌉-k -1) sends its value to processor i2 k , which reduces it with its own value: at most 2 ⌈log 2 n⌉-k+1 processors are involved in round k. Note that rounds are not synchronized throughout the platform: each communication starts as soon as the involved processors are available and have terminated the previous round. We can notice that the communication graph induced by this strategy is a binomial tree [6, Chapter 19], hence the name of the algorithm. This strategy is illustrated on Figure 2(a).

The second algorithm, called Fibonacci-stat, is constructed in a way similar to Fibonacci numbers. The schedule constructed for order k, denoted by FS k (k > 0) first consists in two smaller order schedules FS k-1 and FS k-2 put in parallel. Then, during the last computation of FS k-1 , the root of FS k-2 (that is, the processor that owns its final value) sends its value to the root of FS k-1 , which then computes the last reduction. A schedule of order -1 or 0 contains a single processor and no operation. This process is illustrated on Figure 2(b).

Obviously, the number of processors involved in such a schedule of order k is F k+2 , the (k + 2)th Fibonacci number. When used with another number n of processors, we compute the smallest order k such that F k+2 ≥ n and use only the operations corresponding to the first n processors in the schedule of order k.

The previous two schedules were proposed in [START_REF] Canon | Scheduling Associative Reductions with Homogeneous Costs when Overlapping Communications and Computations[END_REF], where their optimality is proved for special homogeneous cases: Binomial-stat is optimal both when the computations are negligible in front of communications (c i = 0 and d i,j = d) and when the communications are negligible in front of computations (c i = c and d i,j = 0). Fibonacci-stat is optimal when computations and communications are equivalent (c i = c = d i,j = d). In the non-commutative case, both algorithms build a tree such that only neighboring partial values are reduced. In the commutative case, any permutation of processors can be chosen.

Then, we move to the design of dynamic reduction algorithms, i.e., algorithms that take communication decisions at runtime. The first dynamic algorithm, called Tree-dyn, is a simple greedy algorithm. It keeps a slot (initially empty), and when a processor is idle, it looks into the slot. If the slot is empty, the processor adds its index in the slot, otherwise it empties the slot and starts a reduction with the processor that was in the slot (i.e., it sends its value to the processor that was in the slot, and the latter then computes the reduced value). It means that a reduction is started as soon as two processors are available. Since in the obtained reduction tree, any two processors may be paired by a communication, this can only be applied to commutative reductions.

Finally, Non-Commut-Tree-dyn, is an adaptation of the previous dynamic algorithm to non-commutative reductions. In this algorithm, when a processor 2. Schedules for Binomial-stat of order 3 and Fibonacci-stat of order 4, both using 8 processors. For Fibonacci-stat, the two schedules of order 2 and 3 used in the recursive construction are highlighted in green and red.

is idle, it looks for another idle processor with a neighbor value (as described above). Now, we keep an array of idle processors rather than a single slot. If there is an idle neighbor processor, a communication is started between them, otherwise the processor waits for another processor to become idle.

Worst-Case Analysis. We analyze the commutative algorithms in the worst case, and we provide some approximation ratios, focusing on communication times. A λ-approximation algorithm is a polynomial-time algorithm that returns a solution whose execution time is at most λ times the optimal execution time. We let ∆ = D d , where d = min i,j d i,j and D = max i,j d i,j . We consider that c i = c, i.e., all computation costs are identical. Results are summarized in Table 1. Due to lack of space, the complete set of theorems and proofs is available in the companion research report [START_REF] Benoit | Non-clairvoyant reduction algorithms for heterogeneous platforms[END_REF]. We discuss a subset of the results below.

Theorem 1. Without computation cost, Binomial-stat and Tree-dyn are ∆approximation algorithms, and this ratio can be achieved.

The proof is in [START_REF] Benoit | Non-clairvoyant reduction algorithms for heterogeneous platforms[END_REF], and we exhibit here an instance on which the ratio is achieved. Let d i,j = d for 1 ≤ i < j ≤ n and d i,j = D for all remaining 1 ≤ j < i ≤ n. With both Binomial-stat and Tree-dyn, we consider that any processor P i sends its element to a processor P j such that i > j, which takes a time D⌈log 2 (n)⌉. The optimal solution, however, consists in avoiding any communication of size D (with total time of d⌈log 2 (n)⌉).

Theorem 2. Without computation cost, Fibonacci-stat is a (∆/ log 2 ϕ+∆/⌈log 2 n⌉)-approximation algorithm, where ϕ = 1+ √ 5 2
is the golden ratio (1/ log 2 ϕ ≈ 1.44).

Proof. Since computation costs are negligible, the makespan of Fibonacci-stat schedule in the worst case is kD, with k the order of the Fibonacci schedule [START_REF] Canon | Scheduling Associative Reductions with Homogeneous Costs when Overlapping Communications and Computations[END_REF]. We know that n > F k+1 and, by definition of the Fibonacci numbers, we have

F k+1 = 1 √ 5 (ϕ k+1 -(1 -ϕ) k+1 ). Since -1 < 1 -ϕ < 0, it follows that n > F k+1 ≥ 1 √ 5 (ϕ k+1 -(1 -ϕ)
2 ) (as soon as k ≥ 1). We therefore have ϕ k+1 ≤ √ 5n + (1 -ϕ) 2 , and

(k + 1) log 2 ϕ ≤ log 2 n + log 2 √ 5 + (1-ϕ) 2 n . Thus, k ≤ log 2 n log 2 ϕ + log 2 √ 5 + (1-ϕ) 2 n log 2 ϕ -1 ≤1 when n≥1 c = 0 any c c = d Binomial-stat ∆ (Th. 1) ∆ + 1 (Th. 3) (∆ + 1)(1 + 1 log 2 n ) log 2 ϕ (Th. 4) Tree-dyn ∆ (Th. 1) ∆ + 1 (Th. 3) (∆ + 1)(1 + 1 log 2 n ) log 2 ϕ (Th. 4) Fibonacci-stat ∆ log 2 ϕ + ∆ ⌈log 2 n⌉ (Th. 2) ∆ log 2 ϕ + 2∆ ⌈log 2 n⌉ (Th. 5) ∆ (Th. 5)
Table 1. Approximation ratios for commutative algorithms. Theorem numbers in parenthesis refer to the theorems in [START_REF] Benoit | Non-clairvoyant reduction algorithms for heterogeneous platforms[END_REF].

Thus, k ≤ log 2 n log 2 ϕ + 1 ≤ ⌈log 2 n⌉ log 2 ϕ + 1. Recall that the lower bound is d⌈log 2 n⌉, hence the approximation ratio of ∆ log 2 ϕ + ∆ ⌈log 2 n⌉ . ⊓ ⊔

Markov Chain Analysis

In this section, we assume that communication and computation costs are exponentially distributed (i.e., each d i,j for 1 ≤ i, j ≤ n follows an exponential law with rate λ d and each c i for 1 ≤ i ≤ n follows an exponential law with rate λ c ). With this model, both dynamic approaches, Tree-dyn and Non-Commut-Tree-dyn, may be analysed using memoryless stochastic processes. Intuitively, each state of those processes is characterized by the number of concurrent communications and computations. A state in which there is neither communication nor computation is an absorbing state that corresponds to the termination of a reduction algorithm. In the initial state, there are ⌊ n 2 ⌋ concurrent communications (and one idle machine that is ready to send its value if n is odd). The completion time of an algorithm is then the time to reach the final state. To determine this duration, we use the first-step analysis.

Formally, let P be the transition rate matrix of a continuous-time Markov chain and s ∈ S be each state. Let φ(s) be a function on S taking real values (it will be the expected duration spent in state s or its variance). Let w i be the sum of the application of φ to each state taken by the process until the final state is reached starting from state s i . Then, w i is determined using the first-state analysis:

w i = φ(s i ) if s i is an absorbing state φ(s i ) + n j=1 p i,j w j otherwise
We apply this analysis to determine the expected duration and the variance of Tree-dyn with λ c = 0. For clarity, n is assumed to be even (the following analysis can be performed to the cases where n is odd by adapting the initial state).

In this case, each state s i,j is characterized by the number of concurrent communications i and whether the slot containing a ready processor is empty (j = 0) or not (j = 1). The initial state is s n 2 ,0 and the final state is s 0,1 . There are two kinds of transitions: from state s i,0 to state s i-1,1 (a communication terminates and the intermediate result of the local reduction is ready to be sent to the next available processor) and from state s i,1 to state s i,0 (a communication terminates and a new one is initiated with the available processor identified by the slot) for 1 ≤ i ≤ n 2 . In both cases, the rate of the transition is determined by the number of concurrent communications, that is iλ d .

In order to determine the expected completion time of Tree-dyn (noted C Tree-dyn ), we define φ(s i,j ) as 1 iλ d , the expected time spent is state s i,j (zero for state s 0,1 ). Therefore,

C Tree-dyn = w n 2 ,0 = n 2 i=1 φ(si,0)+φ(si-1,1) = n 2 i=1 1 iλ d + n 2 -1 i=1 1 iλ d = 1 λ d 2H n 2 -1 + 2 n ,
where H(n) is the nth harmonic number.

Similarly, we compute the variance of the completion time (noted V Tree-dyn ) by defining φ(s i,j ) as zero for state s 0,1 and 1 i 2 λ 2 d otherwise:

V Tree-dyn = 1 λ 2 d   2 n 2 -1 i=1 1 i 2 + 4 n 2   .

Simulation Results

In this section, we consider that the d ij and c i costs are distributed according to a gamma distribution, which is a generalization of exponential and Erlang distributions. This distribution has been advocated for modeling job runtimes [START_REF] Lublin | The workload on parallel supercomputers: modeling the characteristics of rigid jobs[END_REF][START_REF] Feitelson | Workload modeling for computer systems performance evaluation[END_REF]. This distribution is positive and it is possible to specify its expected value (µ d or µ c ) and standard deviation (σ d or σ c ) by adjusting its parameters. Using two distributions with two parameters each is a good compromise between a simulation design with a reasonable size and a general model. Additionally, further simulations with Bernoulli distributions concurred with the following observations, suggesting that our conclusions are not strongly sensitive to the distribution choice. Each simulation was performed with an ad-hoc simulator on a desktop computer.

Cost Dispersion Effect. In this first simulation, we are interested in characterizing how the dispersion of the communication costs affects the performance of all methods. In order to simplify this study, no computation cost is considered (µ c = 0). The dispersion is defined through the coefficient of variation (CV), which is defined as the ratio of the standard deviation over the expected value (this latter is set to 1). The number of processors is n = 64 and the time taken by each method is measured over 1 000 000 Monte Carlo (MC) simulations (i.e., simulations have been performed with 1 000 000 distinct seeds). On a global level, Figure 3 shows the expected performance with distinct CVs. When the heterogeneity is noticeable (CV greater than 1), the performance decreases significantly. In those cases, schedules are mostly affected by a few extreme costs whose importance depends on the CV. Additionally, the variability in the schedule durations is also impacted by the CV (i.e., two successive executions with the same settings may lead to radically different performance depending on the schedule).

Several observations can be made relatively to each method. As expected, Binomial-stat is similar to Tree-dyn for CVs lower than 10%. In this case, the improvement offered by Tree-dyn may not outweigh the advantage of following a static plan in terms of synchronization. For CVs greater than 1, both static approaches perform equally with a similar dispersion. For all CVs, Tree-dyn has the best expected performance while Fibonacci-stat has the worst, and Non-Commut-Tree-dyn has the second best expected performance when the CV is greater than 30%. Finally, when the CV is close to 10, all methods are equivalent as a single communication with a large cost may impact the entire schedule duration. In terms of robustness, we can see that Fibonacci-stat and Non-Commut-Tree-dyn are the two best methods for absorbing variations as their expected durations remains stable longer (until the CV reaches 30%). This is due the presence of idleness in their schedules that can be used when required.

Non-Negligible Computation. When considering nonzero computation costs, we reduce the number of parameters by applying the same CV to the computation and to the communication costs (i.e., σc µc = σ d µ d ). As Fibonacci-stat is designed for overlapping computations and communications, we characterize the cases when this approach outperforms Tree-dyn.

Figure 4(a) shows the improvement of Tree-dyn over Fibonacci-stat when varying the CV and the ratio µc µ d (the overlapping degree between computations and communications). The contour line with value 1 delimits the area for which Fibonacci-stat is better than Tree-dyn on average. This occurs when the computation cost is greater than around half the communication cost and when the variability is limited. When the computation costs are low ( µc µ d = 0.1), the ratio evolution is is consistent with the previous observations. Figure 4(a) is horizontally symmetrical as any case such that µc µ d > 1 is equivalent to the situation where the communication and the computation costs are swapped (and for which µc µ d < 1). These costs can be exchanged because a communication is always followed by a reduction operation.

Non-Commutative Operation. Finally, we assess the performance of Non-Commut-Tree-dyn by comparing it to all other methods that support a noncommutative operation when varying the dispersion and the overlapping degree as in the previous study. shows the method with the best average performance when varying the CV and the ratio µc µ d . We see that Non-Commut-Tree-dyn has the best performance when the cost dispersion is large. Additionally, the transition from Binomial-stat to Fibonacci-stat is when the computation cost reaches half the communication cost (as on Figure 4(a)). With low computation costs ( µc µ d = 0.1), the results are also consistent with Figure 3.

Conclusion

In this paper, we have studied the problem of performing a non-clairvoyant reduction on a distributed heterogeneous platform. Specifically, we have compared the performance of traditional static algorithms, which build an optimized reduction tree beforehand, against dynamic algorithms, which organize the reduction at runtime. Our study includes both commutative and non-commutative reductions. We have first proposed approximation ratios for all commutative algorithms using a worst-case analysis. Then, we have proposed a Markov chain analysis for dynamic algorithms. Finally, we have evaluated all algorithms through extensive simulations to show when dynamic algorithms become more interesting than static ones. We have outlined that dynamic algorithms generally achieve better makespan, except when the heterogeneity is limited and for specific communication costs (no communication cost for Binomial-stat, communication costs equivalent to computation costs for Fibonacci-stat). The worst-case analysis has also confirmed this last observation.

As future work, we plan to investigate more complex communication models, such as specific network topologies. It would also be interesting to design a better

Figure 1 .

 1 Figure 1. Schedule and communication graph for reducing four values. Blue arrows represent communications while red springs stand for computations.

Figure

  Figure2. Schedules for Binomial-stat of order 3 and Fibonacci-stat of order 4, both using 8 processors. For Fibonacci-stat, the two schedules of order 2 and 3 used in the recursive construction are highlighted in green and red.

Figure 3

 3 Figure 3. Average schedule length for each method over 1 000 000 MC simulations with n = 64, µ d = 1, µc = 0 and varying coefficients of variation for the communication costs. The lower part of the ribbons corresponds to the 10% quantile while the upper part corresponds to the 90% quantile for each method.

  -Tree-dyn Effect of the CV on the performances

  with varying overlapping and variability (b) Non-commutative operation

Figure 4 .

 4 Figure 4. Ratio of the average performance of Fibonacci-stat and Tree-dyn (a) and method with the best average performance (b) over 1 000 MC simulations for each square with n = 64, µ d = 1, σc µc = σ d µ d , varying coefficients of variation for the costs and varying µc µ d .

Figure 4 (

 4 Figure 4(b) shows the method with the best average performance when varying the CV and the ratio µc µ d . We see that Non-Commut-Tree-dyn has the best performance when the cost dispersion is large. Additionally, the transition from Binomial-stat to Fibonacci-stat is when the computation cost reaches half the communication cost (as on Figure4(a)). With low computation costs ( µc µ d = 0.1), the results are also consistent with Figure3.

dynamic algorithm for non-commutative reductions, which avoids the situation when many processors are idle but cannot initiate a communication since no neighboring processors are free.


Stability analysis for systems with saturation and backlash in the loop

Sophie Tarbouriech, Isabelle Queinnec, Christophe Prieur

Abstract-This paper deals with the stability analysis problem for linear systems with saturation and backlash in the loop. The resulting system controlled by a static output feedback is a dynamical model with nested backlash and saturation operators. Uniform ultimate boundedness stability is tackled in a regional (local) or global context depending on the stability property of the open-loop system. Suitable regions of the state space in which the closed-loop trajectories can be captured are characterized, together with estimates of the basin of attraction of such regions. Convex optimization problems are proposed in order to give a constructive solution (in the sense of theoretical conditions associated to numerical procedures).

I. INTRODUCTION

Physical control systems are often subject to magnitude limitation in the input and/or the output. This type of nonlinearity may drastically reduce the performance of the closed-loop system or even lead to instability. Therefore, for these reasons the stability analysis or stabilization problems of control systems with saturation in the input attracted research efforts for several decades (see, e.g., [START_REF]Actuator saturation control[END_REF], [START_REF] Hu | Control systems with actuator stauration: Analysis and design[END_REF], [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF]). Another common nonlinearity that limits performance for many applications is the backlash. These operators are nonlinear elements with memory. They are usual phenomena in many physical systems, such as electrical inductors, piezoactuators, gear trains and mechanical friction systems (see e.g. the survey [START_REF] Macki | Mathematical models for hysteresis[END_REF] on such memory elements and other operators as Krasnolsel'skii-Pokrovskii hysteresis). They still constitute an open source of both theoretical and practical problems [START_REF] Nordin | Controlling mechanical systems with backlash: a survey[END_REF], [START_REF] Barreiro | Input-output stability of systems with backlash[END_REF].

Few papers have considered both nonlinearities in view of stability analysis or control design purpose. One can cite [START_REF] Del | Controller structure for plants with combined stauration and deadzone/backlash[END_REF] in which the authors study the presence of saturation and backlash in series at the input of the system. The technique proposes to inverse the backlash in order to cancel its undesirable effects. Then the resulting system corresponds to a system with a new saturation for which the level is depending on the magnitude constraint of the previous saturation and on the parameters of the backlash. The current paper has chosen another route and then does not search for inverting the backlash element. Stability analysis results are studied from the use of Lyapunov theory arguments. More precisely, generalized sector conditions using properties of the saturation and backlash operators, associated to Lyapunov arguments allow us to analyze the stability in a regional (local) or global S. Tarbouriech and I. Queinnec are with CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France and Univ de Toulouse, LAAS, F-31400 Toulouse, France. tarbour@laas.fr, queinnec@laas.fr C. Prieur is with Gipsa-lab, Department of Automatic Control, Grenoble Campus, 11 rue des Mathématiques, BP 46, 38402 Saint Martin d'Hères, France christophe.prieur@gipsa-lab.fr context. The regional result (see Theorem 1 below), presents sufficient conditions to prove that a compact set (inner set) is a finite-time attractor, for all admissible initial conditions belonging to another compact set (outer set). Then, the second contribution of the paper (see Theorem 2) adapts the technique of Theorem 1 to the global case. In this case, the outer set becomes all the state space, provided that the openloop system is not strictly exponentially unstable. Furthermore, at the opposite of conditions developed in [START_REF] Tarbouriech | Global stability for systems with nested backlash and saturation operators[END_REF], we do not need to use the time-derivative version of the system. Moreover, while [START_REF] Jayawardhana | The circle criterion and input-to-state stability[END_REF] addresses the characterization of the global finite-time attractor (no saturation being considered), our approach uses the knowledge of the Lyapunov function and, in this sense, can be considered as bringing a new brick to complement the method proposed in [START_REF] Jayawardhana | The circle criterion and input-to-state stability[END_REF]. Both in the local and global cases, the conditions are given in terms of bilinear matrix inequalities and associated algorithms based on LMI constraints are proposed to numerically implement the theoretical conditions.

Notation. For two vectors x, y of ℜ n , the notation x y means that x (i)y (i) ≥ 0, ∀i = 1, . . . , n. 1 and 0 denote the identity matrix and the null matrix of appropriate dimensions, respectively.

x ∈ ℜ n + means that x 0. The Euclidian norm is denoted • . A ′ and trace(A) denote the transpose and the trace of A, respectively. For two symmetric matrices, A and B, A > B means that A -B is positive definite. In partitioned symmetric matrices, the symbol ⋆ stands for symmetric blocks. λmax(A) (respectively, λmin(A)) denotes the maximal (respectively, minimal) eigenvalue of the matrix A.

II. PROBLEM FORMULATION

The class of systems under consideration consists of systems with saturation and backlash in the loop. The system is described by:

ẋ = Ax + Bu y = Cx + DΦ[w] w = Ex u = sat(Ky) (1)
where x ∈ ℜ n is the state, u ∈ ℜ m is the input of the plant, y ∈ ℜ p is the measured output of the plant and w ∈ ℜ q is the input of the nonlinearity Φ. A, B, C, D, E and K are known matrices of appropriate dimensions. sat denotes the classical vector-valued saturation function defined as ∀i = 1, ..., m: (sat(z)) (i) = sat(z (i) ) = sign(z (i) ) min(u 0(i) , |z (i) |) with u 0(i) > 0 the ith level of the saturation. Φ is a componentwise backlash operator (see, for example, [START_REF] Macki | Mathematical models for hysteresis[END_REF], [START_REF] Paré | A KYP lemma and invariance principle for systems with multiple hysteresis non-linearities[END_REF], [START_REF] Taware | Control of sandwich nonlinear systems[END_REF], [START_REF] Brokate | Hysteresis and phase transitions[END_REF]). We denote the set of continuous, piecewise differentiable functions f :

[0, +∞) → ℜ q by C 1 pw ([0, +∞); ℜ q ),
that is the set of continuous functions w being, for some unbounded sequence (t j ) ∞ j=0 in [0, +∞) with t 0 = 0, continuously differentiable on (t j-1 , t j ) for all j ∈ N. Given the vector ρ in ℜ q + and L = diag(ℓ (i) ), with positive values ℓ (i) , i = 1, . . . , q, the operator Φ is defined as follows, for all f ∈ C 1 pw ([0, +∞); ℜ q ), for all j ∈ N, for all t ∈ (t j-1 , t j ) and for all i ∈ {1, . . . , q}:

( ˙ Φ[f ](t)) (i) =            ℓ (i) ḟ(i) (t) if ḟ(i) (t) ≥ 0 and (Φ[f ](t)) (i) = ℓ (i) (f (i) (t) -ρ (i) ) ℓ (i) ḟ(i) (t) if ḟ(i) (t) ≤ 0 and (Φ[f ](t)) (i) = ℓ (i) (f (i) (t) + ρ (i) ) 0 otherwise
(2) where 0 = t 0 < t 1 < . . . is a partition of [0, +∞) such that f is continuously differentiable on each of the intervals (t j-1 , t j ), j ∈ N. Thus, Φ is a time-invariant nonlinearity with slope restriction, as in [START_REF] Park | The asymptotic stability of nonlinear (Lur'e) systems with multiple slope restrictions[END_REF]. Note however that it is a memory-based operator, since to compute it, we need to have information about the past values of w (this is not the case in [START_REF] Park | The asymptotic stability of nonlinear (Lur'e) systems with multiple slope restrictions[END_REF]).

Throughout the paper, we define the set of admissible initial conditions (w(t = 0)), from which we want to guarantee the stability, as follows

L(w(0) + ρ) Φ[w](0) L(w(0) -ρ) (3)
According to [START_REF] Nordin | Controlling mechanical systems with backlash: a survey[END_REF], [START_REF] Corradini | Robust stabilization of nonlinear uncertain plants with backlash or dead zone in the actuator[END_REF], that means that the nonlinearity Φ is active. Then, with (2), one gets

L(w(t) + ρ) Φ[w](t) L(w(t) -ρ), ∀t ≥ 0 (4)
Ψ = Φ -Lw = Φ -LEx (6)
Hence, the closed-loop system (1) reads:

ẋ = (A + BK(C + DLE))x + Bϕ + BKDΨ w = Ex (7)
It is important to note that both nonlinearities are nested since, by definition of y, ϕ depends on Ψ:

ϕ = sat(K 0 x + KDΨ) -(K 0 x + KDΨ) (8)
with

K 0 = K(C + DLE).
The presence of the backlash operator Φ may induce the existence of multiple equilibrium points or a limit cycle around the origin. Furthermore, in a neighborhood of the origin, system [START_REF] Barreiro | Input-output stability of systems with backlash[END_REF] operates in open loop. In parallel, the presence of the saturation asks for characterizing suitable regions of the state space in which the stability of the closedloop saturated system can be ensured [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF]. The regional (local) uniform ultimate boundedness stability (see [9, Section 4.8]) of the system is then carried out. Hence, the proposed approach aims at characterizing two sets such that the closedloop trajectories initialized in the outer set are ultimately bounded in the inner set. It is important to emphasize that the proposed technique does not require the open-loop system to be stable. Nevertheless, depending on the open-loop stability, the global uniform boundedness stability context is also carried out. In this case, the outer set corresponds to the whole state space. This is stated in the following problems:

Problem 1 (Regional Case): Considering A 0 = A + BK(C +DLE) = A+BK 0 being Hurwitz, characterize the regions S 1 and S 0 of the state space such that the closed-loop trajectories of the system (7) remain confined in S 1 and are uniformly ultimately bounded in the set S 0 , when initialized as in [START_REF] Corradini | Robust stabilization of nonlinear uncertain plants with backlash or dead zone in the actuator[END_REF].

Imposing A 0 = A + BK 0 Hurwitz is classical since it corresponds to system (1) when the saturation and the backlash Φ are neglected. That corresponds to consider u = Ky and to replace Φ by Lw (which corresponds to take ϕ = 0 and Ψ = 0 in system [START_REF] Jayawardhana | The circle criterion and input-to-state stability[END_REF]).

Problem 2 (Global Case): Assume that A is Hurwitz. Considering A 0 = A + BK(C + DLE) = A + BK 0 being Hurwitz, characterize the region S 0 of the state space in which the solutions of system (1) are uniformly ultimately bounded when initialized as in (3).

III. MAIN RESULTS

A. Preliminary results

From Lemma 3.1 in [START_REF] Tarbouriech | Stability analysis for linear systems with input backlash through sufficient LMI conditions[END_REF] and (4), we can formulate the following properties with respect to Ψ.

Lemma 1: For any diagonal positive definite matrices

N 1 , N 2 , N 3 in ℜ m×m , with N 3 ≥ 1, we have, for all w ∈ C 1 pw ([0, +∞); ℜ q ), for all t ∈ (t i-1 , t i ) ( Ψ + L ẇ) ′ N 1 Ψ ≤ 0, (9)
( Ψ + L ẇ) ′ N 2 ( Ψ + (1 -N 3 )L ẇ) ≤ 0, ( 10
) -Lρ Ψ Lρ (11)
where 0 = t 0 < t 1 < . . . is a partition of [0, +∞) such that w is continuously differentiable on each of the intervals (t j-1 , t j ), j ∈ N.

B. Stability analysis

The following result relative to Problem 1 can be given. Theorem 1:

Consider A 0 = A + BK(C + DLE) = A + BK 0 being Hurwitz.
If there exist a symmetric positive definite matrix P ∈ ℜ n×n , four diagonal positive definite matrices N 1 ∈ ℜ q×q , N 2 ∈ ℜ q×q , T 3 ∈ ℜ q×q and T ∈ ℜ m×m , a matrix G ∈ ℜ m×n , two positive scalars η, τ 1 satisfying the following conditions

M 1 < 0 (12) ρ ′ LT 3 Lρ -τ 1 ≤ 0 (14)
P G ′ (i) ⋆ ηu 2 0(i) ≥ 0, i = 1, ...., m (15)
η < 1

with M 1 defined in [START_REF] Park | The asymptotic stability of nonlinear (Lur'e) systems with multiple slope restrictions[END_REF] (see at the top of the next page), then, for any initial admissible conditions (x(0), Ψ(0)), with

M 1 =     A ′ 0 P + P A 0 + τ 1 P P B -K ′ 0 T -G ′ T P BKD -A ′ 0 E ′ LN 1 -A ′ 0 E ′ LN 2 ⋆ -2T -B ′ E ′ LN 1 -T KD -B ′ E ′ LN 2 ⋆ ⋆ -T 3 -N 1 LEBKD -D ′ K ′ B ′ E ′ LN 1 -N 1 -D ′ K ′ B ′ E ′ LN 2 ⋆ ⋆ ⋆ -2N 2     (13)
x(0) ∈ S 1 , the resulting trajectories of the closed-loop system (1) are uniformly ultimately bounded in the set S 0 , where the sets S 1 and S 0 are defined as follows:

S 1 = {x ∈ ℜ n ; x ′ P x ≤ η -1 } (17) S 0 = {x ∈ ℜ n ; x ′ P x ≤ 1} (18)
Proof: Consider a quadratic Lyapunov function candidate V defined by V (x) = x ′ P x, P = P ′ > 0, for all x in ℜ n . Regarding the nonlinearity ϕ, the proof relies on the application of the generalized sector condition [START_REF] Tarbouriech | Stability analysis and stabilization of systems presenting nested saturations[END_REF] stating that for any diagonal positive define matrix T , one satisfies:

-2ϕ ′ T (sat(Ky) + Gx) ≥ 0 for any x ∈ S(G, u 0 ) = {x ∈ ℜ n ; -u 0(i) ≤ G (i) x ≤ u 0(i) , i = 1, ..., m}.
Furthermore, the satisfaction of relation [START_REF] Tarbouriech | Stability analysis and stabilization of systems presenting nested saturations[END_REF] ensures that the set S 1 , defined in ( 17) is included in the region S(G, u 0 ) in which the sector condition on ϕ is satisfied. Hence, in S 1 , we can write V (x) ≤ V (x) -2ϕ ′ T (sat(Ky) + Gx). Note that the satisfaction of relation [START_REF] Tarbouriech | Stability analysis for linear systems with input backlash through sufficient LMI conditions[END_REF] guarantees that S 1 contains S 0 .

We want to verify that there exists a class K function α such that V (x) ≤ -α(V (x)), for all x such that x ′ P x ≥ 1 and x ′ P x ≤ η -1 (i.e. for any x ∈ S 1 \S 0 ), and for all nonlinearities Ψ satisfying Lemma 1 (i.e. satisfying relations (9), ( 10), ( 11)).

By using the S-procedure, it is sufficient to check that L < 0, where

L = V (x) -τ 1 (1 -x ′ P x) -Ψ ′ T 3 Ψ + ρ ′ LT 3 Lρ -2( Ψ + L ẇ) ′ N 1 Ψ -2( Ψ + L ẇ) ′ N 2 ( Ψ + (1 -N 3 )L ẇ) -2ϕ ′ T (sat(Ky) + Gx) (19)
with τ 1 a positive scalar and T 3 a positive diagonal matrix.

Choosing N 3 = 1, noting that ẇ = EA 0 x + EBϕ + EBKDΨ, from the definitions of ϕ and of Ψ in ( 8) and in [START_REF] Hu | Control systems with actuator stauration: Analysis and design[END_REF], noting that

V (x) = x ′ (A ′ 0 P + P A 0 )x + 2x ′ P Bϕ + 2x ′ P BKDΨ, it follows that L = L 0 + ρ ′ LT 3 Lρ -τ 1 with L 0 =     x ϕ Ψ Ψ     ′ M 1     x ϕ Ψ Ψ    
and M 1 defined in [START_REF] Park | The asymptotic stability of nonlinear (Lur'e) systems with multiple slope restrictions[END_REF]. The satisfaction of relations ( 12) and ( 14) implies both L 0 < 0 and ρ ′ LT 3 Lρ -τ 1 ≤ 0, and then L < 0, for all (x, ϕ, Ψ, Ψ) = 0.

Therefore, the satisfaction of relations ( 12)-( 16) ensures that there exists ε > 0, such that L ≤ -ε (x ′ ϕ ′ Ψ ′ Ψ′ ) ′ 2 ≤ -εx ′ x. Hence, since by definition one gets V (x) ≤ V (x)τ 1 (1 -x ′ P x) ≤ L, one can also verify

V (x) ≤ -εx ′ x , ∀x such that η -1 ≥ x ′ P x ≥ 1 (20)
Consider now a solution of (1) starting from any admissible initial condition at t 0 such that η -1 ≥ x(t 0 ) ′ P x(t 0 ) ≥ 1. According to [START_REF] Taware | Control of sandwich nonlinear systems[END_REF], there exists a time T ≥ t 0 +(x(t 0 ) ′ P x(t 0 )-1)λ max (P )/ǫ such that x(t) ∈ S 1 , ∀t ≥ T . Furthermore, S 1 is an invariant set for the trajectories of system (1). Hence, in accordance with [START_REF] Khalil | Nonlinear Systems[END_REF], it follows that the trajectories are uniformly ultimately bounded in S 0 . That concludes the proof of Theorem 1.

Using the same tools and assuming that the matrix A is Hurwitz we can propose the following solution to the global case, i.e., to Problem 2.

Theorem 2: Assume that A is Hurwitz. Consider A 0 = A + BK(C + DLE) = A + BK 0 being Hurwitz. A positive scalar τ 1 being given, if there exist a symmetric positive definite matrix P ∈ ℜ n×n , four diagonal positive definite matrices N 1 ∈ ℜ q×q , N 2 ∈ ℜ q×q , T 3 ∈ ℜ q×q and T ∈ ℜ m×m satisfying the following LMI conditions

M 2 < 0 (21) ρ ′ LT 3 Lρ -τ 1 ≤ 0 (23)
with M 2 defined in (22) (see at the top of the next page), then, for any initial admissible conditions (x(0), Ψ(0)), the resulting trajectories of the closed-loop system (1) are uniformly ultimately bounded in the set S 0 defined as follows

S 0 = {x ∈ ℜ n ; x ′ P x ≤ 1} (24)
Proof: The proof follows the same stages as in that one of Theorem 1. In this case one wants to verify that there exists a class K function α such that V (x) ≤ -α(V (x)), for all x such that x ′ P x ≥ 1 (i.e. for any x ∈ ℜ n \S 0 ), for all nonlinearities Ψ satisfying Lemma 1 (i.e. satisfying relations ( 9), ( 10), [START_REF] Nordin | Controlling mechanical systems with backlash: a survey[END_REF]) and any ϕ satisfying -2ϕ ′ T sat(Ky) ≥ 0, ∀x. In the global context, one then chooses G = 0.

Theorem 1 guarantees the finite-time convergence of the solutions to system (1) inside S 0 for any initial condition satisfying (3) and belonging to S 1 . Therefore the property is locally satisfied. Theorem 2 assures the same type of convergence but in a global sense, since it is valid for any admissible initial condition. In both cases, the region S 0 , which captures the closed-loop trajectories, contains the possible equilibrium points or limit cycles induced by the backlash operator. The equilibrium points or limit cycles induced by the saturation operator are outside the set S 1 .

Remark 1: In the case without backlash (Ψ = 0), in the definition of M 1 or M 2 , one removes the columns 3 and 4 and the corresponding lines. One removes relations [START_REF] Tarbouriech | Stability and stabilization of linear systems with saturating actuators[END_REF] and

M 2 =     A ′ 0 P + P A 0 + τ 1 P P B -K ′ 0 T P BKD -A ′ 0 E ′ LN 1 -A ′ 0 E ′ LN 2 ⋆ -2T -B ′ E ′ LN 1 -T KD -B ′ E ′ LN 2 ⋆ ⋆ -T 3 -N 1 LEBKD -D ′ K ′ B ′ E ′ LN 1 -N 1 -D ′ K ′ B ′ E ′ LN 2 ⋆ ⋆ ⋆ -2N 2     (22)
(16) (which, roughly speaking, corresponds to take τ 1 = 0 and η = 1) and therefore one retrieves the results of [START_REF] Tarbouriech | Stability analysis and stabilization of systems presenting nested saturations[END_REF]. Remark 2: In the case without saturation (ϕ = 0), in the definition of M 2 one removes the second column and line, and one sets T = 0. The results of [START_REF] Tarbouriech | Stability analysis and stabilization of systems with input backlash[END_REF] are recovered.

Remark 3: The kind of systems studied in [START_REF] Tarbouriech | Global stability for systems with nested backlash and saturation operators[END_REF] could be addressed by using the framework of the current paper. In [START_REF] Tarbouriech | Global stability for systems with nested backlash and saturation operators[END_REF], the following system is considered:

ẋ = Ax + BΦ[w] w = sat(K 1 y + K 2 Φ[w]) y = Cx + DΦ[w] (25)
which can be rewritten, by noting ∆ = (1 [START_REF] Tarbouriech | Global stability for systems with nested backlash and saturation operators[END_REF], as follows:

-(K 1 D + K 2 )L) -1 from Assumption 2 in
ẋ = (A + BL∆K 1 C)x + BL∆ϕ 2 +B(1 + L∆(K 1 D + K 2 ))Ψ w = ∆ϕ 2 + ∆K 1 Cx + ∆(K 1 D + K 2 )Ψ (26)
with ϕ 2 defined by mimicking ( 5) and Ψ defined in [START_REF] Hu | Control systems with actuator stauration: Analysis and design[END_REF]. According to our framework, in the global case (i.e. matrix A is Hurwitz) we can relax the condition [START_REF] Macki | Mathematical models for hysteresis[END_REF] in [START_REF] Tarbouriech | Global stability for systems with nested backlash and saturation operators[END_REF] linking up K 1 and K 2 .

IV. NUMERICAL PROCEDURES

The objective related to Theorem 1 is to maximize the region S 1 , while minimizing the region S 0 . This double objective needs to consider a trade-off between two conflicting goals: maximize the region of attraction of the closed loop (namely select a small trace for ηP ) and minimize the region of ultimate boundedness stability (namely select a large trace for P ). Thanks to the use of a unique matrix P to shape both sets S 0 and S 1 , the objective may be expressed as maximizing the distance between those two sets, and a solution to Problem 1 is searched with the following (bilinear) optimization problem: min η under conditions ( 12), ( 14), ( 15), [START_REF] Tarbouriech | Stability analysis for linear systems with input backlash through sufficient LMI conditions[END_REF] (

In Theorem 1, relation [START_REF] Paré | A KYP lemma and invariance principle for systems with multiple hysteresis non-linearities[END_REF] admits two kinds of nonlinearities: the products τ 1 P and G ′ T . These nonlinearities may be managed relatively easily. Actually, T ∈ ℜ m×m is a diagonal positive definite matrix. In the particular case of single-input systems (m = 1), T becomes a scalar and the optimal solution of (12) can be obtained from an iterative line search. For m = 2, a search for the optimal solution over a bi-dimensional grid (composed by the 2 elements of T ) can be considered. For systems presenting high number of inputs (m > 2), a relaxation scheme, which translates the problem into a sequence of iterative LMI problems fixing T or G at each step, can be considered. The convergence of the procedure is always ensured, but not necessarily to the global optimal value. Furthermore, the convergence value depends on the initialization of T or G in the iterative procedure. A grid search may also be used to select τ 1 . This latter parameter being issued from the use of the S-procedure, its influence remains rather limited on the optimal cost of the optimization problem (27) as soon as a feasible value has been selected.

In Theorem 2, the condition (21) becomes a LMI when τ 1 is fixed. The same treatment as in Theorem 1 may be applied.

V. ILLUSTRATIVE EXAMPLES

Example 1: Consider first an illustrative example borrowed from [START_REF] Fong | State feedback stabilization of single input systems through actuators with saturation and deadzone characteristics[END_REF] (Example 6) or [START_REF] Tarbouriech | Ultimate bounded stability and stabilization of linear systems interconnected with generalized saturated functions[END_REF] and adapted to our context (the deadzone is replaced by a backlash). This example mimics the case of a nonlinear actuator involving a backlash and a saturation (see Figure 1). A preliminary solution is first obtained on a rough grid then a more precise grid on a thiner interval is applied to get a solution to Problem 1. To illustrate the influence of parameter τ 1 and matrix T (actually in this example of dimension 1), the cost function η solution to the optimization problem ( 27) is plotted in Figure 2.

The optimal ellipsoids S 0 and S 1 obtained for τ 1 = 0.65 and T = 7.5 are plotted in Figure 3 in solid lines. Trajectories starting from several initial states (for any given initial state x(0) one selects Φ(0) = Φ[w(0)] verifying (3)) are also plotted in the figure. They illustrate that the outer ellipsoid S 1 obtained is a reasonable approximation of the domain of attraction. Moreover, to illustrate this aspect, the ellipsoidal approximation of the basin of attraction focusing on the case of saturation only is plotted in Figure 3 (dashed ellipsoid). It may also be compared to the outer ellipsoid corresponding to the maximal volume issued from the grid search (dashdot ellipsoid). A zoom of the inner ellipsoid S 0 and, by the way, of the state-space trajectory of the system issued from x(0) = [2 0] ′ (with Φ(0) = 1.46) is plotted in Figure 4. It illustrates that the system does not converge to the origin but to a limit cycle. This limit cycle remains however confined inside the inner ellipsoid S 0 .

Example 2: Let us consider again the system of Example 1 but with a reorganization of the nonlinear element such as to mimic a system with input saturation and output backlash, such as given in Figure 5.

System (1) is defined by the following data:

A = 0.5 -1 1 0.5 ; B = 0.5 1 ; C = 1 0 0 0 D = 0 1 ; E = 0 1 ; K = 0.93 -3.84
with the same u 0 = 4, L = 1 and ρ = 0.4. The linear behavior of both systems of Examples 1 and 2 is exactly the same. However, the presence of the backlash in the output makes the problem much more difficult to solve and feasibility may not be obtained depending on the organization of matrices C, D and E. In the current case, the problem is made feasible by considering that only the second state is affected by a backlash element in its sensor. A solution, similar to the one obtained in Example 1, is obtained for τ 1 = 0.65 and T = 0.55. The feasibility set with respect to τ 1 and T is displayed in Figure 6, so as the associated cost level η (smallest costs in black, largest in white). The time evolution of w = x 2 and backlash Φ[w] is plotted in Figure 7, for the system initialized in x(0) = 1.5 -0.9 ′ (with Φ(0) = -1.3). It illustrates the local behavior and the limit cycle induced by the backlash. The input signal u and its value Ky before the saturation are plotted in Figure 8. It may verified that u is initially saturated.

Example 3: Let us now consider the following system, with backlash in the output and saturation in the input: The open-loop matrix is Hurwitz and Theorem 2 may be applied to solve the global case described in Problem 2. The ellipsoid S 0 solution to the optimization problem (27) with τ 1 = 0.3 selected with a grid search is plotted in Figure 9. Several trajectories are also plotted which enter the ellipsoid and, then, remain inside. Systems presenting a backlash in the loop and a saturated input have been considered in this paper. The stability has been studied for this class of nonlinear systems by exploiting generalized sector conditions of both nonlinearities. By means of a Lyapunov function, some numerically tractable sufficient conditions have been derived for the regional (local) or global stability of a neighborhood of the origin (where the solutions are captured). The main results have been illustrated with different examples, in particular for systems with saturation in the input and backlash in the output. Several generalizations may be considered. Following the research line suggested in this paper, the design of antiwindup compensator using the knowledge on the backlash operator could be pursued in order to alleviate the bad effects due to both saturation and backlash operators.

A =   -2 -1 -0.

For

  conciseness, throughout the paper, we denote Φ instead of ˙ Φ[w], and Φ instead of Φ[w]. Let us define the nonlinearities ϕ and Ψ ϕ = sat(Ky) -Ky (5)

Fig. 1 .D = 1 ;

 11 Fig. 1. Closed-loop system with backlash and saturation in the actuator

Fig. 2 .Fig. 3 .

 23 Fig. 2. Example 1 -influence of parameters τ 1 and T on the cost η

Fig. 4 .Fig. 5 .

 45 Fig. 4. Example 1 -Zoom on the inner ellipsoid S 0 and state space trajectory issued from x(0) = [2 0] ′ (inside S 1 )

Fig. 6 .

 6 Fig. 6. Example 2 -Feasibility set of the optimization problem (27) with respect to τ 1 and T and associated cost level η. The star at τ 1 = 0.65, T = 0.55 corresponds to the optimal solution

2 Fig. 7 .Fig. 8 .

 278 Fig. 7. Example 2 -Time-evolution of the state of w = x 2 and backlash Φ[w] for the system initialized in x(0) = [1.5 -0.9] ′

Fig. 9 .

 9 Fig. 9. Example 3 -Ellipsoidal set S 0 and several trajectories of the closed-loop system