A Pouya

M Ghoreychi

G 3s -

DETERMINATION OF ROCK MASS STRENGTH PROPERTIES BY HOMOGENIZATION

Keywords: rock mass, fracture, strength criterion, homogenization, elastoplasticity

A method for determining fractured rock mass properties is presented here on the basis of homogenisation approach. The rock mass is considered to be a heterogeneous medium composed of intact rock and of fractures. Its constitutive model is studied numerically using Finite Element Method and assimilating the fractures to joint elements (Coste 1 ). The method has been applied to a granite formation in France. Geological data on different families of fractures have been used for the statistical representation of the fractures. A mesh-generating tool for the medium with high density of fractures has been developed. The mechanical behaviour of the rock mass (elasticity, ultimate strength and hardening law) has been determined assuming linear elasticity and Mohr-Coulomb strength criterion both for the intact rock and the fractures. Evolution of the mechanical strength in different directions has been determined as a function of the mean stress, thanks to various numerical simulations. The mechanical strength appears to be anisotropic due to the preferential orientation of the fractures. The numerical results allowed us to determine an oriented strength criterion for the homogenized rock mass. A 2D constitutive law for the homogenized medium has been deduced from numerical data. A 3D extension of this model is also presented.

INTRODUCTION

Mechanical behaviour a of fractured rock mass depends on the properties of the intact rock and the fractures. The properties of the intact rock are determined by laboratory experiments. These experiments overestimate the strength of the rock mass at large scale whereas its deformability is underestimated. Rock mass classification methods provide empirical estimations for mechanical properties of fractured rock masses (Barton et al. [START_REF] Barton | Engineering classification of rock masses for the design of tunnel support[END_REF] , Bieniawski [START_REF] Bieniawski | Rock mass classifications in rock engineering[END_REF] ). Semi analytical calculations may provide exact solutions in some special case, for instance when the fractures are parallel to one or two different directions (Bekaert and Maghous [START_REF] De Buhan | Comportement élastique non linéaire macroscopique d'un matériau comportement un réseau de joints[END_REF] , de Buhan and Maghous [START_REF] Bekaert | Three-dimensional yield strength properties of jointed rock mass as a homogenized medium[END_REF] ). For the study of underground facilities in a rock mass with a great number of fractures, numerical models based on Discrete Elements Method (DEM) can be used. However, in presence of a very high density of fractures DEM becomes very heavy to handle and requires very long calculations. An example is given farther ( §3) for a rock mass in a granite located in the West of France. According to geological data, this granite contains a high density of fractures (see Figure 5), mostly filled with recrystalized materials (Calcite), but corresponding to potential discontinuities. One can imagine that due to the very high density of fractures, an attempt to modelling of underground facilities taking account of each fracture may be very heavy, leading to very long calculations. On the other hand, if the dimensions of the underground facilities are great enough compared to the heterogeneity size, a homogenized macroscopic model of the rock becomes easier to handle and more interesting. Many attempts and investigations have been made in this area for materials and rocks (Kachanov 6 , Cai and Horii [START_REF] Cai | A constitutive model of highly jointed rock masses[END_REF][START_REF] Cai | A constitutive model and FEM analysis of jointed rock masses[END_REF] , Le Ravalec and Guéguen [START_REF] Ravalec | High and low frequency elastic moduli for a saturated porous/cracked rock -Differential self consistent and poro-elastic theory[END_REF] , Oda [START_REF] Oda | A method for evaluating crack geometry on the mechanical behaviour of cracked rock masses[END_REF] , Budianski and O'Connell [START_REF] Budianski | Elastic Moduli of cracked solid[END_REF] , Hashin [START_REF] Hashin | Analysis of composite materials[END_REF] , Huet [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF] , …). Most of these attempts have been made for some idealized cases without considering crack interaction. However thus aspect has been taken under account for the study of flow in fractured rocks (Long et al. [START_REF] Long | Porous media equivalent for networks of discontinuous fractures[END_REF] , Long and Withespoon [START_REF] Long | The relation of the degree of interconnection to permeability in fracture networks[END_REF] , Billeux et al. [START_REF] Billeux | Long 'Three dimensional statistical modelling of a fractured rock mass -An example from the Fanay-Augères mine[END_REF] ). More recently homogenization of fractured rocks has been done using Finite Element numerical models. Coste [START_REF] Coste | Comportement Thermo-Hydro-Mécanique des massifs rocheux fracturés[END_REF] has used such a technique for the study of the homogenized hydromechanical behaviour of rock masses using the non linear elastic model of Goodman [START_REF] Goodman | Methods of geological engineering in discontinuous rocks[END_REF] for the fractures. In the present paper, the numerical Finite Element Method has been improved by using an averaging lemma allowing to deduce the homogenized volume quantities from boundary values of nodal forces and displacements. An automatic procedure for mesh generation is also presented for a fractured medium with high density of fractures. Then, the method is applied to the determination of the mechanical properties of a fractured rock mass. The originality of this work is the extension of the previous approach beyond elastic domain, since not only homogenized elastic properties but also mechanical strength, its evolution under bulk stress and as a function of rock deformation (so-called hardening) are determined numerically.

NUMERICAL METHOD 2.1. General scheme

The rock mass is considered as a heterogeneous medium constituted by intact rock and different families of fractures with a statistical distribution of size and orientation and with a given mechanical behaviour. The homogenized strain ε hom and stress σ hom of this medium are defined as the average values of the local strain and stress distributions ε and σ in the medium (Figure 1). The homogenized behaviour of this medium, i.e., the relationships between ε hom and σ hom , are determined by numerical calculations. For this purpose, displacements or tractions corresponding to homogeneous macroscopic strains or stresses are prescribed on the boundary of a sufficiently large domain (Representative Elementary Volume, REV). Strains and stresses in this domain are calculated using Finite Element Method and representing the fractures by joint elements (Coste 1 ). The average values of strains and stresses in the domain are then calculated. The REV size is determined by a preliminary investigation. For this purpose, values of elastic parameters, especially those of the Young modulus, are calculated on domains of increasing sizes. The REV size corresponds to the minimum size beyond which the values obtained for elastic parameters are stable (not fluctuating). In such a condition, the average values of displacements and strains do not show any significant fluctuation under a constant stress (or force) applied on the boundary. Also, the average stress does not change significantly when the boundaries are submitted to constant displacements.

σ hom ε hom σ ε x y Representative Elementary Volume
Figure 1. The fractured rock mass considered as an heterogeneous medium constituted of the intact rock and fractures

Mesh generation

One major difficulty is the mesh generation for Finite Element computation in presence of a high density of fractures. The elements must not contain any fracture inside them ; the fractures have to be located only on the boundary of the elements in order to be assimilated to joint elements. A numerical tool was developed able to generate such a mesh automatically. The automatic procedure for mesh generation is the following : the domain containing fractures is first divided in a number of closed sub-domains by introducing additional lines in the domain which connect the end of the fractures but do not intersect with the existing fractures or lines (figure 2). These closed sub-domains are then discretized individually using classical methods. The meshes corresponding to the individual domains are finally assembled to form the mesh for the whole domain. The interface between different sub-domains is modelled by joint elements only if the interface corresponds to a fracture. Indeed, the division of the whole domain in a set of sub-domains can be made in different ways, therefore, the mesh generation solution is not unique. But, as it is the case of any other mechanical problem studied by FEM, it can be shown that when the mesh is sufficiently fine, the final result does not depend upon the mesh.

Figure 2 : Method of mesh making for the fractured domain by an automatic procedure : subdivision of the domain in closed sub-domains not containing fractures, discretizing subdomains individually and assembling the individual meshes to form the whole mesh.

Averaging method

The transition between "microscopic" and "macroscopic" variables requires averaging strain and stress fields over the whole domain. For this purpose, the following basic lemma (Hill 18,[START_REF] Hill | The essential structure of constitutive structure of constitutive laws for metal composites and polycrystals[END_REF] , François et al. [START_REF] François | Elasticité et Plasticité[END_REF] ) is used :

Let the straightforward volume averages of Cartesian components of Cauchy stress σ and infinitesimal strain ε be denoted by σ and ε . In relation to a given domain, the surface tractions or displacements that would be generated by uniform internal fields σ or ε will be spoken of as uniform loading or uniform constraint respectively. Then, let σ be any bounded and self equilibrated field of stress and ε any field of infinitesimal strain derivable from a continuous displacement. Then, when the boundary conditions for one or the other are uniform, the volume average of their scalar product can be precisely evaluated as the scalar product of their volume average (Hill 18 ) :

σ ε : = σ : ε (1)
Uniform displacement boundary conditions on the boundary ∂V of a volume V are defined by : ∀x∈∂V,

u i (x) = B ij x j ( 2
)
where B is a constant symmetric tensor. Consider now a displacement field u corresponding to the uniform boundary conditions (2) and the strain field ε deriving from u :

ε ij = (∂ i u j +∂ j u i )/2. Integration of ∂ i u j over V gives : V ∫ ∂ i u j dv = u j n i ds = B ik x k n i ds = B ik x k n i ds = V B ij ∂V ∫ ∂V ∫ ∂V ∫
where the following identity has been used :

∂V ∫ x k n i ds = ∂ i x k dv = V δ kj (3) V ∫
This allows to show that : ε = (B+B T )/2

Uniform traction boundary conditions on ∂V are defined by a surface traction density T on ∂V given by : ∀x∈∂V,

T i (x) = C ij n j (x) (4)
Consider a self-equilibrated stress field σ corresponding to the uniform boundary conditions (4). Writing :

V ∫ σ ij dv = σ ik δ kj dv = σ ik ∂ k x j dv = σ ik x j n k ds -(∂ k σ ik ) x j dv V ∫ V ∫ ∂V ∫ V ∫
and taking account of the equilibrium equation in V , ∂ k σ ik = 0, the boundary condition, σ ik n k = C ik n k on ∂V, and of (3), one finds :

σ = C (5)
Consider now a domain V and a loading on ∂V consisting in a surface traction density T or in displacement increments δU. Furthermore, let σ and δε denote the stress and strain fields resulting from this loading in V. The virtual work principle allows to express :

δW e = δW i ⇒ Τ.δ ∂V ∫ U ds = : δε dv σ V ∫
Assuming that the prescribed boundary conditions T or δU are uniform, the above lemma allows to write :

σ ε = V : d V δ ∫ v σ ε : δ = V σ : δ ε
and so the following equality is derived :

1 V ∂V ∫ Τ .δU ds= σ : δ ε (6)
If uniform displacement boundary conditions are prescribed then δU and δ ε are known, and ( 6) allows to calculate σ as a function of the boundary tractions T. In the same way, if traction conditions are prescribed, (6) allows to calculate δ ε as a function of the boundary displacement δU. It can be showed that the above equality holds for combined force and displacement boundary conditions, i.e., when the boundary conditions correspond to u i (x) = B ij x j on a part of ∂V and to T i (x) = C ij n j (x) on the complementary part.

Numerical experiments

Let us assume that the considered domain is a rectangle with side L x in x direction and L y in y direction, and let us number the sides from 1 to 4 as shown in Figure 3. As well, let the following loading cases be considered : a) Compression in x direction with different values of lateral stress σ yy : A constant displacement δU y is prescribed on the side ( 1) and ( 2) with δU y (1)=0, and a constant surface traction T x is prescribed on the sides 3 and 4, with T x (4)=-T x (3). b) Compression in y direction with different values of lateral stress σ xx : A constant displacement δU x is prescribed on the sides ( 3) and ( 4) with δU x (3)=0, and a constant surface traction T y is prescribed on the sides 1 and 2, with T y (2)=-T y (1). c) Shear in xy direction with different values of the mean stress σ m =σ xx =σ yy : A constant displacement δU x is prescribed on the sides ( 1) and ( 2), and a constant displacement δU y is prescribed on the sides ( 3) and ( 4) with δU x (1)= δU y (3) = 0, and δU x (2)= δU y (4), and a constant normal surface traction T is prescribed on all the sides (T = T y (1)= -T y (2)=-T x (4)=T x (3)). These loadings are special cases of combined uniform loadings, the particularity of which is that the prescribed displacement or surface traction is constant on each side of the domain. In these cases, the calculation of the first member of (6) and so the determination of the homogenized quantities are very easy. In fact let δ u (n) designate the average displacement on the nodes of the side number n (n=1,...,4), and by F (n) the total force on the nodes of the side number n (n=1,...,4). Following different loading cases, these quantities are deduced either from the prescribed conditions or from the calculation results. The equation ( 6) and the fact that boundary conditions are constant along each side, allow to obtain different components of δ ε and of σ by means of the following formulas :

