Mohamed Hedi Riahi
email: riahi28@yahoo.fr

Hend Ben Ameur
email: hbenameur@yahoo.ca

Jérôme Jaffré
email: jerome.jaffre@inria.fr

Rachida Bouhlila
email: rjbouhlila@yahoo.fr

Refinement indicators for estimating hydrogeologic parameters

Keywords: Inverse problem, parameter estimation, parameterization, refinement indicators, storage coefficient, hydraulic transmissivity, zonation

We identify simultaneously the hydraulic transmissivity and the storage coefficient in a ground water flow governed by a linear parabolic equation. Both coefficients are assumed to be functions which are piecewise constant in space and constant in time. Therefore the unknowns are the coefficient values as well as the geometry of the zones where these parameters are constant. The identification problem is formulated as the minimization of a misfit least-square function. Using refinement indicators, we refine the parameterization locally and iteratively. We distinguish the cases where the two parameters have the same parameterization and different parameterizations.

Indicateurs de raffinement pour une paramétrisation optimale dans l'estimation de coefficients hydrogéologiques

Résumé :

On identifie simultanément les coefficients d'emmagasinement et de transmissivité hydraulique dans un écoulement souterrain gouverné par une équation parabolique linéaire. Ces deux paramètres sont supposés être des fonctions constantes par morceaux en espace. Les inconnues du problème sont non seulement les valeurs de ces coefficients mais aussi la geométrie des zones dans lesquelles ces coefficients sont constants. Le problème est formulé comme la minimisation d'une fonction de moindres carrés calculant la différence entre les mesures et les quantités correspondantes évaluées avec la valeur courante des paramètres. L'objectif principal de cet article est la construction d'une technique de paramétrisation adaptative guidée par des indicateurs de raffinement.

Mots-clés : Problème inverse, estimation de paramètres, coefficient d'emmagasinement, transmissivité hydraulique, paramétrisation, indicateurs de raffinement, zonation

Introduction

For many hydrogeological problems, such as management and protection of water resources, we need to study the processes of subsurface flow and of transfer of solutes in the subsurface. The corresponding mathematical models involve the hydrological coefficients which characterize different geological zones in the aquifer. Identifying these coefficients by experimental measurements is difficult and have a high cost. In this work, we are interested on solving the inverse problem of identifying these coefficients. There may be one or several unknown coefficients.

The resolution of such inverse problems has been the subject of many studies [START_REF] William | Review of parameter identification procedures in groundwater hydrology: The inverse problem[END_REF][START_REF] Sun | Inverse problems in groundwater modeling[END_REF][START_REF] De Marsily | 40 years of inverse problems in hydrogeology[END_REF].

The assumption that parameters are constant in each zone is already mentioned in [START_REF] Cooley | Incorporation[END_REF] but the set of zones (zonation) was determined by a priori information. The method proposed in this paper consists to build the parameterization in a progressive and iterative way. The first iterative parameterization method is the multiscale method [START_REF] Chavent | Multiscale parameterization for the estimation of a diffusion coefficient in elliptic and parabolic problems[END_REF][START_REF] Liu | A multiresolution method for distributed parameter estimation[END_REF]. The drawback of this method, where the parameterization is uniformly refined during the iterations, is that it may lead to an overparameterization. The adaptive parameterization method proposed in [START_REF] Ben Ameur | Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities[END_REF] avoids this drawback since the refinement is local and depends on some "refinement indicator". This method was later used for a scalar parameter identification in [START_REF] Ben Ameur | Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities[END_REF][START_REF] Grimstad | Adaptive multiscale permeability estimation[END_REF][START_REF] Hayek | An adaptive subdivision algorithm for the identification of the diffusion coefficient in two-dimensional elliptic problems[END_REF][START_REF] Hayek | A new refinement indicator for adaptive parameterization: Application to the estimation of the diffusion coefficient in an elliptic problem[END_REF]. In [START_REF] Ben Ameur | The multidimensional refinement indicators algorithm for optimal parameterization[END_REF], the parameterization method was extended to the identification of a vector parameter, the colour (Red Green Blue), in a problem of image segmentation .

In this work, we apply the adaptive parameterization method guided by refinement indicators to identify simultaneously the storage coefficient and the hydraulic transmissivity in a confined aquifer which are both unknown. The two coefficients are assumed to be functions piecewise constant in space and constant in time. Therefore they present discontinuities at the interfaces between different geological zones where they are supposed to be constant. The unknowns of our inverse problem are the zonation and the values of the parameters in the zones. We develop different strategies to deal with the cases where the two parameters have the same zonation or different zonations. The developed strategies can be used to identify more then two parameters.

We formulate our inverse problem as a minimization problem of a least-squares function defined as a misfit between measurements and the corresponding quantities computed with "current" parameters. We like to keep as small as possible the number of unknowns in the minimization problem in order to ease the minimization procedure. The data are measurements of the piezometric head in a confined aquifer which is the solution of a linear parabolic equation governing the flow.

The idea of the adaptive parameterization method is to refine a current parameterization according to refinement indicators. A refinement indicator indicates the effect on the optimal data misfit of adding some degree of freedom. The unknown parameterization is estimated, through an iterative process, by refining the scale at which the parameter distribution is described and the process is stopped when the refinement of the scale does not induce a significant decrease of the misfit function. An advantage of this method is to avoid overparameterization. At a given refinement level the parameter is estimated by minimizing the least-squares misfit function. Then we compute refinement indicators and we use these indicators to choose the next refinement level.

This paper is organized as follows. After the introduction in section 1, in section 2 we introduce the direct problem and we set the inverse problems as a minimization problem of a misfit function. Then we present the adjoint method for computing the gradient of that function. In section 3 we present the algorithm of adaptive parameterization and we develop different strategies to identify parameterizations. Finally, in the fourth section we present several numerical experiments corresponding to different versions of the parameterization algorithm including cases where the two parameters have the same zonation or where they have two different zonations. The paper ends with concluding remarks.

Problem setting

We consider the following equations modeling groundwater flow in an isotropic and confined aquifer:

S ∂Φ ∂t -div(T ∇Φ) = Q in Ω × (0, t f ) Φ = 0 on Γ D × (0, t f ) (-T ∇Φ) • n = 0 on Γ N × (0, t f ) Φ(0) = Φ 0 in Ω (1)
where Ω is a bounded connected domain of R 2 , the time variable t belongs to the interval (0, t f ), S is the storage coefficient and T is the hydraulic transmissivity, Φ is the piezometric head and Q is a distributed source terms. n is the outer normal to Ω, Γ D and Γ N are a partition of the boundary of Ω denoting respectively Dirichlet and Neumann conditions. Equations (1) correspond to the direct problem where the unknown is Φ, all other variables, in particular, S and T , are known. We are interested on the inverse problem, where S and T , are unknown and we have measurements of Φ. This inverse problem is formulated as a minimization problem of a least-square function J(S, T ), defined as a misfit between measurements of piezometric head and the piezometric head computed with current parameters S and T :

J(S, T ) = 1 2 Nt i=1 m j=1 (Φ(S, T )(t i , M j ) -d obs ij ) 2 (2)
where d obs ij is the piezometric head measured at time t i and at point M j (x j , y j ) and Φ(S, T )(t i , M j ) is the model output for the current coefficients S and T at the same time and the same point.

Estimating the hydraulic transmissivity T and the storage coefficient S amounts to solving Find(S * , T * ) = arg (S,T )∈U ad min J(S, T )

where U ad = {(S, T ) : S min ≤ S ≤ S max , T min ≤ T ≤ T max } is the set of admissible parameters.

Gradient methods are efficient methods for minimizing J. A crucial step in the minimization process is the computation of the gradient of the function J with respect to S and T . We use the adjoint method [START_REF] Chavent | On the theory and practice of non-linear least-squares[END_REF][START_REF] Chavent | Nonlinear Least Squres for Inverse Problems: Theoretical Foundations and Step-by=Step Guide for Applications[END_REF][START_REF] Sun | Inverse problems in groundwater modeling[END_REF] to compute this gradient and we present it here briefly. The weak formulation of the direct problem 1 is

Find Φ ∈ L 2 (0, t f ; H 1 0 (Ω)) ∩ C 0 (0, t f ; L 2 (Ω))
as

Ω S∂ t Φv + Ω T ∇Φ∇v = Ω Qv ∀v ∈ H 1 0 (Ω) Φ |t=0 = Φ 0 (4)
where L 2 (0, t f ; H 1 0 (Ω)) is the space of square-integrable functions on (0, t f ) with values in H 1 0 (Ω) and C 0 (0, t f ; L 2 (Ω)) is the space of continuous functions on (0, t f ) with values in L 2 (Ω).