δ ε xx = [δ u x (4)-δ u x (3)]/L x , ε yy = [δ u y (2)-δ u y (1)]/L y , δ ε xy = 1 2 {[δ u x (2)-δ u x (1)] /L y + [δ u y (4)-δ u y (3)]/L x } σ xx = 1 2 [ F x (4)-F x (3)]/L y , σ yy = 1 2 [ F y (2)-F y (1)]/L x σ xy = 1 2 {[ F x (2)-F x (1)]/L x + [ F y (4)-F y (3)]/L y }
One can also notice that equilibrium conditions imply the following equalities :

n=1,4

∑ F x (n) = n=1,4 ∑ F y (n) = 0
The homogenized stress and strain are calculated for each load case using the method presented above. Strain-stress curves can be calculated in this way for different loading paths (see Figure 8). The constitutive model for the material (homogenized medium) can be determined through these curves.

FRACTURED ROCK MASS 3.1. Fracture distribution

A 2D representation of a granitic fractured rock mass situated in the West of France is considered here. The x axis in this representation corresponds to the horizontal North-South direction, and the y axis, to the vertical direction. The fractures size, orientation, density and aperture or thickness (filled fractures) introduced in this model have been deduced from the geological data . The possible propagation of existing fractures that may locally occur is not considered. This question is of second order regarding the intensity of existing fractures and also the insufficient data available specially on the extension of fractures. The extension of the existing fractures varies from 20cm to 30m and their slope (angle with the x axis) varies from 0 to 90°. The statistical distribution of the fractures extension is assumed to obey to an exponential law with an average value of 10m. The orientation obeys to a normal distribution law with an average value of 60° (sub-vertical fractures) and a 10° variance. The fractures density is uniform and equal to 0.4 m -2 . The value of the fractures thickness, required for determination of the joint behaviour, is assumed to be about a few millimetres on the basis of some core samples observation. A statistical representation of these fractures is given in Figure 5 for a square domain with 50m side.

Constitutive models of intact rock and fractures

According to the results of laboratory tests on the rock without any fracture and on single fractures, the intact rock has a linear and isotropic behaviour in the elastic stage, with values of Young modulus E = 72000 MPa and Poisson ratio ν = 0.25. The strength properties of the intact rock are given by a Mohr-Coulomb criterion with cohesion C = 17 MPa and internal friction angle φ = 57°. The fractures behaviour in the elastic domain is assumed to be linear and characterized by normal and tangential stiffness k n and k t respectively. Beyond the elastic domain, Mohr-Coulomb criterion is used with cohesion c = 1.51 MPa and internal friction angle ϕ = 27°. These simple models (linear elasticity and Mohr-Coulomb criterion) allows to compare easily the homogenized rock mass behaviour with the behaviour of its constituents (intact rock and fractures). More representative models of fracture behaviour (non linear elasticity) can easily be taken into account if needed. The values of fractures parameters considered here have been obtained by laboratory measurements on some natural fractures with a few millimetres of thickness and filled with recrystallized calcite. Because of the high value of the intact rock cohesion compared to the fractures cohesion (C=17MPa, c=1.51MPa), the fractures reach their elastic limit before the intact rock. Therefore, the damage behaviour of the rock mass is largely determined by the fractures parameters, and in the loading paths presented in the following, The intact rock essentially remains elastic. In order to simplify the calculation, the elastic-brittle behaviours of the intact rock and the fractures are replaced by an elastic-perfectly plastic behaviour. For monotonous increasing loading paths studied here, when the unloading does not take place, this simplification does not change significantly the macroscopic results.

NUMERICAL STRESS-STRAIN CURVES

A preliminary investigation of the REV size for the considered fractured rock mass, by the method described above ( §2.1), led to a side value of about 8m. A square domain of 8m side of was submitted to compressive or shear loading paths described in §2.4. The REV and the mesh generated for this REV are presented in Figures 6 and7 respectively. Homogenized stresses and strains were calculated using the averaging method presented above. The calculations were made under plane stress assumption. This assumption makes the determination of the strength properties easier since this assumption leads to a stress tensor with only three independent components.

To simplify the notations, the homogenized stresses and strains are denoted by σ and ε in what follows. Stress-strain plots obtained for different loading cases are given in Figure 8.

CONSTITUTIVE MODELLING

The stress-strain plots given in Figure 8 provide somehow the "experimental data" on the basis of which the constitutive model of the rock mass is to be determined. These data are to be completed, as it is the case of real experimental data, by further assumptions concerning the expression of the constitutive model.

The analysis of the results shows up three different stages of deformation. The stress-strain plots begin with a linear elastic stage. During the following stage, the slopes of the curve decrease but the stress is always increasing (positive hardening stage). The curves then tend to a horizontal asymptote corresponding to a maximum value of the stress which cannot be exceeded (asymptotic perfect plastic behaviour). In the following, first the results corresponding to the elastic stage, then those characterizing the ultimate plastic behaviour, and finally the results corresponding to the transient hardening stage will be analysed. 4 for the loading paths).

Elastic stage

During the elastic stage, stress-strain curves have a linear trend. This stage extends up to the strain values of about 10 -5 . The ratios obtained for different variables during this stage are given in Table 1. This Table shows that when a compressive load is applied in x or y direction, very small values are obtained for ε xy /ε xx or ε xy /ε yy compared to ε yy /ε xx or ε xx /ε yy . The shear strain ε xy is therefore negligible in these cases. Moreover, the elastic moduli in x and y directions have very close values. For these reasons, the elastic behaviour can be assumed to be isotropic with a good approximation. The average values of Young modulus and Poisson are (Table 1

Elastic limit

Determination of the initial elastic limit is important because it characterizes the initiation of damage which has a big effect on the permeability of the rock mass (important parameter of the disposal safety). The initial elastic limit in different directions can be deduced from the numerical plots. It can also be determined by theoretical calculation based on some simplifying assumptions.

The overall elastic limit corresponds to that of macroscopic (homogenized) stress for under which one of the constituents of the heterogeneous medium reaches its elastic limit. As mentioned above, since the cohesion value of the intact rock is greater than that of the fractures (C=17MPa, c=1.51MPa), the fractures reach their elastic limit prior to the intact rock. Let us consider a fracture with an angle θ with respect to the x direction. Normal and tangent unit vectors are designated respectively by n =(-sinθ, cosθ) and =(cosθ, sinθ). Assuming that the shear and normal stresses on this fracture are given by the projection of the homogenized stress tensor σ on the fracture's plane is equivalent to neglect the interaction between fractures during the initial elastic stage and to assume that they have infinite size. Normal and shear stresses on the fracture, σ n and τ, are then respectively given by σ n = t n .σ. n and τ = .σ. . The elastic domain for the fracture is given by (compressive stress is supposed as negative) :

t n ⎢τ ⎢< c -σ n tgϕ
This inequality leads to the following condition for the homogenized stress σ : ⎢ t .σ. n ⎢ < c -n .σ. n tgϕ (7) Considering that under plane stress assumption, the components of the stress tensor are σ xx , σ yy and σ xy , one can define the mean stress σ m and the deviatoric stress in the plane by the following expressions :

σ m = 1 2 (σ xx +σ yy ) s = 1 2 (σ xx -σ yy ) (8)
Let us consider the domain defined by the inequality (7) in the plane (s, σ xy ) for a given value of σ m and let us define in this plane S = (s, σ xy ), U =(cos2θ, sin2θ) and = (-sin2θ, cos2θ). The inequality ( 7) can be written as :

V ⎢ . ⎢< c -σ S V m tgϕ + S . U tgϕ
This condition defines a domain limited by two straight lines D 1 and D 2 making an angle ϕ with the direction (Figure 9). The elastic domain for the set of fractures is the intersection of the domains defined in this way for different values of θ. For the rock mass considered here θ varies between 0 et 90°, and so, the intersection of the considered domains is a circle centred at the origin O and with a radius R given by (see Figure 9 The elastic domain of the rock mass is the intersection of this domain with the elastic domain of the matrix (intact rock). The latter, defined by the Mohr-Coulomb criterion with cohesion C and an internal friction angle φ, corresponds, for a given value of σ m , to a circle in the plane (s, σ xy ) centred at O and with a radius R' = C cosφ -σ m sinφ. For the values of the parameters c, ϕ, C and φ given above, the second circle contains always the first one (R'>R when R >0), and so their intersection is equal to the first one.

In conclusion, the elastic domain of the rock mass in the stress space (σ m , σ xy ,s) is the circular cone defined by :

Asymptotic stage -strength criterion

As shown in Figure 8, stress-strain curves tend to horizontal asymptotes. The stress value corresponding to this asymptote represents the ultimate strength of the rock mass in the corresponding direction. This ultimate strength depends on the lateral stresses for compressive loads or on the mean stress for shear loads. The results show up a significant anisotropy related to the plastic behaviour of the rock mass. The strength in x direction is in fact about 2 times greater than that of y direction. The ultimate stress versus lateral or mean stress has been plotted in Figure 10. This figure shows that the ultimate strength depends linearly on the mean stress. The equations corresponding to the straight lines are given in Figure 10. The equation of the line σ xy versus σ m in Figure 10 can be written as : σ xy = -σ m tgφ 1 + C 1 . The two other lines of this figure define linear dependency of s upon σ m which can be written as s =-σ m tgφ i + C i with i=2,3. The parameters φ i and C i (i=1,2,3) define the values of the homogenized internal friction angle φ hom and cohesion C hom in different directions. They are given in Table 2 One can notice that the values obtained for homogenized cohesion and internal friction angle are intermediate between those of the intact rock (C =17 MPa, φ = 57°) and the fractures (c = 1.51 MPa, ϕ=27°).

These results can be used to elaborate a model for the ultimate strength of the rock mass. The strength domain in the stress space (σ m , σ xy , s) can be described by a yield function f ; the strength domain corresponding to domain of stresses verifying f(σ m , σ xy , s) ≤ 0. The function f has to be determined using the numerical results given above in Figure 10. According to these results, the intersection of the boundary of the strength domain, corresponding to the surface f(σ m , σ xy , s) = 0, and the planes σ xy = 0 and s = 0 in the stress space, gives the straight lines represented in Figure 10. For a given value of σ m , these results, represented also by the data given in Table 2, define 3 points in the plane (σ xy , s). For σ m = 0, these points are respectively (0, 4.58), (0, -3.16) and (5.02, 0). They can be completed by the point corresponding to the shear strength in -xy direction (shear loading with σ xy < 0, not represented in Figure 8). A numerical simulation led to the value of 2.17 MPa for this strength. The four points obtained in this way have been presented in Figure 11. They can be supposed to form a circle. The position of the centre and the radius of this circle vary with σ m .

Furthermore, an expression for the yield criterion f has to be found verifying these properties :

(a) : the intersections of the surface f = 0 with the planes σ xy = 0 and s = 0 define straight lines in these planes, (b) : the intersection of the surface f = 0 with the plane σ m = 0 defines a circle. An expression for f (a simple one which is not a unique one) can be proposed as following : f(σ m , σ xy , s) = (σ xy -D) 2 + (s -D) 2 -(-Aσ m +G) 2 -D 2 As yield criteria are usually supposed to have the dimension of stress, this expression for f will be changed into the following :

f(σ m , σ xy , s) = ( D) (s D) (-A G) D xy 2 2 m 2 σ σ - + - - + + 2 (11)
Indeed, these two expressions define the same domain in the stress space. The parameters A, D and G must be determined using the numerical data. The comparison between the expressions of f in the cases σ xy = 0 or s = 0 allows to identify A = tgφ hom . The expression given for f assumes that the internal friction angle have the same value in all the (loading) directions. This provides a satisfactory approximation since φ hom has very similar values in different directions (Table 2). So, φ hom can be supposed to be constant and equal to the average value of 31.3°, and then A = tg 31.3° = 0.608. Now, considering the case σ m = 0, the equation f = 0 where f is given by (11), defines a circle in the plane (σ xy , s) centred at the point (D,D) and having a radius R = (G 2 + D 2 ) 1/2 . The values of the two constants D and G are obtained by fitting the four points known in the plane (σ xy , s) (Figure 11).

Hardening stage

The hardening stage begins when the representative stress point reaches the boundary of the initial elastic domain. As seen above, the initial elastic domain is a cone with axis parallel to the σ m axis. For a given value of σ m , the elastic domain in the plane (σ xy , s) is a circle centred at the origin and with a radius R = c cosϕ -σ m sinϕ, where c and ϕ are respectively the cohesion and the internal friction angle of the fractures. This circle is located within the circle corresponding to the ultimate strength (Figure 12). During the hardening stage, the representative stress point follows a path beginning at the boundary of the initial elastic domain and tending to the ultimate strength boundary. This path corresponds to the transient stage described above in the stress-strain plots. Since the initial elastic boundary given by (10) corresponds to the same expression (11) for adequate values of parameters A, D and G, one can assume that the elastic boundary at any stage of deformation can be given by the same expression in which the parameters A,D and G may change. The initial values of these parameters are : A 0 = sin 27° = 0.454, D 0 = 0 and G 0 = 1.51 cos 27° = 1.35 MPa. The final (or asymptotic limit) values are those obtained above : A ∞ = 0.61, D ∞ = 1.07 MPa and G ∞ = 3.73 MPa. Regarding the shape of the strain-stress plots, the evolution of these parameters with strain can be assumed to be given by an exponential expression :

A = A 0 + (A ∞ -A 0 ) (1-e -βξ ) (12)
In this expression, the hardening parameter ξ is defined by ξ ε = p with an initial value ξ = 0, and β is a constant parameter which has to be determined by fitting the stress-strain plots in the transient stage. The fitting of different stress-strain plots corresponding to different loading directions, leads to different values of β varying between 547 and 2447. This means that β depends on the stress orientation. However, this parameter describes the evolution during the transient hardening stage and has no effect on the ultimate strength properties (see Figure 13). Therefore, in case one is interested in the ultimate strength properties of the rock masse, a constant value can be assumed for this parameter in order to simplify the model. A constant value of β = 1000 for all the directions has been considered in the following, and an example of curve fitting using this value of β is shown in Figure 13. Hardening law :

A = (A A ) - 0 e β βξ ξ ∞ - D = (D D ) - 0 e β βξ ξ ∞ - G = (G G ) - 0 e β βξ ξ ∞ - ξ ε = p Initial state : ξ = 0, A = A 0 , D = D 0 , G = G 0 Numerical Values : β = 1000 A 0 =0.454 , D 0 = 0 , G 0 =1.35 MPa A ∞ =0.608 , D ∞ =1.07 MPa, G ∞ =3.73 MPa
Even for 2D modelling, under plane strain or plane stress assumptions, numerical codes usually use 3D constitutive laws. Elaboration of a real 3D constitutive law requires specific investigations and simulation of 3D loading paths. However, it would be possible, thanks to simplified assumptions, to give a 3D extension of the law elaborated above. This extension would not be unique, but it would be very useful for implementation in numerical codes and for modelling of underground structures.

It is worthy of mention that in the model defined above, y axis corresponds to the vertical direction in the rock mass. This axis will be denoted by z in the 3D model. Then it will be assumed that the rock mass has isotropic properties in horizontal planes : all the horizontal directions are equivalent. Therefore an attempt is made to give a 3D expression in which x and y directions (horizontal directions) have equivalent properties. Moreover, the 3D law must give exactly the same 2D law for the plane stress case : substituting σ y = 0 in this expression, the same expression (11) for the yield function must be found with y replaced by z. To obtain this 3D expression one can proceed as following : (σ xx -σ yy ) in ( 8) is substituted first by (σ xx -σ zz ) (substitution of y axis by z) and then by [(σ xx +σ yy )/2-σ zz ] (symetrisation of x and y directions). In the same way σ xy in (11) is substituted first by σ xz and then by (σ xz +σ yz )/2. The σ m is defined by σ m = (σ xx + σ yy + σ zz )/3 , and in order to obtain the same 2D expression for the case of plan stresses (σ y = 0), Aσ m in (11) The values of the parameters A, D, G and β are exactly defined by the same relationships as in the 2D case.

The 3D constitutive model obtained in this way can be introduced in numerical codes in order to study stress and strain fields around underground structures in a rock mass under study. Indeed, such a modelling can be applied in case the characteristic dimension of the structure is greeter than that of the REV. In this case, stress and strain fields can be calculated using the homogenized behaviour model. Only, such a constitutive model is introduced in the numerical code and not the boundary conditions which are used for the computations which one made for the specific purpose of elaboration of such a model.

CONCLUSION AND PERSPECTIVES

The method presented above allows to elaborate a homogenized constitutive model for a rock mass with a high density of fractures. This method is based on a few simplifying assumptions. Its very interesting feature is its capability to take account of the whole available data and information on the rock mass : distribution of fractures and constitutive behaviours of fractures and of intact rock. The authors believe that the suggested methodology provides a much more accurate model than more or less empirical models usually proposed in the framework of « rock mass classification ». An interesting prospect and application of the present work would be to determine the strength properties of typical configurations of fractured rock masses in order to provide classification tables based on numerical calculations instead of empirical estimations. This may require a significant numerical work and a high number of long calculations. Nevertheless, such an effort may provide a very helpful tool for the estimation of fractured rock mass strength and for the prediction of damage risks and the related hydromechanical disturbances around underground facilities.

Figure 3 .

 3 Figure 3. Sides numbers of the sides of the VER rectangle

Figure 4 .

 4 Numerical experiments : (a) Compression in x direction, (b) in y direction, (c) shear in xy direction

Figure 5 .

 5 Figure 5. Fracture distribution for a fractured rock mass deduced from geological data

Figure 6 :

 6 Figure 6 : REV corresponds to a square domain of 8m side

Figure 7 :

 7 Figure 7 : Discretization of the REV

Figure 8 .

 8 Figure 8. Stress-strain curves for a homogenized fractured rock mass corresponding to different loading paths: compression in x direction under different values of lateral stress σ yy , compression in y direction under different values of lateral stress σ xx , and shear in xy direction under different values of the mean stress σ m = σ xx = σ yy (see Figure 4 for the loading paths).

  ) : R = |OM| sinϕ. For the point M, the equation of D U 1 leads to : ⎢ . ⎢= 0 ⇒ S V S . U = -c cotg ϕ + σ m ⇒ |OM|= c cotg ϕσ m and then : R = c cosϕ -σ m sinϕ (9)