Inria

To minimize J under the constraint (1) we introduce the Lagrangien L: L(S, T ; Φ, Φ * ) = J (S, T )+

t f 0 Ω S∂ t ΦΦ * dt + t f 0 Ω T ∇Φ∇Φ * dt - t f 0 Ω QΦ * dt. (5)
Note that if Φ(S, T ) satisfies (4) then Equation (5) becomes

L(S, T ; Φ, Φ * ) = J(S, T ) ∀Φ * ∈ H 1 0 (Ω). (6)
From Equation ( 6) it follows

∂J ∂T δT + ∂J ∂S δS = ∂L ∂Φ δΦ + ∂L ∂Φ * δΦ * + ∂L ∂T δT + ∂L ∂S δS. (7)
We choose

Φ * ∈ H 1 0 (Ω) such that ∂L ∂Φ (Φ(S, T ), Φ * ; S, T )δΦ = 0 ∀δΦ ∈ H 1 0 (Ω). (8)
This choice leads to a simplification of Equation [START_REF] Cooley | Incorporation[END_REF]. Φ * is the adjoint state and Equation ( 8) is the adjoint equation which we can rewrite as

-S ∂Φ * ∂t -div(T ∇Φ * ) = Nt i=1 m j=1 (Φ(t i , x j ) -d obs ij ) dans Ω × (0, t f ) Φ * = 0 in Γ D × (0, t f ) (-T ∇Φ * ) • n = 0 on Γ N × (0, t f ) Φ * (t f ) = 0 in Ω (9)
Taking into account Equations ( 7), [START_REF] De Marsily | 40 years of inverse problems in hydrogeology[END_REF] the gradient of J can be written as

∂J ∂T δT + ∂J ∂S δS = ∂L(S, T ; Φ, Φ * ) ∂T δT + ∂L(S, T ; Φ, Φ * ) ∂S (H, p; S, T )δS = t f 0 Ω δT ∇Φ∇Φ * dΩdt + t f 0 Ω δS ∂Φ ∂t Φ * dΩdt. (10)
It follows that

∂J ∂T = t f 0 Ω ∇Φ∇Φ * dΩdt ∂J ∂S = t f 0 Ω ∂Φ ∂t Φ * dΩdt (11)
Equations ((1)) and (( 9)) are discretized using the finite element method in space. For the numerical solution of the direct problem 1 and the adjoint problem 9, we use the software SUTRA [START_REF] Clifford | A finite-element simulation model for saturated-unsaturated, fluid-densitydependent ground-water flow with energy transport or chemically-reactive single-species solute transport[END_REF] for simulations of groundwater models in two or three dimensions. The minimization of the misfit function J is obtained by the code N2QN1 [START_REF] Lemaréchal | Modulopt, N2QN1[END_REF]. This optimizer implements a quasi-Newton (BFGS) method with line search. sectionParameterization Choosing the parameterization is a crucial point in the solution of inverse problems. In hydrogeology many techniques werere developped. As examples close to our problem we cite [START_REF] Cooley | Incorporation[END_REF] which uses piecewise constant coefficients with a zonation obtained from a priori information, [START_REF] William | Aquifer parameter identification with optimum dimension in parameterization[END_REF] presenting an interpolation method and [START_REF] Chavent | Multiscale parametrization for the estimation of a diffusion coefficient in elliptic and parabolic problems[END_REF][START_REF] Liu | A multiresolution method for distributed parameter estimation[END_REF] proposing a multiscale parameterization in which the refinement process is iterative and uniform.

In this work for each coefficient and at each iteration we build a zonation which is a partition of the domain into zones where the parameter is constant and we compute the value of the parameter in each zone. The proposed method will produce an optimal parametrization, in the sense that it minimizes the number of zones, and so the number of unknowns to explain the available data. The zonations can be viewed as discretization grids for the coefficients. If we had used the computational mesh for the coefficient discretization, the inverse problem would have a very large number of unknowns resulting in an optimization problem of very large dimension. There would not be any balance between the amount of data and the number of unknowns and we would not be able to obtain meaningful results.

The idea of our adaptive parameterization method is to refine the zonation iteratively by adding at each iteration only one degree of freedom obtained by dividing one zone of the current zonation into two zones The choice of the zone to be refined and the choice of the refining cut are given by the refinement indicators. The new zone is added where the refinement indicators indicate that it should induce a significant decrease of the misfit function. So the refinement of the parametrization is not arbitrary and not uniform. In this work we suppose that interfaces between zones are carried by the edges of the computational mesh. The refining cuts are thus carried by the edges of mesh cells and belong to a predefined set of curves.

Refinement indicators

In this section, following [START_REF] Ben Ameur | Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities[END_REF], we recall what is a refinement indicator and how it allows to select the refinement which is likely to induce the largest decrease of the misfit function. Let us consider an example where the initial zonation includes only one rectangular zone Z 0 (see figure 2) and we will proceed to its refinement by cutting it into two rectangular zones Z 1,1 and Z 1,2 .

Let us suppose that the vectorial parameter is constant in the hole domain Ω = Z 0 , we denote by m 0 this parameter. In our case m 0 = S 0 T 0 , where T 0 is the hydraulic transmissivity and S 0 is the the storage coefficient. The optimization problem (3) corresponding to this one-zone parameterization cZ 0 is:

m opt 0 = arg min m0∈U ad J Z0 (m 0 ). ( 12
)
Inria Z 0
Z S1 , Z T 1 . We devide Z 0 into two rectangular zones Z S1,1 et Z S1,2 for the component S and two others Z T 1,1 et Z T 1,2 for the component T , Z 0 = Z S1,1 ∪ Z S1,2 = Z T 1,1 ∪ Z T 1,2 .
New zones are separated by the cuts C S1 and C T 1 . The new zonations are We denote by c = c S c T the discontinuities of S and T respectively across C S1 and C T 1 :

Z S1 = {Z S1,i } i=1,2 and Z T 1 = {Z T 1,i } i=1,2 . We denote by Z 1 = {Z S1 , Z T 1 }. The unknown of the optimization problem (3) corresponding to the refined zonation is m 1 = (S 1,1 , S 1,2 , T 1,1 , T 1,2 ) t . C S1 Z S1,1 Z S1,2
c S = S 1,1 -S 1,2 , c T = T 1,1 -T 1,2 .
Introducing the matrix A = 1 -1 0 0 0 0 1 -1 and assuming that c is known the estimation of m1 comes to solve an optimization problem under the constraint Am

1 = c m opt 1 = arg min Am1=c J Z1 (m 1 ) (13)
where J Z1 is the misfit function corresponding to the new zonation Z 1 . We note that the minimisation problem ( 12) is equivalent to (13) when c = 0.

To the minimization under constraint problem [START_REF] Liu | A multiresolution method for distributed parameter estimation[END_REF] we associate the Lagrangian function

L Z1 (m 1 , λ) = J Z1 (m 1 ) + λ Z1 , Am 1 -c (14)
where

λ Z1 = λ S,Z1 λ T,Z1
is the Lagrange mutiplier associeted to the contraint Am 1 = c. The optimality conditions associated to (13) are then :

∂L Z1 (m opt 1 , λ opt Z1 ) ∂m 1 = ∇J Z1 (m opt 1 ) + A t λ opt Z1 = 0, ∂L Z1 (m opt 1 , λ opt Z1 ) ∂λ = Am opt 1 -λ opt Z1 = 0. (15)
Taking into account the optimality conditions [START_REF] Clifford | A finite-element simulation model for saturated-unsaturated, fluid-densitydependent ground-water flow with energy transport or chemically-reactive single-species solute transport[END_REF] and deriving the Lagrangian function L Z1 with respect to c, we obtain

∂L Z1 (m opt 1 , λ opt Z1 ) ∂c c=0 = ∂L Z1 (m opt 1 , λ opt Z1 ) ∂m 1 c=0 ∂m 1 ∂c + ∂L Z1 (m opt 1 , λ opt Z1 ) ∂λ c=0 ∂λ ∂c = 0
Using ( 14) we deduce

∂L Z1 (m opt 1 , λ opt Z1 ) ∂c c=0 = ∇ c J Z1 (m opt 1 ) c=0 -λ opt Z1 = 0. ( 16
)
Then we have

∇ c J Z1 (m opt 1 ) c=0 = λ opt Z1 . The Taylor expansion of J Z1 as a function of the discontinuity c in a neighborhood of c = 0 is J Z1 (m opt 1 ) c = J Z1 (m opt 1 ) c=0 + c, ∇ c J Z1 (m opt 1 ) c=0 + o(||c||). Since J Z1 (m opt 1 ) c=0 = J 0 (m opt 0 ) we conclude that c, ∇ c J Z1 (m opt 1 ) c=0 = c, λ opt Z1 modelize,
at a first order, the variation between J Z1 (m opt 1 ) and J Z0 (m opt 0 ) in a neighborhood of c = 0.