Figure 9 :

 9 Figure 9 : Elastic domain of the homogenized rock mass

σ

  xx versus σ yy : σ xx = 4.70 σ yy + 26.13 σ yy versus σ xx : σ yy = 3.49 σ xx + 14.19

Figure 10 :

 10 Figure 10 : Effect of lateral or of mean stress on the strength in different directions

Figure 11 .

 11 Figure 11. The strength domain in (σ m , s) plane is bounded by two straight lines making an angle ±31.3° with -σ m axis, whereas the strength domain in (σ xy , s) plane is approximately a circle.

Figure 12 .

 12 Figure 12. During the hardening the representative stress point follows a path going from the initial elastic boundary to the ultimate strength boundary.

Figure 13 .

 13 Figure 13. Comparison between the numerical plot (corresponding to compression in y direction), the individual fitting of this plot, leading to β = 2447, and fitting of this plot with a fixed value of β = 1000.

  is replaced by 3Aσ m /2. The final expression obtained for the 3D law is :f(σ, ξ) (σ xx + σ yy )/2 -σ zz ]/2 , σ m = (σ xx + σ yy + σ zz )/3

Table 1 :

 1 ) : G = σ xy /2ε xy = 24.7 GPa). This result confirms the validity of the hypothesis of linear and isotropic elastic behaviour.

	Young modulus:	E = 62.6 GPa
	Poisson ratio :	ν = 0.267
	These values lead to the shear modulus value of G =E/2(1+ ν) = 24.7 GPa. This value
	corresponds very accurately to that obtained numerically through a shear loading path (see the
	last line of

Table 1 :

 1 Ratios obtained during the elastic stage

	Load direction xx	σ xx / ε xx (GPa)	ε yy / ε xx	ε xy / ε xx
		63.2	0.269	0.015
	Load direction yy	σ yy / ε yy (GPa)	ε xx / ε yy	ε xy / ε yy
		62.0	0.264	0.005
	Load direction xy	σ xy / ε xy (GPa)	ε xx / ε xy	ε yy / ε xy
		49.4	-0.014	-0.03

Table 2 :

 2 : Homogenized cohesion and internal friction angle in different directions

	Load direction	C hom	φ hom
		(MPa)
	Xx	4.58	33.0°
	Yy	3.16	29.0°
	Xy	5.02	31.8°


Christophe Dony
email: christophe.dony@lirmm.fr

Sylvain Ductor
email: sylvain.ductor@lip6.fr

Zahia Guessoum
email: zahia.guessoum@lip6.fr

Selma Kchir
email: selma.kchir@lip6.fr

Chouki Tibermacine
email: chouki.tibermacine@lirmm.fr

Christelle Urtado
email: christelle.urtado@mines-ales.fr

Sylvain Vauttier
email: sylvain.vauttier@mines-ales.fr

Combining Exception Handling and Replication for Improving the Reliability of Agent Software

Exception handling and replication are two complementary mechanisms that increase software reliability. Exception handling helps programmers in controlling situations in which the normal execution flow of a program cannot continue. Replication handles system failures through redundancy. Combining both techniques is a first step towards building a trustworthy software engineering framework. This paper presents some of the results from the Facoma project. It proposes the specification of an exception handling system for replicated agents as an adaptation of the Sage proposal. It then describes its implementation in the Dimax replicated agent environment.

Introduction

Exception handling and replication are two well-known mechanisms to enhance software reliability, and more particularly fault-tolerance. Replication handles system fail-stops. Exception handling enables programmers to dynamically handle situations that prevent software from running normally. These two mechanisms handle complementary causes of untrustworthiness. Associating exception handling techniques and replication is a first interesting step towards building a trustworthy software engineering framework. Such association is one of the objectives of the Facoma project 4 where it is applied to making reliable multi-agent software.

The Dimax replicated agent platform [START_REF] Marin | Darx-a framework for the fault-tolerant support of agent software[END_REF][START_REF] Guessoum | Adaptive replication of large-scale multi-agent systems: towards a fault-tolerant multi-agent platform[END_REF] is a multi-agent platform that can recover from fail-stops. To do so, the platform seamlessly manages agent replicas. It is thus able to replace an agent that failed by one of its replicas. This replacement is as seamless as possible to software users and does not require any additional code from programmers. A failure is generally detected when an agent fails to answer messages for a given amount of time, either because network connections are lost or because the machine on which the agent ran is down. Active replication systems such as Dimax include algorithms capable of identifying the most critical agents to automatically replicate them. Messages sent to a replicated agent are transmitted to all its replicas, which process the same message in parallel. Responses are all filtered except for those from a single replica, the leader.

Exceptions are those events signaled when the standard control flow of a running service cannot continue. An exception is a kind of answer from an agent running the service. It indicates that the agent is unable to continue its task in the standard way but that it is still alive. The primitives offered by the exception handling system (Ehs) enable agent programmers to signal the problem, to search for handlers and to write handlers that propose alternative behaviors to put the system back into a consistent running state or smoothly terminate its execution.

Exception handling and replication do not apply to the same situations and are of different nature: replication is preventive and exception handling is curative. Both mechanisms are obviously very complementary. If studies can be found on how to combine n-versionning and exception handling [START_REF] Mancini | Exception handling in replicated systems with voting[END_REF][START_REF] Issarny | An exception handling model for parallel programming and its verification[END_REF][START_REF] Romanovsky | An exception handling framework for n-version programming in object-oriented systems[END_REF], combining exception handling and replication appears to be an original and promising solution. This paper describes X-Sage, a specification and concrete running implementation, achieved in the Dimax [START_REF] Guessoum | Adaptive replication of large-scale multi-agent systems: towards a fault-tolerant multi-agent platform[END_REF][START_REF] Marin | Darx-a framework for the fault-tolerant support of agent software[END_REF] context of a programming language running on top of a replication middleware, of such a combination. It extends, adapts and provides a running implementation to our primary specification described in [START_REF] Dony | Specification of an exception handling system for a replicated agent environment[END_REF]. The main issue raised by this combination are: i) to provide agent programmers with an exception handling system that behaves correctly in presence of seamless active replication, and ii) to improve replication strategies, on the base of information conveyed while propagating exceptions from replicas.

Chouki A relire

The remainder of this paper is structured as follows. Section 2 briefly describes the Dimax replicated agent platform which is the context of the implementation. In Section 3, we introduce the Api of our X-Sage exception handling system and an example of a user-defined program that uses X-Sage within Dimax. This example will support the case study discussion in section 5 and shows that replication is seamlessly managed for end-programmers. Section 4 discusses key issues of the implementation of the combination exception/replication : i) how to view a replication manager as a broadcasting agent, ii) how to adapt handlers search to the presence of replication managers and iii) how to deal with exceptions at the replication manager level. In Section 5 we present a case study that highlights the interest of exception concertation made at the replication manager level. Section 4.2 describes the key-aspects of implementation of the X-Sage Ehs within the DimaX platform. Before concluding with open perspectives to this work in Section 7, we discuss related works in Section 6.

Implementation context : the DimaX Replicated Multi-Agent System

The context of this work is the programming of reactive, collaborating agents that are deployed over a middleware which handles agent replication. The concepts exposed in this section are derived from the Dimax software that combines the Dima multi-agent system [START_REF] Guessoum | Adaptive replication of large-scale multi-agent systems: towards a fault-tolerant multi-agent platform[END_REF] and the Darx fault-tolerant middleware [START_REF] Marin | Darx-a framework for the fault-tolerant support of agent software[END_REF].

A Dimax agent is an active and autonomous computation entity that executes in its own thread. Agents interact by exchanging asynchronous messages. Each agent holds a message box and a communication interface to send and receive messages. A specific semantics is associated to messages in order to set up a request / response interaction protocol between agents. This protocol describes peer-to-peer collaborations in which a client agent asks a server agent for a service via a request message. Conforming to a contract-based approach of software, whenever a server agent accepts a request, it commits to send a result, either standard or exceptional, back to the client agent, via a response message. Response messages are correlated with request messages.

Agents are executed on a middleware which provides a fault-tolerant execution context thanks to a replication mechanism [START_REF] Marin | Darx-a framework for the fault-tolerant support of agent software[END_REF]. The execution context consists of a set of distributed replication servers which manage the execution of tasks. Every task belongs to a replication group that identifies the set of tasks which are replicas of a same logical task. All tasks within a replication group thus have the same behavior. They only differ by the environments in which they are running (machines of different kinds that have their own resources). Logical tasks are identified by logical names that are used to send them messages. The replication middleware is in charge of the location and delivery of messages to the corresponding replicas. More precisely, messages are delivered first to the leader of the corresponding replication group. The leader is a replica which has the specific role to control the replication group. The leader forwards the messages sent to the task to its active replicas so that they do the same computation and reach the same new state.

Conversely, all messages sent by replicas are filtered by the replication middleware. Only messages sent by the leader are actually delivered to other tasks. This makes replication transparent to other tasks. Whatever the number of replicas of a task is, a unique message is sent to invoke a computation and a unique message is received as a response. The number and type of replicas is determined by the replication policy, regarding the criticality of the task and the availability of resources (memory, Cpu). In case of failures, new replicas can be dynamically created in order to maintain the redundancy required to provide an expected level of fault-tolerance.When the leader fails, its responsibility is transferred to another replica.If no replica still exists, the task has finally been destroyed by the failure. Every agent executes inside a task. As such, agents can be replicated by the middleware and benefit from this fault-tolerance mechanism. This section describes an exemple of exception handling programming with X-Sage in the Dimax context that will support the following sections discussions. The rationale for the Sage exception handling system used in this example can be found in [START_REF] Dony | Exception handling and asynchronous active objects: Issues and proposal[END_REF][START_REF] Souchon | A proposition for exception handling in multi-agent systems[END_REF]. Its main requirements are:

to enforce agent encapsulation, -to take into account collaborative concurrent activities [START_REF] Romanovksy | Advances in Exception Handling Techniques[END_REF] and to provide mechanisms for their coordination and control [START_REF] Randell | From recovery blocks to concurrent atomic actions[END_REF], -to look for handlers in the runtime history and to execute handlers in their lexical definition context; -to handle concurrent exceptions with resolution functions [START_REF] Issarny | An exception handling model for parallel programming and its verification[END_REF], -and, to support asynchronous signaling and handler search, in order to preserve agent reactivity.

As an illustration, we use the canonical Travel Agency example (cf. Figures 1 and2) in which a Client can send to a Broker a reservation message in order to request a bid for a travel. The contacted broker sends in turn a bid request to several travel providers and airline companies and collects their responses. Then, the Broker selects the best offer and requests the Client and the selected Provider to contract. The underlying request / response interaction pattern of the agent model highlights the role of three key entities: the request, the service and the active agent. The system makes it possible to define handlers at the request, service and agent level and to define resolution functions at the service and agent groups level. This is illustrated in Figure 3 by the colored squares. The position of these squares in the figure represents the order in which these handlers are searched for. The details of handler search, augmented by the adaptation required to take replication into account, are given in Section 4.

Figure 4 shows some parts of the java code of the three agents that contain various exampls of handler definitions.

-Exception handlers can be attached to requests using the @requestHandler annotation. Such handlers can, for example, specify two distinct reactions to the occurrence of two identical exceptions raised by two invocations of the same service.

A handler is attached to a request message within the service contractBroker of class Client (lines 4-6). It states that, for this spcific request, if there is no place left at the given date, it is possible to change the date. -Exception handlers can be attached to services. Such handlers can exceptions that are raised in the dynamic scope of the service's execution.

As an example, a handler is attached to the pollProvider service within the Broker class. It explicitly propagates to the client for him to adapt his request, the exception NoAirportForDestination raised by providers. -X-Sage enables resolution functions to be defined at places where concurrent activities are launched and have to be coordinated, i.e., at the service level) A programmer can define exception resolution functions (based on a model inspired by [START_REF] Issarny | An exception handling model for parallel programming and its verification[END_REF]) using the @serviceResolutionFunction annotation as shown in the example of Figure 4, lines 20-27. such a function is called, is it exists, before any service level handler, each time an exception is trapped at a service level. Its standard behavior is, once all called agents have replied, to aggregate all the exceptions that occurred into a new concerted excep-Fig. 3. Agents are potentially equiped with three kind of handlers tion that conceptually abstract them. Other possible behavior and uses are depicted in [START_REF] Dony | Exception handling and asynchronous active objects: Issues and proposal[END_REF]. -Finally, exception handlers can be associated to agents (see Figure 4, lines 32-38). Such handlers act as if they were repeatedly attached to all of the agent's services. They can be used, for example, to uniformly maintain the consistency of the agent's private data.

We think that these capabilities are powerful enough to encompass most frequent cases the agent programmer is confronted to, and simple enough to be easy to learn and use. The overall contribution of the present is that this example will continue to compute correctly and transparently for the users in presence of replication. The following sections discuss this.

Combining Exception Handling and Replication

EHSs leverage the elements of execution models to provide adapted control structures for managing exceptions. The previous section has presented the X-Sage exception handling system, initially designed to manage exceptions in multiagent systems. It provides means to associate exception handlers to the main elements of the execution model of agents (the requests, the services and the agents themselves). To support replication, X-Sage has been extended to take into account replication groups as a new kind of elements in the execution model. The exception signaling scheme is modified so that the exceptions signaled by replicas are sent to the managers of the replication groups they belong to. A handler and a resolution function are added to the definition of the replication group managers so that they are able to concert and handle exceptions which reflect the status of whole groups of replicas. Handling exceptions at replication manager level enhances both the dependability and performance of applications.

Chouki Je pense qu'il faut preciser en quoi cela ameliore les performances.

Handler Search in a Replicated Multi-Agent Environment