Components of λ opt

Z1 give two refinement indicators I Z1 S , I Z1 T corresponding respectively to the cuts C S1 and C T 1 inducing the zonation Z 1 . Without solving the optimization problem [START_REF] Liu | A multiresolution method for distributed parameter estimation[END_REF] computing λ opt Z1 allows to conclude if the zonation refinement induced by the cuts C S1 et C T 1 would produce a significant decrease of the misfit function or not. The Lagrange multiplier λ opt is computed using Equation [START_REF] Clifford | A finite-element simulation model for saturated-unsaturated, fluid-densitydependent ground-water flow with energy transport or chemically-reactive single-species solute transport[END_REF]. Then λ opt Z1 = (λ opt S,Z1 , λ opt T,Z1 ) where

λ opt S,Z1 = ∂J Z1 ∂S 1,1 (m opt 1 ) = - ∂J Z1 ∂S 1,2 (m opt 1 ), λ opt T,Z1 = ∂J Z1 ∂T 1,1 (m opt 1 ) = - ∂J Z1 ∂T 1,2 (m opt 1 ). (17)
We note that for the computation of the refinement indicator we use the gradient of J Z1 with respect to the values of the coefficients in the zones. This gradient can be easily calculated from "the fine gradient" computed on the computing mesh:

∂J Z1 ∂S 1,i (m opt 1 ) = K⊂Z S1,i ∂J Z1 ∂S K (m opt 1 ), ∂J Z1 ∂T 1,i (m opt 1 ) = K⊂Z T 1,i ∂J Z1 ∂T K (m opt 1 ), i = 1, 2 (18)
where K denotes a cell of the computing mesh. Thus using the adjoint method (see Section 2) is crucial for an efficient computation of the gradient and of the refinement indicators. In this work, to refine a zone by dividing it into two zones we use only vertical and horizontal cuts carried by the computational mesh edges (see Figure 4). Other cuts are possible (see [START_REF] Ben Ameur | Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities[END_REF]).

Inria

Algorithm

We are adapting an algorithm initially proposed in [START_REF] Ben Ameur | Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities[END_REF] for the estimation of a scalar parameter in porous media and extended in [START_REF] Ben Ameur | The multidimensional refinement indicators algorithm for optimal parameterization[END_REF] to image segmentation with a vector parameter whose components are the three colors Red, Green, Blue. In image segmentation the direct problem is just the identity model, while for the hydrogeological problem, the model is given by ( 1) which require the solution of a PDE. Therefore the application of the adaptive parameterization algorithm is much more complicated.

We choose an initial zonation Z 0 = Ω. A step of the algorithm is as follows:

1. Given a zonation Z k-1 = (Z S,k-1 , Z T,k-1 ) with k zones for each parameter, minimize J Z k-1 to obtain an optimal parameter m opt k . The stopping criteria for the algorithm are:

Build a new zonation

Z k = ( Z S,k , Z T,
• The value of the misfit function is close to zero.

• The misfit function stops decreasing during a number of iterations.

• The values of the refinement indicators are small.

• The norm of the gradient of the misfit function is small.

RR n°8877

Strategies for the computation of the refinement indicators

In this section, inspired by [START_REF] Ben Ameur | The multidimensional refinement indicators algorithm for optimal parameterization[END_REF], we develop different strategies for the computation of refinement indicators I C k from the Lagrange mutiplier λ opt 17), (18).

Z C k computed in (

Strategy 1: euclidean norm for refinement indicators

In this strategy we assume that the two coefficients S and T have the same zonation, we define the refinement indicators by

I C k = λ opt S, Z C k 2 + λ opt T, Z C k 2 .
Since we have a scalar refinement indicator we can apply the adaptive parametrization algorithm defined in the scalar case [START_REF] Ben Ameur | Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities[END_REF]. At each iteration we perform the following steps:

1. Compute I max k = max C∈C ad I C k .
2. Introduce the subset of cuts This strategy is efficient if both parameters S and T have the same zonation.

C ad δ = {C ∈ C ad : I C k ≥ δ * I max k } to 0.5 ≤ δ ≤ 1.

Strategy 2: best cut for all parameters

In this strategy the algorithm is guided by two refinement indicators, one for each coefficient S and T . We select two sets of possible cuts according to the two indicators. During an iteration, a temporary selection is made with respect to only one of the two parameters. At the end of the iteration the same cut is applied to both parameters. It is chosen according to the criterion of "best decrease" of the misfit function. At each iteration we start by the following steps: Then, after obtaining a cut for each parameter, we apply the following steps:

1. Minimize the misfit function with the zonation obtained by applying the cut C * S for both parameters S and T . Thus, going from one iteration to the next, we introduce the same cut for S and T .

Strategy 3: best cut for each parameter

In this strategy, we treat the parameters S and T independently. The choice of this strategy is justified by the case where two geological zones for the porosity and the permeabilities are different. The zonations are constructed in the following way.

Numerical experiments

We consider as an example the aquifer of Rocky Mountain Arsenal, Denver, Colorado in USA [START_REF] Clifford | A finite-element simulation model for saturated-unsaturated, fluid-densitydependent ground-water flow with energy transport or chemically-reactive single-species solute transport[END_REF] represented in Figure 5. This aquifer is horizontal and rectangular. The computational domain is discretized with a mesh of 320 square elements and 357 nodes. The thickness of the aquifer is equal to 40f t. In this aquifer, there are four wells, three pumping wells with a constant rate Q out = -0.2f t -3 /s and an injection well with a constant rate Q IN = 10f t -3 /s (see Figure 5). The lateral boundaries are impermeable. On the top boundary, the piezometric head is set constant and equal to 250.0f t. On the bottom boundary, the piezometric head is varying linearly from left to right from 17.5f t to 57, 5f t.

In the first following numerical experiences, data are numerical measurements of the piezometric head at each node of the computing mesh and at each time step. In later experiments the number of measurements will be reduced. These measurements will be observations of the numerical solution obtained with the exact coefficients.

The exact hydraulic transmissivity and storage coefficient are constant in the whole aquifer with values T = 2.5.10 -4 m/s, S = 6.0010 -7 , except in inclusions where their values are T = 2.5.10 -6 m/s, S = 9.9510 -7 . We consider a case with the same zonation for both coefficients, with two inclusions as shown in Figures 5,6. We will also consider a case where each coefficient has its own inclusion as in Figure 11. We will compare results obtained with the three different strategies for computing the refinement indicators presented in section 2.3.

Experiments using the euclidean norm for refinement indicators (strategy 1)

In a first test we consider exact coefficients S and T with the same zonation shown In Figure 6 which shows also the computational mesh. Figure 7 shows the piezometric head computed with the exact coefficients which is used as data d obs . We observe two depressions corresponding to the two inclusions with very low permeability compared to the surrounding domain and we notice a peak located at the injection well. We identify the zonation of the coefficients S and T and their values in each zone. We use the algorithm described in section 2.2 with strategy 1 where the refinement indicator is defined as the euclidean norm of the vector indicator I = (I S , I T ).

We start with a parameterization where the two coefficients S and T are constant in the whole domain (Figure 8). With strategy 1 the zonation is the same for both coefficients at all iterations.

Going from one iteration to the next, the current zonation for the two coefficients is refined by introducing one cut corresponding to the largest decrease of the misfit function. of color corresponds to a change of zone. We remark that in the third iteration, we already identified edges limiting a large zone sourrounding the two inclusions that we are looking for. At the 9th iteration, we observe that both sought inclusions start to appear and they are identified at iteration 14. The exact parameter values are obtained at the 16th iteration. Note that the exact T is obtained faster than the storage coefficient. This behavior will be observed in all experiments and is due to a smaller sensitivity of the system to S than to T . We initialize the algorithm with a parametrization where the two coefficients are constant in the whole domain (Figure 12). Figure 13 shows the values of the coefficients obtained during the iterations. The algorithm identifies T at iteration 17 but then fails to identify S . Figure 14 shows the final zonation for the two coefficients after 20 iterations and the variation of the normalized misfit function during the iterations.

Experiments using the best cut for the two coefficients (strategy 2)

In this section we take the same data as in tests 1 and 2 but we apply strategy 2: going from one iteration to the next, we add the same cut for both coefficients corresponding to the largest decrease of the misfit function. Again strategy 2 calculates the same zonation for both coefficients. Test 3 corresponds to the case where the two exact coefficients have the same zonation (Figure 6) and test 4 corresponds to the case where the two exact coefficients have different zonations (Figure 11). Initializations are the same as in test 1 and test 2, with constant coefficients over the whole domain (Figures 8 and12).

For test 3 Figure 15 gives the evolution of the computed coefficients during the iterations. It took 24 iterations to obtain the exact coefficients, a few more iterations than in test 1. Again T was identified first. Figure 16 shows the final zonation obtained after 24 iterations and the variation of the normalized misfit function during the iterations. We apply now strategy 2 to the case where the two coefficients have different zonations. As before, the algorithm is initialized with a parameterization in which the two coefficients are constant. The initial values of the two coefficients are shown in 12. 17 shows the variation of the computed coefficients during the iterations. At iteration 10 we identify the zonation and the exact values for the hydraulic transmissivity. For the storage coefficient, it took more iterations to find its zonation and its value. Like for other tests before this is due to the smaller sensitivity of the system to this coefficient. Figure 18

Experiments using the best cut for each parameter (strategy 3)

The main feature of this strategy is that during the algorithm iterations we treat the two coefficients independently. Each coefficient has its own zonation constructed with its own refinement indicator as described in section 2.3.3. We start with the case where the two coefficients have the same zonation 6. Figure 19 shows the variation of the computed coefficients during the iterations.

Figure 20 shows zonations for S and T after 28 iterations and the variation of the normalized misfit function during the iterations. The algorithm identifies T after 24 iterations but fails to identify S. We apply now strategy 3 to the case with exact coefficients having different zonations (see Figure 11). As before, the algorithm is initialized with a single zonation for each coefficient (figure 12). Figure 21 shows the variation of the computed coefficients during the iterations. Figure 19 shows the variation of the computed coefficients during the iterations. Figure 22 shows zonations for S and T after 12 iterations and the variation of the normalized misfit function during the iterations. The algorithm identifies S and T after 12 iterations.

Summary of results

Results for convergence for different strategies and different exact zonations are summarized in Table 1. They show that strategy 1 works well when both exact coefficients have the same zonation while strategy 3 works well when the two coefficients have different zonations. Strategy 2 works well in both cases. In all cases the smaller sensitivity of the system to S than to T is observed by the higher number of iterations for convergence for S than for T , and even in some instances the failure to converge for S while the algorithm was converging for T .

Influence of the number of measurements

Because in more realistic cases, we do not have measurements in each cell of the computing mesh, we study the behaviour of the algorithm with respect to a diminution of the number of measurements. We consider the two cases where the parameters have the same or different zonations. For each case we apply the more efficient strategy that are strategy 1 for the first case and strategy 3 for the second case. Two tests are investigated for each case: in a first test we take one observation point out of 2 cells and in a second test we take one observation point out of four cells as shown in Figure 23. We do not decrease the time sampling of data, because the physical properties of our model correspond to a situation of a saturated porous media where the parameters are not time dependent. In order to evaluate the fit to the data achieved by the algorithm with different number of Inria observation points, we introduce the percentage of explained data

1 - i,j |Φ cal (x j , t i ) -Φ obs ij | 2 i,j |Φ obs ij | 2 .
Case of strategy 1 with exact coefficients having the same zonation (Fig. 6)

In test 7 we use one observation point out of 2 cells as in Figure 23 left while in test 8 we use one observation point out of 4 cells as in Figure 23 right. Figure 24 and Figure 25 show for tests 7 and 8 the values of the obtained coefficients and the decrease of the misfit function when iterating. As before it was easier to recover T than S, but the zonations for both coefficients was well identified though the exact values for exact S were not obtained as expected. The results obtained for these two tests are summarized in Table 2. 3.As before it was easier to recover T than S, but in test 9 the inclusions for each coefficient were well identified with values for close to the exact ones. But for test 10 the situation deteriorates significantly, especially for S.

Conclusion

We presented a method for the simultaneous estimation of the storage coefficient and the hydraulic transmissivity in a confined aquifer. These two coefficients are assumed to be functions piecewise constant in space and constant in time. The unknowns of the inverse problem are the zonations for each coefficient and their values in each zone. The two coefficients may or may not have the same zonation. The method is a generalization of the adaptive parameterization algorithm guided by refinements indicators in [START_REF] Ben Ameur | Refinement and coarsening indicators for adaptive parametrization: application to the estimation of hydraulic transmissivities[END_REF] introduced for scalar parameters to the case of a vector parameter. We consider three variants corresponding to three strategies for choosing the refinement indicators. The first two strategies construct the same zonation for the two coefficients while the last one construct different zonations for the coefficients. Strategy 1 is efficient for the case when the two exact coefficients have the same zonation while strategy 3 is efficient when the two exact coefficients have different zonations. Strategy 3 was working well in both cases while a little slower than the two other strategies when they work well. Obviously the results depend on the sensitivity of the system to each coefficient.

These first results should be completed with more experiments exploring more configurations of the inclusions we are looking for. They could be improved by combining the method that we presented with coarsening indicators in order to decrease the number of zones produced, thus making minimization steps easier. One should also note that our method can be applied for a problem with more than two coefficients to identify.

Figure 1 :

 1 Figure 1: Parametrization and computation mesh (on the left) parametrization and zonation (on the right).

Figure 2 :

 2 Figure 2: Initial zonation (on the left) and computing mesh (on the right)

Figure 3 :

 3 Figure 3: Left: a two-zone zonation Z S1 for the coefficient S. Right: computational mesh for solving the direct and adjoint problems. The cut C S1 is marked in red.

Figure 4 :

 4 Figure 4: horizontal cut (on the left) and vertical cut (on the right).

) 3 .

 3 k ) with k +1 zones for each parameter, by introducing a cut dividing a zone into two zones following four steps: (a) Consider C ad the set of cuts C giving a possible zonation Z C k having k + 1 zones for each parameter obtained by dividing into two zones only one zone in Z S,k-1 and only one zone in Z T,k-1 . (b) For each C ∈ C ad we compute λ opt using (17), (18). (c) For each zone Z C k , C ∈ C ad , compute a refinement indicator I C k from vectors λ opt Z C k according to a choosen strategy proposed in the following Section 2.3. (d) Select a set of zonations Z C k , C ∈ C k ⊂ C ad corresponding to the largest values of I C k and for each of these zonations minimize J Z k . Update the zonation by taking fo each Z k the zonation Z C k which produces the largest decrease of the misfit function.

3 . 4 .

 34 Minimize the misfit function for each partition corresponding to a cut C ∈ C ad δ . Keep only the cut C * which induces the largest decrease of the misfit function.

1 . 2 . 3 .

 123 Compute the refinement indicators: I max S,k = max Introduce the subset of cuts C ad Sδ = {C ∈ C ad S : I C S,k ≥ δ * I max S,k } to 0.5 ≤ δ ≤ 1. Introduce the subset of cuts C ad T δ = {C ∈ C ad T : I C T,k ≥ δ * I max T,k } to 0.5 ≤ δ ≤ 1. 4. Minimize the misfit function freezing the zonation for T at Z T,k-1 and using as zonation for S, the zonations Z C S,k where C ∈ C ad S . Keep only the cut C * S that maximizes the decrease of the misfit function. 5. Minimize the misfit function by freezing the zonation for S at Z S,k-1 and using as zonation for T , the zonations Z C T,k where C ∈ C ad T . Keep only the cut C * T that maximizes the decrease of the misfit function.

Inria

2 .

 2 Minimize the misfit function with the zonation obtained by applying the cut C * T for both parameters S and T . 3. Keep for the two parameters the cut C * = C * S or C * T which induces the largest decrease of the misfit function.