The handler search process for replicated agents is presented on figure 6. As replication is transparent, a replica is a standard agent. So, when an exception signaled to a replica, the handler search process begins with the standard Sage exception handling schemes proposed for agents (see Section 3). A request then a service and finally agent-level handlers are sought in the replica. When the exception is not caught by the handlers defined in the replica, the exception is propagated by the EHS to the manager of the replication group. The handling Fig. 5. The example with Replicated Agents of the exception carries on by the execution of the resolution function associated with the replication manager. The resolution function determines whether it is worth waiting for the results of other replicas before signaling an exception to the client of the replicated agent. Indeed, the failure of one replica does not imply the failure of the whole replication group. This depends on the contextual nature of the exception, as explained in the following subsection. The resolved exception is caught by a generic handler defined in the replication manager. Its role is to signal the resolved exception first to the client of the replicated agent but also to trigger the replication strategies that determine how the organization of the replication group must evolve to maintain the dependability of the replicated agent, depending on the signaled exceptions. Without such a resolution, the exceptions signaled by the leader of the replication group would be immediately sent to the client, taking no advantage of the other replicas to provide a more dependable agent.

Handling Exceptions at the RM Level

When an agent's replica raises an exception, the Ehs invokes the replication manager. The RM will either, as described in the following section, put the system back into a consistent state, signal a new exception to the request caller or propagate one of those trapped by the resolution function. In this latter case, the handler search will continue as explained previously.

Typology of Exceptions. In order to adapt its behavior, the replication-level resolution function of an agent has to determine whether an exception raised by one of the agent's replicas will also be raised by the others. Two categories of exceptions thus have to be identified: replica-specific exceptions (i.e., exceptions raised when some resources specific to a given replica become unavailable) and replica-independent exceptions (i.e., a bad parameter in the request sent to the We thus have designed our algorithms based on a classification of exceptions. Many of such classifications exist. Goodenough's seminal paper [START_REF] Goodenough | Exception handling: Issues and a proposed notation[END_REF] has proposed a first classification in domain, range and monitoring exceptions. This classification is based on "the reason why an exception is raised". It however appears that we have no way to know whether a range exception (for example) is replica-specific or replica-independent. A classification in terms of Error (serious problem that should not be handled) and Exception (business problem that can be handled) as in Java, inherited from the Flavors system, highlights exception gravity but does not answer our question.

A more appropriate classification relies on exception semantics and distinguishes business (also called domain or applicative) exceptions from system (or resource) exceptions. Business exceptions are direct consequences of a programmer's code. Assuming all replicas of a given agent have the same deterministic behavior, exceptions identified as business exceptions can be considered replicaindependent: they are raised by either all replicas of a given agent or by none of them. System exceptions are raised by the runtime environment and are likely to reflect a specific communication or resource access problem. In our context, system exceptions can be raised by the Java virtual machine, the agent language interpretor or the replication middleware. They can be considered as replicaspecific.

In our system, the question of knowing how to distinguish these two exception categories has been resolved statically, by typing exception classes with two distinct interfaces.

Controlling a Set of Replicas. Once the category of the exception is known, the RM reacts adequately.

Handling a replica-independent exception. When a replica-independent exception is raised by a replica, it is immediately propagated to the client agent since all other replicas are expected to raise the same exception. The RM handler does not stop the execution of the request in the other replicas but filters the exceptions they raise (in order not to send the same exception to the client agent several times). This strategy makes it possible to check afterwards (when all requests are processed) if all replicas have actually raised the same exception and thus are in the same consistent state.

Handling a replica-specific exception. When a replica-specific exception is raised by a replica, it is logged by the resolution function of the RM. All remaining replicas pursue their computation. If a normal response is computed by any replica, it is forwarded to the client agent. The other subsequent normal responses are discarded. When all replicas have processed the request (responded either with a normal value or with a replica-specific exception), the RM destroys all the faulty (logged) replicas. If the leader is destroyed, a new leader is chosen among the remaining replicas. If all replicas are destroyed (they all ended with an exception), a service failure exception is signaled to the client agent.

Overview of the Implementation. To implement our Ehs, we created a new SaGEManagerAgent class. All agents of our example application are instances of this class (see Figure 4). All the annotations introduced in Section 3 (@service, @serviceHandler, @agentHandler, @requestHandler and @serviceResolutionFunction) are implemented as Java annotation interfaces. Figure 7 depicts a class diagram of this implementation.

To distinguish request messages from response messages, we have specialized the sendMessage() method of the SaGEManagerAgent class to directly send request messages through Darx classes to agent addresses and to redirect response messages to the DarxCommunicationComponentX class where they will be handled before passing through Darx classes. In the DarxCommunicationComponentX class, response messages sent by replicas are encapsulated in a kind of response messages that we call shadow messages. SageShadowMessage is a subclass of the SageResponseMessage class. Once replicas' responses have been encapsulated in shadow messages, they are transfered to the leader's mailbox using the receiveAsyncMessage(m) method from the AgentManager class.

To summarize, all responses computed by an agent's replicas, either they be normal or exceptional, are sent back to their leader agent as shadow messages. The leader is then able to record, filter, dismiss or concert responses so as to send a unique response to the client (if needed).

When an exception is signaled by a service of a replica, searching for a relevant handler is done thanks to the SaGEManagerAgent class for which the localHandlerSearch(e) method is defined. If no handler is found, the search process is delegated to the Replication Manager of the leader (represented by an instance of the SageReplicationManager class). Depending on the excep-tion category (it being replica-specific or replica-independant), the exception is logged or signaled to the client agent thanks to the callerHandlerSearch(e) method of the SaGEManagerAgent class5 . We used our exception handling system in the implementation of the example introduced in the previous section (the "search for flights" application). The execution of the application passes through the Darx layer to create a name server and to launch the server on local and remote machines. Darx passes then the control to Dimax for agent creation on the local machines. Each agent is encapsulated in a TaskShell which contains specific informations to each agent (identifier, task name, ...) and can manage replication groups. When an agent of a specific task is created for the first time, it creates a replication group and it registers itself as a leader of this task. As mentioned previously in Section ??, the Travel agency and Airline company agents are duplicated.

Case Study: Search for Flights

Once all agents (leaders and replicas) have been created, the service lookForService() of the Client is executed and sends a request message to the service getService() of the Broker. The latter transmits the client's request to the travel agency and airline company agents and invoke their getPrice() service. As those agents are replicated, they forward the message they received from the broker to their replicas (method deliverAsyncMessage() of ActiveReplicationStrategy). Each replica executes its service and returns its response to its leader encapsulated in a shadow message. While executing its service, one of the replicas of the agent Travel Agency signaled the exception NoFlightAvailableException which implements the interface ReplicaIndependant. The execution of the service getPrice() is so suspended and a handler is searched in the agent class. As no handler has been found locally, the search process continues as described in Section 3. The resolution function detects that the type of the exception is ReplicaIndependant and the handler of the replication manager sends the message callerHandlerSearch(e) to the broker which traps the exception in its request handler.

Regarding the agent Airline Company, one of its replicas signals a ReplicaSpecific exception: TemporaryTechnicalProblem. The search for a handler passes through the steps described previously. The resolution function detects that the exception is ReplicaSpecific, log it and wait for other answers of the replicas. All replicas have signaled a Replica-Specific exception then the handler of the replication manager is called. The handler of this RM sends the message callerHan-dlerSearch(e) to the broker like in the previous case (e is the exception the most signaled by replicas, which is TemporaryTechnicalProblem).

Without the Ehs, the exception signaled by the replica of the airline company could not be trapped and logged. Also, the agent travel agency could not send a response to the broker. Unknown behavior in case of exceptions generates thus a failure of the system, while it was possible to avoid this situation.

Related Work

Multi-agent systems are inherently concurrent, because of the autonomy of agents [START_REF] Wooldridge | An Introduction to MultiAgent Systems -Second Edition[END_REF]. Specific EHS have been designed to handle these concurrency situations when multiple exceptions are asynchronously signaled [START_REF] Campbell | Error recovery in asynchronous systems[END_REF][START_REF] Keen | Exception handling during asynchronous method invocation[END_REF][START_REF] Carlsson | Erlang: Exception handling revisited[END_REF][START_REF] Iliasov | Exception handling in coordination-based mobile environments[END_REF][START_REF] Miller | The guardian model and primitives for exception handling in distributed systems[END_REF][START_REF] Mostinckx | Ambientoriented exception handling[END_REF][START_REF] Platon | A definition of exceptions in agent-oriented computing[END_REF][START_REF] Cacho | Exception handling in context-aware agent systems: A case study[END_REF]. Most of these EHS are based on a dedicated entity (called a monitor, guardian, action, ...) which monitors a group of active concurrent entities and handles the exceptions signaled by any member of the group This kind of EHS is designed to handle exceptions that impact collectively a group of entities which share an execution context. Fewer works study the integration of EHS to asynchronous service (method) invocation schemes. This kind of EHS are best suitable for managing applicative (business) exceptions pertaining to peer-to-peer collaborations. However, they are basic EHSs provided by middlewares which lack concurrent exception resolution and propagation. Our EHS can be considered as an original hybrid of these two approaches, able to seamlessly manage the execution of complex collaborations spawning over groups of agents, while preserving agent behavior encapsulation and reactivity thanks to the support for contextualized exceptions handlers defined as part of the applicative behavior of agents.