1 . 2 .

 12 Compute the refinement indicators: I max S,k = max Introducing the subset of cuts C ad Sδ = {C ∈ C ad S : I C S,k ≥ δ * I max S,k } to 0.5 ≤ δ ≤ 1.3. Introduce the subset of cutsC ad T δ = {C ∈ C ad T : I C T,k ≥ δ * I max T,k } to 0.5 ≤ δ ≤ 1. 4.Minimize the misfit function by freezing the zonation for T at Z T,k-1 and we use as zonation for S, the zonation Z C S,k for C ∈ C ad S . Keep only the cut C * S that maximizes the decrease of the misfit function. 5. Minimize the misfit function by freezing the zonation for S at Z S,k-1 and we use as zonation for T , the zonation Z C T,k for C ∈ C ad T . Keep only the cut C * T that maximizes the decrease of the misfit function. 6. The new zonation is Z k = ( Z

Figure 5 :

 5 Figure 5: Geometry of the Rocky Mountain aquifer .

Figure 6 :

 6 Figure 6: Test 1: exact coefficients S (left) and T (right).

Figure 7 :

 7 Figure 7: Test 1: piezometric head at 3 times calculated with the coefficients shown in Figure 6.

Figure 8 :

 8 Figure 8: Test 1: optimal values of S (left) and T (right) after the first minimization.

14 Figure 9 :

 149 Figure 9: Test 1: coefficients S (left) and T (right) computed during the iterations with strategy 1, when exact S and T have the same zonation.

Figure 10

 10 Figure 10 shows the final zonation with 16 zones and the variation of the normalized misfit function during the iterations.

Figure 10 :

 10 Figure 10: Test 1: strategy 1 when exact S and T have the same zonation. Final zonation for S and T and variation of the misfit function during the iterations with strategy 1.

Inria 11 :

 11 Test 2: exact coefficients S (left) and T (right) with different zonations.

Figure 12 :Figure 13 :

 1213 Figure 12: Test 2: optimal values of S (left) and T (right) after the first optimization.

Figure 14 :

 14 Figure 14: Test 2: strategy 1 when exact S and T have different zonations. Zonation after 20 iterations (left) and variation of the misfit function during the iterations (right).

Figure 15 :

 15 Figure 15: Test 3: coefficients S (top row) and T (bottom row) computed at iteration 24 using strategy 2, when exact S and T have the same zonation.

Inria

Figure 16 :

 16 Figure 16: Test 3: zonation after 24 iterations for the two coefficients and variation of the misfit function during the iterations.

Figure 17 :

 17 Figure 17: Test 4: coefficients S (on the left) and T (on the right) computed during the iterations using strategy 2, when exact S and T have different zonations.

Figure 18 :Figure 19 :Figure 20 :

 181920 Figure 18: Test 4: strategy 2 when the two exact coefficients have different zonations. Zonation for S and T after 29 iterations (left) and variation of the misfit function during the iterations (right).

12 Figure 21 :

 1221 Figure 21: Test 6: coefficients S (on the left) and T (on the right) computed during the iterations using strategy 3, when exact S and T have different zonations.

Figure 22 :

 22 Figure 22: Test 6: strategy 3 when the two exact coefficients have different zonations. Zonations for S and T after 12 iterations (left and center) and variation of the misfit function during the iterations (right).

Figure 23 :

 23 Figure 23: Points of observations: every other point (left), one point out of four (right).

Figure 24 :

 24 Figure 24: Test 7: strategy 1, when exact S and T have the same zonation. Coefficients S (left) and T (center) computed after 24 iterations and variation of the misfit function (right). Observation points shown in Figure 23 left.

Figure 25 :

 25 Figure 25: Test 8: strategy 1, when exact S and T have the same zonation. Coefficients S (left) and T (center) computed after 26 iterations and variation of the misfit function (right). Observation points shown in Figure 23 right.

Case of strategy 3

 3 with exact coefficients having different zonations (Fig. 11) In strategy 3 we use different zonations for each coefficient. In test 9 we use one observation point out of 2 cells as in Figure 23 left while in test 10 we use one observation point out of 4 cells as in Figure 23 right.Figure 26 and Figure 27 show for tests 9 and 10 the values of the obtained coefficients and the decrease of the misfit function when iterating. Results are summarized in Table

Figure 26 :

 26 Figure 26: Test 9: strategy 3, when exact S and T have different zonations. Coefficients S (left) and T (center) computed after 14 iterations and variation of the misfit function (right). Observation points shown in Figure 23 left.

Figure 27 :

 27 Figure 27: Test 10: strategy 3, when exact S and T have different zonations. Coefficients S (left) and T (center) computed after 15 iterations and variation of the misfit function (right). Observation points shown in Figure 23 right.

Table 1 :

 1 Number of iterations for convergence of S and T for different strategies and different exact zonations.

Table 2 :

 2 Results with a reduced number of measurements, strategy 1 and case of exact coefficients with the same zonation

		1 observation point 1 observation point
		out of 2 cells	out of 4 cells
	Final number of iterations	24	26
	Final number of zones	24	26
	J opt	6.67 10 -5	1.45 10 -2
	Percentage of explained data	0.99	0.97

Table 3 :

 3 Results with a reduced number of measurements, strategy 1 and case of exact coefficients with the same zonation

	Inria

RR n°8877

InriaRESEARCH CENTRE PARIS -ROCQUENCOURTDomaine de Voluceau, -Rocquencourt B.P. 105 -78153 Le Chesnay Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

This work was partially supported by the INRIA project Euromed 3+3 Hydrinv.


Damien Martin-Guillerez

Jérémie Guiochet

David Powell

Experience with a Model-based Safety Analysis Process for Autonomous Service Robot

Safety is a major concern for autonomous systems that physically interact with humans, such as service robots. However, modeling dynamics of such systems is hard so classical safety analysis methods need to be adapted. In this paper, we propose an approach based on a combination of well-known safety analysis techniques. We propose to describe scenarios of use with the common Unified Modeling Language. Risk analysis is then performed using a Preliminary Hazard Analysis, an adaptation of the HAZOP method and the classical Fault Tree Analysis. This paper explains the overal process and illustrates it through the exemple of the MIRAS projects which aims to develop a robotic strolling assistant that will help disabled persons to stand, sit and walk.

Safety assessment process, risk assessment, autonomous systems Hazards List Recommendations Preliminary Hazard list Operational Hazards list 4

I. INTRODUCTION

Safety is now a major concern in many computer-based systems and more particularly for autonomous systems such as service robots in physical contact with human. The traditional approach to analyze safety of such systems is to use methods such as Fault Tree Analysis (FTA) or Failure Mode, Effects, and Criticality Analysis (FMECA). Those methods are usually based on models of the systems such as block diagrams or functional decomposition and automata for dynamics. For autonomous systems it is impossible to represent dynamics with automata as they evolve in an unstructured environment, including humans. A functional decomposition of the decision architecture is also impossible. Moreover, the fact that environment is not structured means that the number of operating conditions is essentially infinite.

We propose an approach to cope with these issues through the combination and adaptation of several well-known techniques. We consider that the earliest models of a system usually describe scenarios of use of the system. We claim that the analysis of deviations of such scenarios allows the identification of major risks. We propose to describe scenarios of use with the common Unified Modeling Language (UML [START_REF]OMG Unified Modeling Language (OMG UML), Superstructure, V2.1.2[END_REF]), and to analyze risks using with a Preliminary Hazard Analysis (PHA), the guideword-based collaborative method HAZOP (HAZard OPerability) [START_REF]Hazard and operability studies (HAZOP studies) -Application guide[END_REF] and Fault Tree Analysis (FTA) [START_REF]Fault tree analysis (fta)[END_REF]. A major advantage of using UML as input model is that it is now a de facto standard for system description, and non-experts can easily understand diagrams such as sequence and use-case diagrams. HAZOP analysis is also well-adapted to the initial steps of the development as it is easily understandable, and through the use of guidewords, it enables a systematic analysis.

This process has been successfully applied to robotic projects and we illustrate each section with the MIRAS project [START_REF] Miras | Multimodal Interactive Robot for Assistance in Strolling[END_REF]. The objective of this project is to develop an assistive robot for standing up, sitting down and walking, and also capable of health state monitoring. It is designed to be used in elderly care centers by people suffering from gait and orientation problems. It is composed of a mobile base and a moving handlebar (Figure 1). This paper is structured as follows. We present our general process in Section II. This approach is detailed step by step with its application to the MIRAS project in Section III. In Section IV, we discuss the validity of our approach. Section V gives an overview of a tool we developed to support this approach. We present related work in Section VI. Section VII concludes this paper.

II. METHOD OVERVIEW

In safety critical systems, safety assessment is usually performed using a safety assessment process [START_REF]Safety aspects -Guidelines for their inclusion in standards[END_REF] where the objective is to reduce the risk to an acceptable level. This process and its terminology is now quite stable in industrial standards [START_REF]Risk management -Vocabulary -Guidelines for use in standards[END_REF]. It is based on a decomposition of activities into risk analysis, risk evaluation and risk reduction [START_REF]Medical devices -Application of risk management to medical devices[END_REF]. Risk analysis aims to identify hazards and estimate the risk. Risk evaluation is a step for comparing the estimated risk against given risk criteria to determine the acceptability of the risk.

Achieving tolerable risk

The following procedure (see Figure 1) should be used to reduce risks to a tolerable level: a) identify the likely user group(s) for the product, process or service (including those with special needs and the elderly), and any known contact group (e.g. use/contact by young children); b) identify the intended use and assess the reasonably foreseeable misuse of the product, process or service; c) identify each hazard (including any hazardous situation and harmful event) arising in all stages and conditions for the use of the product, process or service, including installation, maintenance, repair and destruction/disposal; d) estimate and evaluate the risk (see Figure 1) to each identified user/contact group arising from the hazard(s) identified;