The work presented in this paper is an extension of Sage which addresses replicated agents. Replication is a special case of n-versionning, which consists in executing multiple versions of a given functionality for robustness and faulttolerance. All the results returned by the multiple versions are collected and compared. The effective result of the execution is resolved as a majority vote. [START_REF] Mancini | Exception handling in replicated systems with voting[END_REF][START_REF] Issarny | An exception handling model for parallel programming and its verification[END_REF][START_REF] Romanovsky | An exception handling framework for n-version programming in object-oriented systems[END_REF] have studied exception handling in n-versionning frameworks. A main difference in our work is that the responses returned by a group of agent replicas are not resolved by votes. As all the replicas are identical in a group, they should return the same responses, as far as their executions are predictable (are not interrupted by contextual exceptional situations). The responses of the group leader are thus returned primarily. When the leader undergoes a failure (returns a contextual runtime exception), the other replicas are used as failovers. The first replica which succeeds in returning a valid response (whether a result or an expected business exception) becomes the leader. This resolution strategy promotes execution performance. A perspective is to compare all the valid responses of the replicas in order to detect byzantine faults [START_REF] Kotla | Zyzzyva: Speculative byzantine fault tolerance[END_REF] and also address robustness issues.

Conclusion and Discussion

In this paper, we have presented an exception handling system combined with a replicated agent environment for improving software reliability. This EHS firstly offers agent programmers an exception signaling and handling mechanism that works transparently with replicated agents. It second provides to the replication system programmers the capability to control internal exceptions raised by the replication algorithms, as proposed in [START_REF] Mancini | Exception handling in replicated systems with voting[END_REF]. This last point is not developed in this paper. We have proposed an original and lightweight programming Api, using java annotations to define handlers and resolution functions. We detailed then a handler search and a handler invocation algorithms that relies on service execution traces. This works asynchronously to improve agent reactivity. All the propositions have been implemented on the DimaX platform and instrumented through a case study.

As a future work, we wish to further analyze the interactions between the replication mechanism and the exception handling system in order to refine replication strategies (to deal with passive and semi-active replication). Besides, we will consider comparing the responses (normal responses or exceptions) to a same request given by several replicas of an agent in order to learn some agent's execution behavior. Indeed, if an exception is repeatedly raised by replicas of a given agent, it might probably indicate that the exception is a business exception. In such a situation, after a determined number of occurrences of the same exception, the system might consider it is useless to wait for other replicas' responses. This allows us to dynamically distinguish between replica-specific and replica-independent exceptions, without typing statically exception classes.

Fig. 1 .

 1 Fig. 1. An e-commerce application example with 3 agents : contractual requests and responses.

Fig. 2 .

 2 Fig. 2. An e-commerce application example with 3 agents : contractual responses.

Fig. 4 .

 4 Fig. 4. Examples of request, service, agent-level handlers and resolution function definitions in X-SaGE

Fig. 6 .

 6 Fig. 6. Flowchart for handler search

Fig. 7 .

 7 Fig. 7. EHS integrated to DimaX

The Facoma project (http://facoma.lip6.fr/) is partly financed by the french national research agency (Anr).

Further details on the implementation are available at the Url: http://www.lirmm.fr/∼tibermacin/facoma/DimaxX/

We plan at last to use replication as a support to give the core implementation of an EHS that supports a resumption policy. Indeed, even if handler search is stack destructive, as in most systems, a replica of an agent could restart the computation where it has been stopped in the original one.

Acknowledgments. Authors wish to thank the French Research Agency (Anr) that supported this work through the Facoma project of the SetIn 2006 program. They also want to thank A. Romanovsky for his help and all colleagues from the Facoma project -J.-P. Briot, O. Marin and J.-F. Perrot -for fruitful and inspiring discussions.


Michel Fliess
email: michel.fliess@polytechnique.edu

Cédric Join
email: cedric.join@cran.uhp-nancy.fr

Samer Riachy
email: samer.riachy@ensea.fr

Revisiting some practical issues in the implementation of model-free control

Keywords: Model-free control, intelligent PID controllers, estimation, identification

This paper simplifies several aspects of the practical implementation of the newly introduced model-free control and of the corresponding intelligent PID controllers (M. Fliess, C. Join, "Model-free control and intelligent PID controllers: towards a possible trivialization of nonlinear control?," 15 th IFAC Symp. System Identif, Saint-Malo, 2009). Four examples with their computer simulations permit to test our techniques.

INTRODUCTION

Let us start with a brief review of the general principles of model-free control, introduced by Fliess & Join [2008[START_REF] Lévine | Analysis and Control of Nonlinear Systems[END_REF], and of the corresponding intelligent PID controllers (see d 'Andréa-Novel et al. [2010b] for the connections with "classic" PIDs), which already led to a number of exciting applications in various fields: d' Andréa-Novel et al. [2010a], [START_REF] Choi | Model-free control of automotive engine and brake for Stop-and-Go scenarios[END_REF], [START_REF] Formentin | Model-free control for active braking systems in sport motorcycles[END_REF], [START_REF] Gédouin | Experimental comparison of classical pid and model-free control: position control of a shape memory alloy active spring[END_REF], [START_REF] Join | Étude préliminaire d'une commande sans modèle pour papillon de moteur[END_REF][START_REF] Join | Vers une commande sans modèle pour aménagements hydroélectriques en cascade. 6 e[END_REF], [START_REF] Michel | Modelfree control of dc/dc converters[END_REF], [START_REF] Villagra | Robust stop-and-go control strategy: an algebraic approach for nonlinear estimation and control[END_REF][START_REF] Villagra | Control basado en PID inteligentes: aplicación al control de crucero de un vehículo a bajas velocidades[END_REF][START_REF] Villagra | A model-free approach for accurate joint motion control in humanoid locomotion[END_REF]. For simplicity's sake we are restricting ourselves to singleinput single-output systems. The input-output behavior of the plant is assumed to be well approximated within its operating range by an ordinary differential equation E(y, ẏ, • • • , y (a) , u, u, • • • , u (b) ) = 0, which is nonlinear in general and unknown, or at least poorly known. Replace it by the "ultra-local" model, which is continuously updated,

y (n) (t) = F (t) + αu(t) (1)
where

• the order n, 1 ≤ n ≤ a, of derivation has always been chosen to be equal to 1 or 2, and 1 in all concrete situations, 1 • the constant coefficient α is chosen by the practitioner, such that αu and y (n) are of the same order of magnitude, • the time-varying function F (t), which is estimated thanks to the knowledge of u and y, subsumes the structural properties of the unknown system.

Close the loop, if n = 2 in Equation ( 1), with an intelligent PID controller, or i-PID,

1 See [START_REF] Fliess | Rien de plus pratique qu'une bonne théorie: la commande sans modèle[END_REF] for an explanation.

u = - F + ÿ + K P e + K I e + K D ė α (2)
where

• y is the output reference trajectory,

• e = y -y is the tracking error,

• K P , K I and K D are the usual tuning gains.

The above control strategy was put into practice until now via an estimate of the n th -order derivative y (n) in Equation (1), which yields an estimate of F in Equations ( 1) and (2).

In spite of recent significant advances on the numerical differentiation of noisy signals by [START_REF] Mboup | Numerical differentiation with annihiators in noisy environment[END_REF], this task remains quite complex and sometimes difficult to implement. We replace it by recent, but quite simple, algebraic and non-asymptotic techniques [START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008]) for online parameter identification of linear systems. They have been utilized in several concrete case-studies (see, e.g., [START_REF] Abouaïssa | Fast parametric estimation for macroscopic traffic flow model[END_REF], [START_REF] Becedas | Adaptive controller for single-link flexible manipulators based on algebraic identification and generalized proportional integral control[END_REF], [START_REF] Pereira | Adaptive input shaping for manoeuvring flexible structures using an algebraic identification technique[END_REF][START_REF] Sira-Ramírez | A fast on-line frequency estimator of lightly damped ibrations in flexible structures[END_REF]). Let us summarize this new viewpoint by considering the equation

L d dt z = φ + αu (3)
where

• φ ∈ R is an unknown constant, • d ν z dt ν = y, for some ν ≥ 0, • L( d dt ) ∈ R[ d dt
] is a linear differential operator with constant coefficients.

φ is linearly identifiable according to [START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008]. We thus approximate an unknown function like F by a piecewise constant one. We show that this new setting, which is easier to grasp and to implement, possesses excellent robustness properties.

Another important modification with respect to Fliess & Join [2008[START_REF] Lévine | Analysis and Control of Nonlinear Systems[END_REF] is related to the case of a partially known model. If this partially known model happens to be flat [START_REF] Fliess | Flatness and defect of non-linear systems: introductory theory and examples[END_REF], [START_REF] Lévine | Analysis and Control of Nonlinear Systems[END_REF], [START_REF] Sira-Ramírez | A fast on-line frequency estimator of lightly damped ibrations in flexible structures[END_REF]), it might facilitate the choice of a reference trajectory and of a corresponding nominal control. The stabilization around this reference is nevertheless achieved in Section 3.4 in a straightforward way, thanks to the model-free i-PID (2).

Section 2 explains our identification procedure. Four examples accompanied by several computer simulations are discussed in Section 3 in order to test our implementation, even with quite noisy measurements. Some concluding remarks are presented in Section 4

ONLINE PARAMETER IDENTIFICATION

Linear identifiability

Rewrite Equation (3) via the classic rules of operational calculus (see, e.g., [START_REF] Yosida | Operational Calculus[END_REF])

L(s)Z = φ s + αU + I(s) (4)
where I ∈ R[s] is a polynomial associated to the initial conditions. For N ≥ 1 sufficiently large, d N I ds N ≡ 0. Multiplying both sides of Equation ( 4) by d N ds N permits to get rid of the initial conditions. It yields the linear identifiability [START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008]) of φ thanks to the formula

d N ds N L(s)Z = (-1) N N ! s N +1 φ + α d N U ds N (5)
Multiplying both sides of Equation ( 5) by s -M , where M ≥ 0 is sufficiently large, permits to get rid of positive powers of s, i.e., of derivatives with respect to time. The corresponding formulae in the time domain are easily deduced thanks to the correspondence between d ds and the multiplication by -t in the time domain.

Identification scheme

Assume that n = 1 in Equation (1). Close the loop, like in d' Andréa-Novel et al. [2010a], via an i-P, i.e., an i-PID (2) where K I = K D = 0. Assume that F is approximatively constant during the short time window [T -δ, T ]. The above algebraic manipulations lead to the following estimate of F :

F ≈ 1 δ T T -δ ẏ - α δ T T -δ u - K P δ T T -δ e. (6)
The robustness with respect to noises is ensured by the integrals which are the simplest low-pass filters. 2 Note moreover that the estimator (6) may be easily implemented in the form of a discrete linear filter.

FOUR COMPUTER SIMULATIONS

A perturbed pendulum

Bring and maintain an actuated simple pendulum around its upright unstable equilibrium position. Numerical simulations are performed via 2 See the explanations in [START_REF] Fliess | Analyse non standard du bruit[END_REF].

J θ = mgl sin θ + τ -ksign( θ) -c θ (7) where θ is the angular coordinate, and m = 0.2, g = 10, l = 0.7, k = 0.01, c = 0.4 and J = ml 2 are physical parameters. The control torque is τ . Take n = 1 in Equation ( 1) and estimate F by a procedure similar to the one in Section 2.2. Figures 1 and2 show a successful test where the pendulum starts near to its downward position and rejoin its upright position. Noises are added in Figures 3 and4. The control input, displayed in Figure 4, reflects the excellent filtering provided by the estimator of F especially when the pendulum is close to its upright (zero) position where the signal to noise ratio is very small (negative in dB). The control input is actualized each 10 -3 second and a Runge-Kutta method is used in order to simulate Equation ( 7) during each time increment of 10 -3 second. 1). The model-free controller provides very good tracking performances in the absence of noise (see Figures 5 and6), as well as in a noisy situation as displayed in Figures 7 and8.