Risk reduction policies application

Sequence diagrams Fig. 2. Adaptation of the standard safety assessment process [START_REF]Medical devices -Application of risk management to medical devices[END_REF] and its generated artifacts Risk reduction is a process in which decisions are made and measures implemented by which risks are reduced to, or maintained within, specified levels. In order to complete risk analysis many techniques have been developed (Fault Tree Analysis, Failure Mode Effects and Criticality Analysis, Event Tree, etc.) and applied in very different safety critical domains from nuclear power-plants to medical robots.

In this context, we based our approach on the classical safety assessment process as described in Figure 2. The same cycle is repeated until the designed system achieves tolerable risk. When applying this process, it is important to:

• describe the target of evaluation at the right level of abstraction; • facilitate communication and interaction between different stakeholders involved in the safety assessment process (e.g., in our case, the stakeholders are the patients, the medical staff, and the robotics experts); • manage the combinatorics of risk analysis, which often results in an excessive number of documents and models; • document safety analysis results and the assumptions on which these results depend to support reuse and maintenance. We chose a subset of UML to describe the system, communicate with the stakeholders and organize safety analysis documents. UML is a standard general-purpose modeling language that includes a graphical notation enabling the representation of an abstract model of a system [START_REF]OMG Unified Modeling Language (OMG UML), Superstructure, V2.1.2[END_REF]. The UML model of a system is composed of different UML diagrams, each of which is a partial graphical representation of the system that concentrates on a particular viewpoint. Two diagrams are commonly used for description of the system usage: use case and sequence diagrams. Use cases represent intended use of the system and are linked with the actors that can trigger scenarios of the use case. Each use case is further documented by fields such as pre and post conditions. Each sequence diagram represents one particular scenario of one use case. Those two diagrams, as presented in the process view of Figure 2, are also the input for hazard identification step.

Hazard identification is then performed using a Preliminary Hazard Analysis (PHA [START_REF] Rausand | System Reliability Theory: Models, Statistical Methods and Applications[END_REF]) and an adaptation of HA-ZOP and UML [START_REF] Martin-Guillerez | A UMLbased method for risk analysis of human-robot interactions[END_REF]. Through the HAZOP method, a system is analyzed by holding a review of the systematic generation of deviations defined by the conjunction of parameters of the system (e.g., pressure, temperature...) and guidewords (e.g., no, more, less...). We apply the PHA in a standard way but adapt the HAZOP method to apply it to UML use cases and sequence diagrams.

Risk estimation is then carried out after the HAZOP analysis. Quantitative risk estimation is however impossible to obtain, mainly because of the impossibility to evaluate probabilities of identified hazards. However, it is important to propose risk reduction means even if probability is not calculated. We cope with this problem by defining different objectives regarding the iterations of the development process. For instance, during the first iteration, only severity (see Table I) is considered, and when the cost (in term of design changes impact, development efforts) is acceptable by robotic experts, risk reduction is carried out. As a result of the first iteration, the document Hazards List is produced (see Figure 2). The next step is based on the PHA and HAZOP-UML hazards, which are analyzed with the Fault Tree Analysis (FTA). This leads to an estimation of the final risks (document Risks List in Figure 2) and also to new recommendations.

III. SAFETY PROCESS STEPS

In this section, we detail the various steps of the safety process. These steps are illustrated by the application in the MIRAS project.

A. Definition of intended use

In the first step, the intended use of the system is modeled using UML use cases and scenarios with UML sequence diagrams. The use cases are completed by conditions of use: preconditions to be fulfilled before any action of the use case can be performed, postconditions to be satisfied at the end of the use case and, invariants that must hold during the use case. Interactions are represented by messages in sequence diagrams. Messages can be annotated (e.g., to specify whether the interaction is physical or cognitive).

These diagrams are often used for high-level representations of a system. They are easy to understand even by non experts so they are suitable as a means for presenting the scenarios to the different stakeholders. Furthermore, they can be included as part of the documentation for the certification process. Sequence diagrams are highly expressive yet can remain quite simple when used to describe use scenarios. This simplicity makes them an attractive support for hazard identification by deviation analysis since it helps to keep the combinatorial aspects of such analysis under control. Throughout the modelling process, preliminary safety remarks and recommendations can be issued.

In the MIRAS project, we identified 11 use cases: Strolling (UC01), Standing up operation (UC02), Sitting down operation (UC03), Balance loss handling (UC04), Summoning and autonomous movement of the robot (UC05), End of use detection and movement to a waiting position (UC06), Positioning the robot by hand (UC07), Alarm handling (UC08), Patient profile programming (UC09), Patient profile learning (UC10), and Robot set-up (UC11). By way of an example, Table II list the conditions put on use case UC02 (Standing up operation) and Figure 3 describes the nominal scenario for this use case. This sequence diagram depicts the interaction between the users and the robot as the potential source of harms. Robot decisions are represented as self-messages.

B. Hazard identification

To identify hazards that can arise from the use of the robot, we used two complementary techniques: Preliminary Hazard Analysis (PHA) and UML-HAZOP (UML -HAZard OPerability).

A standard PHA [START_REF] Rausand | System Reliability Theory: Models, Statistical Methods and Applications[END_REF] is applied at an early stage of the design. The various stakeholders of the project meet together for several workshops (two in the MIRAS project).

TABLE II UC02 "STANDING UP OPERATION"

Use case name UC02. Standing up operation Abstract

The patient stands up with the help of the robot Precondition

The patient is sitting down The robot is waiting for the standing up operation Battery charge is sufficient to do this task and to help the patient to sit down The robot is in front of the patient Postcondition The patient is standing up

The robot is in admittance mode

Invariant

The patient holds both handles of the robot The robot is in standing up mode Physiological parameters are acceptable During the workshops, participants try to consider all the possible causes of hazards in the system (e.g., environmental, electrical, mechanical, hardware/software and human). For each cause, the participants identify the hazards that can arise. In the MIRAS project, the PHA led to the identification of 45 hazards (Tab. III).

We then apply the UML-HAZOP method. We adapt the HAZOP [START_REF]Hazard and operability studies (HAZOP studies) -Application guide[END_REF] method to analyze deviations of the UML use cases and sequence diagrams. According to the Defence Standard 00-58 [START_REF] Defstan | HAZOP studies on systems containing programmable electronics[END_REF], HAZOP analysis is the systematic identification of every deviation of every attribute of every entity. Each deviation is a potential hazard that can lead to a harmful event. We adapted the guideword lists to apply them to attributes of use cases and sequence diagrams. The guideword list we use for the use case entity is given in Table V and an extract of the analysis of UC02 (Standing up operation) is presented in Table VII. All guidewords are applied to generate deviation. The analyst then establishes the effect at the use case level, and the result in the real world. The other columns of the table guide the analyst to establish a severity level, to deduce requirements and otherwise make remarks on that deviation. The complete method is presented in [START_REF] Martin-Guillerez | A UMLbased method for risk analysis of human-robot interactions[END_REF].

In the MIRAS project, the analysis of 297 deviations led to the identification of 13 hazard classes (Table IV). This table presents the main hazardous situations of the system. In the HAZOP analysis, each deviation that potentially leads to a hazard class is labeled with the corresponding number (column "Hazard Number", Table VII). Table VI) gives an extract of the list of recommendations resulting from application of the UML-HAZOP method. This list is derived from the "new safety requirements" column of the UML-HAZOP tables. The conjunction attribute + guideword facilitates the generation of deviations. We adapted the guidewords of Erreur ! Source du renvoi introuvable. according to the two types of entities. They are presented hereafter in tabular form. Each table associates a list of guidewords for an attribute of one entity. An interpretation of the generic deviation is also provided in order to guide the mental process. Table 3 is the adaptation of the HAZOP guidewords for use cases, and Table 4

TABLE V ATTRIBUTES, GUIDEWORDS AND INTERPRETATIONS FOR USE CASE ENTITY IN THE UML-HAZOP METHOD

for sequence diagrams (which is an extension of the results of Lano et al. [LAN02]).

Entity = Use Case

Attribute Guideword Interpretation

No/none The condition is not evaluated and can have any value

Other than

The condition is evaluated true whereas it is false The condition is evaluated false whereas it is true

As well as

The condition is correctly evaluated but other unexpected conditions are true

Part of

The condition is partially evaluated Some conditions are missing

Early

The condition is evaluated earlier than required (other condition(s) should be tested before)

The condition is evaluated earlier than required for correct synchronization with the environment

Preconditions / Postconditions / Invariants