A comparison with a sliding mode controller

Control via sliding modes of perturbed nonlinear, systems where the classical matching condition of [START_REF] Drazenovic | The invariance conditions in variable structure systems[END_REF] is not satisfied, is a challenging problem. It has been the subject of many studies especially by [START_REF] Estrada | Quasi-continuous HOSM control for systems with unmatched perturbations Int. Workshop Variable Structure Systems[END_REF], 2010], where a particular class of nonlinear single-input systems of the form

ẋ1 = f 1 (x 1 , t) + B 1 (x 1 , t)x 2 + ω 1 (x 1 , t) ẋi = f i (x i , t) + B i (x i , t)x i+1 + ω i (x i , t), i = 2, • • • , m -1 ẋm = f m (x, t) + B m (x, t)u + ω m (x, t)
= (x 1 , • • • , x i ).
The scalar functions f i and B i smooth. The unknown perturbations ω i are bounded and at least n -i differentiable. The controller is a cascaded structure where the state x i is controlled through the virtual input x i+1 . Since the states x 2 , • • • , x m , taken as virtual inputs, cannot be discontinuous, higher order sliding mode algorithms were used for each

x i , i = 2 • • • , m -1.
We consider the simulation examples treated in [START_REF] Estrada | Quasi-continuous HOSM control for systems with unmatched perturbations Int. Workshop Variable Structure Systems[END_REF], 2010] and test our model-free controller. We start first by the example treated in [START_REF] Estrada | Quasi-continuous HOSM Control for systems with unmatched perturbations[END_REF] which is given by ẋ1 = 2 sin(x 1 ) + 1.5x 2 + g 1 (x 1 , t)

ẋ2 = 0.8x 1 x 2 + x 3 + g 2 (x 1 , x 2 , t) (8) ẋ3 = -x 2 3 + 2u + g 3 (x 1 , x 2 , x 3 , t) g 1 (x 1 , t) = 0.2 sin(t) + 0.1x 1 + 0.12 g 2 (x 1 , x 2 , t) = 0.3 sin(2t) + 0.2x 1 + 0.2x 2 -0.4 g 3 (x 1 , x 2 , x 3 , t) = 0.2 sin(2t) + 0.2x 1 + 0.3x 2 + 0.2x 3 + 0.3
The goal is to track y (t) = 2 sin(0.15t) + 4 cos(0.1t) -4. The simulations of [START_REF] Estrada | Quasi-continuous HOSM Control for systems with unmatched perturbations[END_REF] are reproduced in Figures 9 and10. For the model-free controller, n = 2 is used in Equation ( 1) to design an i-PD. Estimation of F is computed by adapting the method described in Section 2. The initial conditions are (0.2, 0, 0) T . For the sake of comparisons with the simulations in [START_REF] Estrada | Quasi-continuous HOSM Control for systems with unmatched perturbations[END_REF], noise-free simulations are also displayed.

The output tracking, which is quite similar to the sliding mode one shown in Figure 9, is not shown again. Besides, the control input in Figure 11 is smoother than the one in 10. It seems that our controller provide two advantages with respect to the sliding mode based controller:

• There is only one parameter, α in Equation ( 1), to be tuned. The choice of K P and K D is trivial. • Smoothness of the control input, i.e., no chattering. Remark 1. The chattering in the controller of [START_REF] Estrada | Quasi-continuous HOSM Control for systems with unmatched perturbations[END_REF] could have been removed by considering a new input u = v. But this would have induced more complications in the controller synthesis. Remark 2. We point out that the local model utilized to synthesize the i-PD is of the second order while the relative degree of ( 8) with respect to e = y -y is three. In [START_REF] Estrada | Quasi-continuous HOSM control for systems with unmatched perturbations Int. Workshop Variable Structure Systems[END_REF], trajectory tracking of the system

ẋ1 = 2x 1 + 1.5x 2 + g 1 (x 1 , t) ẋ2 = x 2 + x 3 + g 2 (x 1 , x 2 , t) (9)
ẋ3 = -1.5x 3 + 2u + g 3 (x 1 , x 2 , x 3 , t) was accomplished with the same sliding mode controller of [START_REF] Estrada | Quasi-continuous HOSM Control for systems with unmatched perturbations[END_REF] but with different parameters.

We tested our model-free controller based on the ultralocal model (1) with n = 2 and exactly the same parameters α, K P and K D previously used in (8). The control input can be seen in Figure 12. Note that the control input in [START_REF] Estrada | Quasi-continuous HOSM control for systems with unmatched perturbations Int. Workshop Variable Structure Systems[END_REF] is similar to the one in Figure 12 but with a chattering of amplitude 1. Lack of space prevents us from reproducing it here. Consider as in Fliess & Join [2008[START_REF] Lévine | Analysis and Control of Nonlinear Systems[END_REF] a nonlinear spring-mass system: 1). This is again achieved via an i-PI strategy. We impose a critically damped behavior by using the characteristic equation s 2 + 2ξω n s + ω 2 n = 0, and choosing ξ = 0.707 and a settling time of 1 second. Only 20 samples are needed to get a good estimate of F in (1). A Runge-Kutta algorithm is used to simulate Equation (10) for each time increment of 0.01 second. The control input is updated at a rate of 0.01 second. The simulation results in the noise-free case are shown in Figures 15 and16

mÿ = -k 1 y -k 3 y 3 + F( ẏ) -d ẏ + u (

CONCLUSION

Our model-free control design, which

• permitted to bypass the difficult task of mathematical modeling, • leads to a straightforward gain tuning, has been improved by • replacing numerical differentiation of noisy signals by simple online parameter identification procedures, • using the model-free i-PID (2) even if a submodel is partially known.

See [START_REF] Fliess | Rien de plus pratique qu'une bonne théorie: la commande sans modèle[END_REF] for an up to date survey of model-free control, the valisity of which has already been confirmed by several concrete applications.

Fig. 3 .Fig. 4 .

 34 Fig. 1. Pendulum angular position: θ

Fig. 5 .

 5 Fig. 5. Output & reference angular velocities.

Fig. 8 .

 8 Fig. 6. Control input for the DC motor

Fig. 10 .

 10 Fig. 9. output & reference

Fig. 12 .

 12 Fig. 11. Control input of the i-PD

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Noisy output & reference trajectory for the springmass system with the nominal control.

Fig. 15 .Fig. 19 .

 1519 Fig. 15. Noise-free output & referene trajectory for the spring-mass system