Late

The condition is evaluated later than required (condition(s) depending on this one should have already been tested) The condition is evaluated later than required for correct synchronization with the environment

C. Risk estimation

During the first iteration of the process, a qualitative approach is followed to estimate the risk associated with each hazard class (severity column of Table IV). Safety Integrity Level (SIL [START_REF]Functional safety of electrical/electronic/programmable electronic safety-related systems[END_REF]) requirements are estimated when Revised by: Dev Eval Final

Rec1

The standing-up profile should be validated by a human operator

The profil can be validated only after at least one try X

Rec2

Worst-case elictrical consumption must be evaluated beforehand.

Display of the mean battery time left by the robot.

X

Rec22

Send regularly a network heartbeat from the robot to the medical staff control panel. Launch alarm on time-out.

X

Rec26

Plan a procedure to check that the robot has correct knowledge of the environment.

Display of the map used by the robot on the control panel. performing HAZOP. To do this, an indicative SIL is assigned to each severity level (cf. Table I , SIL column) and thereby to safety-related components associated with each HAZOP deviation (cf. Table IV, integrity level requirements column) according to the most severe hazard induced by this deviation (e.g. the patient position detection system should be SIL1).

During the second iteration of the safety assessment process, we extract FTA top events from the hazard class list of Table IV. The FTA gives top event occurrences and enables the computation of quantitative risk. Fault tree analysis was carried out in the usual way so, given the space limitation, we do not detail it there. The only point we want to stress is the link between FTA and the other artifacts of the safety analysis process.

D. Risk evaluation

As mentioned before, risk evaluation consists in comparing the estimated risk with given risk criteria to determine the acceptability of the risk. Even if some proposals have been made for acceptability criteria [START_REF] Nokata | A safety strategy for rehabilitation robots[END_REF] there is no set of generally accepted risk acceptance principles for service robotics. This is an important issue but it deals with political and ethical concerns. Nevertheless, criteria are needed for engineers to determine which risks have to be reduced. This step is again an iterative process. In the MIRAS study, after the first iteration (coming after a preliminary risk estimation only considering hazard severity), each hazard and possible recommendations were proposed to the robotics experts, and on the basis of our risk analysis results they proposed a classification of risk acceptance according to three versions of the robot : development version for using the robot in the laboratory, evaluation version for the prototype used in hospitals for clinical evaluation, and final version for operational life.

In a second iteration of this activity, the final risks as presented as in Table IV with a severity and an occurrence estimation are classified using the ALARP (As Low As Reasonably Practicable) principle. Three risk level zones are defined : intolerable region, the ALARP region, and the broadly acceptable region. A risk can stay in the ALARP region only if further risk reduction is impracticable or if its cost is disproportionate to the improvement gained. For each risk, the corresponding sources of deviations and hazards, and the associated recommendations are evaluated for judging acceptability. The final result is an argumentation for the final acceptable risk the system.

E. Risk reduction

Recommendations were issued as a result of each step of the safety analysis process. UML modeling, PHA and UML-HAZOP give rises to general recommendations. Integrity level requirements are issued during UML-HAZOP analysis. A specific integrity level requirement leads to specific safety recommendations given by IEC 61508 [START_REF]Functional safety of electrical/electronic/programmable electronic safety-related systems[END_REF]. These recommendations are applied to reduce risks if their current level is not tolerable. Of course, design recommendations should be applied first, then protective devices then, if no other solution exists, information for users [START_REF]Safety aspects -Guidelines for their inclusion in standards[END_REF]. Some risk reduction techniques will reduce the severity of the corresponding hazard (e.g., a bumper added to the robot will reduce the severity of a collision) while others will reduce the probability of occurrence (e.g., infrared sensors to detect the patient's feet will decrease the probability of a collision between the patient and the robot). The next version of the MIRAS robot (Fig. 4) is now under validation. In addition to various ergonomic corrections, the MIRAS robotics experts1 took into account the risk reduction recommendations resulting from the first cycle of our safety analysis.

The new version reduces the severity of hazard classes 8, 9 and 10 and reduces the occurrence rate of hazard classes 1, 6, 7, 8, 9 and 11 (Tab. IV, columns Severity and Occurrence rate). This new version especially applies the recommendations and meets the requirements asked for the evaluation version of the robot (i.e. a robot that will be used in the laboratory and used for clinical evaluation).

IV. EVALUATION OF THE METHOD

We evaluate our approach from four different perspectives: integrability into the development process, and usability, validity and applicability of the method.

Integrability: in our approach, UML design models are shared with the development process. Deviation analysis can be carried out at the same time as design refinement or testing/coding by the development team. The results of risk assessment and of testing can be used to generate another cycle (and modifying either the design or the implementation) or to accept the prototype as a final version. Another noticeable point is that the early integration of the safety assessment process enables the design to be modified as a result of identified risk reduction recommendations.

Usability: the of using this method in the overall process is quite low. For instance, in our first iteration of the cycle, apart from the design that was shared between stakeholders, the only critical path overhead that was induced by the meetings specific to the risk assessment process: two workshops of two hours each for the Preliminary Hazard Analysis and two sessions of two hours to present the outcomes of the UML-HAZOP Analysis and the assessed risks. The time devoted to the safety analysis itself (not to the critical path) consisted of the time to prepare guidelines for PHA workshops (a few hours), the time to carry out the UML-HAZOP analysis (2 weeks) and the time to format the results (a few weeks). UML-HAZOP already proved to be usable even by hand but would be even more so with a tool to assist book-keeping and result formatting [START_REF] Martin-Guillerez | A UMLbased method for risk analysis of human-robot interactions[END_REF].

Validity: Our approach relies on the UML-HAZOP method that identifies a large set of operational hazards. Other hazards, especially environmental hazards should be covered by PHA. Of course, all hazardous situations cannot be forseen however much effort is put into safety analysis. Nevertheless, using these techniques, all classical hazards of robots were identified and several new hazards were discovered (e.g., incorrect position of the patient) together with the corresponding misuse or system failure. This gives use confidence that the coverage of UML-HAZOP and PHA are sufficient to identify most hazards of service robots.

Applicability: the first iteration of our risk assessment process leaded to conclusive evidence that the first prototype of the MIRAS robot was unsafe. The hazard list and the recommendations issued were accepted by the robotic experts in the MIRAS project. Our recommendations were taken into account in the design of the next prototype (which also includes several ergonomic changes).

V. CASE TOOL DESCRIPTION Following two case studies where we use our method to analyse safety of service robots (PHRIENDS [START_REF] Phriends | Physical Human-Robot Interaction: Dependability and Safety[END_REF] and MI-RAS [START_REF] Miras | Multimodal Interactive Robot for Assistance in Strolling[END_REF] projects) we developed a Computer-Aided Software Engineering tool (c.f. Figure 5) to support the method, with the following motivating features:

• Support for UML modelling (Use Case and Sequence Diagrams) • Partial automatic generation of HAZOP tables for systematic analysis • Management of the combinatorial aspects of the HA-ZOP • Guaranteed consistency between UML model and risk analysis tables • Support for building a safety argumentation for risk acceptability • Profiles allowing project-dependent configuration of severity levels, HAZOP table columns and guideword lists The tool is built as an Eclipse plugin (www.eclipse.org) using the Graphical Modelling Framework (GMF). A first prototype has been released implementing most of the steps of our method. Users can draw and document UML use case and sequence diagrams. Based on configurable guideword listings and on UML models, HAZOP tables are partially filled. The user can then define severity levels, and enter all data required for the analysis. Final documents (Hazard List, Recommendations and Required Integrity Levels) are then extracted from HAZOP or PHA Tables, and traceability links are displayed to link hazards to causes. The tool is easy to use because of its simplicity and integration to a common environment (Eclipse). The method and guideword list can be adapted thanks to HAZOP Table templates. Furthermore, the analysis can be exported in CSV to reuse the results outside the tool. However, rich formats like HTML or Excel are not yet available for exportation, which currently limits the integration of our tool with other software. A second prototype is under development integrating visual improvements and a more user-friendly production of documentation.

VI. RELATED WORK

There have been several previous studies aimed at linking model-based development with risk analysis. For instance, in the CORAS project [START_REF] Coras | A platform for risk analysis of security critical systems[END_REF], [START_REF] Bjørn | An approach for modelbased risk assessment[END_REF], developed a framework, exploiting risk analysis and object oriented modelling concepts, for risk assessment of security critical systems. In our case we focus on safety and not on security, but the objectives of our study are quite similar to CORAS. Nevertheless, we do not have the same claims in terms of UML diagrams (we only focus on use case and sequence diagrams) and risk analysis techniques (we only focus on HAZOP and FTA). A major difference is that we strongly interconnect UML models and techniques such as HAZOP whereas that is not the case in CORAS. For instance, they use HAZOP without any real link with UML models (their HAZOP guidewords are not applicable on UML elements). Actually, they identify critical assets and analyse for instance any deliberate/unintentional manipulation of this asset [START_REF] Winther | Security assessments for safety critical systems using hazops[END_REF]. In safety, assets are not elements of the system but they are the system itself or system users. Hence, their approach is hardly applicable for safety analysis [START_REF] Bjørn | An approach for modelbased risk assessment[END_REF].

Our risk analysis approach is based on a re-interpretation of the guidewords for Hazard and Operability Studies (HA-ZOP) in the context of different UML models. The proposal in [START_REF] Lano | Safety and security analysis of object-oriented models[END_REF] followed by a more systematic study in [START_REF] Hansen | Hazop analysis of umlbased software architecture descriptions of safety-critical systems[END_REF], also considers a HAZOP guideword interpretation for the deviations of UML elements such as class, association, classifier role, message, etc. A similar approach was followed in [START_REF] Gorski | Development and validation of a hazop-based inspection of uml models[END_REF] and [START_REF] Jarzebowicz | Empirical evaluation of reading techniques for uml models inspection[END_REF], which also present a statistical analysis of the usability of this method. The guideword interpretation for the static UML diagrams in those studies aims to inspect the model to determine development faults and not to identify operational deviations. Nevertheless, for the UML dynamic diagrams (use case, sequence, activity, and statechart diagrams) many guideword interpretations can be used for exploring deviations during operational life. This is the case in studies presented in [START_REF] Johannessen | Hazard analysis in object oriented design of dependable systems[END_REF] and more formally in [START_REF] Allenby | Deriving safety requirements using scenarios[END_REF], which focus on use cases. The latter study led to a method that has been successfully used in [START_REF] Arlow | Safety specification of the active traffic management control system for english motorways[END_REF] and [START_REF] Frantz | Integrating safety and formal analyses using uml and pfs[END_REF]. This work on use cases also inspired a similar approach for security where new interpretations of guidewords have been proposed [START_REF] Srivatanakul | Security analysis with deviational techniques[END_REF]. Even if this work is more oriented towards malicious behavior of actors, several interpretations can be applied in safety-critical systems with human-machine interactions. We combine and extend the results of those studies, but focus only on use case and sequence diagrams in order to explore deviations during operational life. We also give a particular attention to the integration of HAZOPlike human error analysis techniques as presented in [START_REF] Guiochet | Toward a humancentered uml for risk analysis -application to a medical robot[END_REF]. Indeed, human factors methods [START_REF] Stanton | Human Factors Methods: A Practical Guide for Engineering And Design[END_REF] are a major issue in safety-critical systems but their analysis is often uncorrelated with preliminary system modeling activities. On the contrary, a key point of our approach is to consider human factors from the outset, by including them in model based risk analysis.

VII. CONCLUSION

To tackle safety of autonomous systems, appropriate analysis methods are needed especially when the system physically interacts with humans. Standard risk assessment methods are however limited to simple systems and usual model-based risk assessment methods do not enable to model dynamics. Thus we adapted the classical process with a new model based approach for autonomous systems in physical contact with humans. We model the system using a subset of the well-known standard format UML. We apply PHA and our adaptation of HAZOP to identify hazards. A qualitative method is used to evaluate the risk on the preliminary design allowing the safety process to be integrated early in the development process. FTA is used to evaluate the risk on the other iterations of the safety process. A tool was developed to support the process.

We applied the process to the robot assistant developed in the MIRAS project. The first iteration of the safety process in that project confirmed the needs for high-level design analysis. Furthermore, starting the safety assessment process at the very first step of the design is helpful. In MIRAS, we obtained several results on the first iteration and the recommendations issued in that process enabled integration of safety constraints in the design of the second prototype. Assistant robotic is lacking of standard especially regarding the modeling of humans in the safety assessment process We are currently checking the design for the second prototype with our partners in the MIRAS project and plan to apply a second iteration of the safety assessment process, including a quantitative risk estimation using FTA.

Fig. 1 .

 1 Fig. 1. Robuwalker -First prototype

Figure 1 -

 1 Figure 1 -Iterative process of risk assessment and risk reduction

Fig. 3 .

 3 Fig. 3. Sequence diagram UC02.SD01 giving main scenario of UC02 "Standing up operation"

  of the patient without alarm or with a late alarm. Severe HN5 Physiological problem of the patient without alarm or with a late alarm. Severe HN6 Fall of the patient caused by the robot. Severe HN7 Failure to switch to safe mode when a problem is detected. The robot keeps moving. Severe HN1 Incorrect position of the patient during robot use. Serious HN2 Fall of the patient during robot use. Serious HN3 Robot shutdown during its use. Serious HN8 Robot parts catching patient or clothes Serious HN9 Collision between the robot (or robot part) and the patient. Serious HN10 Collision between the robot and a person other than the patient. Serious HN11 Disturbance of medical staff during an intervention Moderate HN12 Patient loses her balance due to the robot Moderate HN13 Patient fatigue Minor

Figure 4 -

 4 Figure 4 -HAZOP methodology adapted from [MoD Def Stan 00-58:2000]

Fig. 4 .

 4 Fig. 4. Mechanical design of the second prototype of the MIRAS robot

Fig. 5 .

 5 Fig. 5. Main view of the CASE Tool to support the UML-HAZOP method

TABLE I SEVERITY

 I LIST USED IN THE MIRAS PROJECT, DERIVED FROM THE ABBREVIATED INJURY SCALE OF [10])

	Num.	Severity	Type of injury	SIL
	0	None	None	0
	1	Minor	Superficial injury	0
	2	Moderate Recoverable	0
	3	Serious	Possibly recoverable	1
	4	Severe	Not fully recoverable without care 2
	5	Critical	Not fully recoverable with care	3
	6	Fatal	Not survivable	4

TABLE III EXTRACT

 III OF THE HAZARDS IDENTIFIED BY THE PRELIMINARY HAZARD ANALYSIS (45 TOTAL)

				Date:		28/10/09
	Project: MIRAS	PHA Hazards		Prepared by:	Jérémie Guiochet
					Damien Martin-Guillerez
						Version
	Number	Hazard	Categories	Who?	Solution/Answer	Dev	Eval	Final
	H13	Wet ground / Risk of grip loss	Environments					X
	H27	Electrical current drop / Risk of sudden reboot with wrong parameters.	Hardware	Robosoft			X X
					Wheel diameter is reduced to
	H40	Robot runs over patient's foot	Mechanical	Robotsoft	100mm. Foot detection using IR		X X
					sensors.

TABLE VII EXTRACT

 VII OF THE UML-HAZOP ANALYSIS TABLE OF UC02 "STANDING UP OPERATION"

										Date:
										Prepared by:
										Revised by:
										Approved by:
	Line Number	Element	Guideword	Deviation	Use Case Effect	Real World Effect	Severity	Possible Causes	Integrity Level Requirements	New Safety Requirements	Remarks	Hazard Number
											If the robot stops
	15	Battery charge is sufficient to do this task and to help the patient to sit down (precondition)	No/none	Battery charge is too low but the robot starts the standing up operation	The robot interrupts its movement (standing up or walking)	Loss of balance or fall of the patient	Serious	HW/SW Failure Specification error	Battery charge sensors should be SIL2	Worse-case electrical estimation consumption must be of the battery charge evaluated beforehand. Take the lower bound	during standing fall back on the operation, the the patient will most probable scenario is that	2,6
											seat.
	16		Other than	cf L15
	17			Battery charge is high enough but the robot thinks otherwise	Robot refuses to start stand up operation	Patient is confused	None	HW/SW Failure error Specification	None	None

Use case description Use case name: Standing up operation Project

: MIRAS 04/08/09 HAZOP table number: UC02 Damien Martin-Guillerez Entity: UC02 Jérémie Guiochet

TABLE IV HAZARD

 IV CLASSES AND THEIR ASSOCIATED SEVERITIES

Table 3 -

 3 Attributes, guidewords and interpretations for use case entity

TABLE VI EXTRACT

 VI OF RECOMMENDATIONS ISSUED DURING THE UML-HAZOP ANALYSIS(26 TOTAL)

	Date:
	Prepared by:

ISIR -Institut des Systmes Intelligents et de Robotique (http://www.isir.fr) and ROBOSOFT (http://www.robosoft.fr)

ACKNOWLEDGEMENTS

This work was partially supported by the MIRAS Research Project, funded under the TecSan (Technologies for Healthcare) program of the French National Research Agency (French ANR).