D M Haughton

R W Ogden

BIFURCATION OF INFLATED CIRCULAR CYLINDERS OF ELASTIC MATERIAL UNDER AXIAL LOADING-II. EXACT THEORY FOR THICK-WALLED TUBES

, (Part I of this series), for membrane cylinders. The relative importance of prismatic, axisymmetric and asymmetric bifurcations for axial tension and compression combined with internal or external pressure is discussed. In particular, the dependence on wall-thickness, ratio of length to radius, and magnitude of the axial stretch is detailed. Some analytical results are presented for a general form of incompressible isotropic elastic strain-energy function. For the most part, however, it is necessary to solve the equations numerically. This is done for specific forms of the strain-energy function.

INTRODUCTION

IN [START_REF] Chadwick | t Designated (I) in the text. 1972 1977 1972[END_REF], (Part I of this series, henceforth referred to as (I)), detailed results were given relating to the bifurcation from a circular cylindrical configuration of an elastic membrane tube subjected to combined axial loading and internal pressure. In the present paper, Part II of this series, an investigation of the corresponding problems for tubes of finite wall thickness is carried out. The results in (I) serve as a guide through the more difficult analysis required for the nonhomogeneous deformations arising in thick-walled tubes.

As in (I), prismatic, axisymmetric and asymmetric modes of bifurcation are examined. A limited number of analytical results are found for an arbitrary form of incompressible,. isotropic elastic strain-energy function, and specific details then discussed for a particular class of strain-energy functions. Analytical results are not obtainable for most bifurcations even for very simple forms of strain-energy function so the governing differential equations are solved numerically. The numerical results are illustrated for four specific forms of strain-energy function.

For thick-walled circular cylindrical tubes the prismatic type of bifurcation has been examined by several authors. In particular, HILL (1976) has obtained analytical results for the "n = 1" bifurcation mode using a method different from that described here. WANG and ERTEPINAR (1972) consider a thick-walled tube subjected to an external pressure and no axial force. They give specific results for the neo-Hookean form of strain-energy function, as also do NowiNSKI and SHAHINPOOR (1969) who confine attention to infinitely-long tubes. However, for externally-pressurized tubes, bifurcation occurs at "small" strains and the results obtained are essentially independent of the precise form of the strain-energy function. We have confirmed this by calculations carried out for several different forms of strain-energy function. Under external pressure it is found, moreover, that a prismatic mode occurs before any possible axisymmetric mode.

The internally pressurized thick-walled tube, which appears not to have been studied previously, provides a contrast to the above results. When the tube is subjected to an internal pressure combined with axial tension the axisymmetric type of bifurcation is found to dominate and the behaviour exhibited is similar in character to that found in (I) for a membrane.

For a tube subjected to internal pressure and axial compression, results are given which, in the limit of zero pressure, coincide with those obtained by WILKES (1955 ), who examined the stability of a thick tube subjected to an end-thrust and gave specific results for the neo-Hookean form of strain-energy function. For the same problem, with attention restricted to the neo-Hookean form of strain-energy function from the outset, SIERAKOWSKI, SuN and EBCIOGLU (1975) have similarly investigated axisymmetric modes of deformation. PATTERSON (1976), also for the neo-Hookean form of strain-energy function, obtains numerical results for asymmetric deformation modes when axial compression is combined with external pressure.

Results are given here for the combination of internal pressure and axial loading (tensile and compressive). The relative importance of the prismatic, axisymmetric and asymmetric bifurcations is discussed in relation to (i) the thickness of the cylinder wall, (ii) the ratio of the length to radius of the cylinder, and (iii) the magnitude of the axial stretch.

THE FINITELY-DEFORMED CIRCULAR CYLINDRICAL CoNFIGURATION

We suppose that the thick-walled cylindrical shell is defined by

A~ R ~ B, 0 ~ e ~ 2n, O~Z~L (1)
in its undeformed configuration, where R, 0, Z are cylindrical polar coordinates.

The cylinder is subjected to an axial loading and an internal pressure so that the circular cylindrical shape is maintained, and the shell defined by

a~ r ~ b, 0 ~ e ~ 2n, 0~ z ~ /,
(2) where r, e, z are cylindrical polar coordinates in the current (deformed) configuration.

Under the assumption that the material is incompressible, the deformation is described by the equations

(3)
where Az is the axial extension ratio.

Let e 1 , e 2 , e 3 denote the unit basis vectors corresponding to the coordinates e, z, r respectively and let A. 1 , A. 2 , A. 3 be the corresponding principal stretches.

From the incompressibility constraint we have 212223 = 1, while from (3 ), where we have used (4) and introduced the azimuthal stretch 2.

We regard Az as a given constant and 2 as a function of r (orR) from (3).

(4)

(5)

For an incompressible isotropic elastic solid the strain-energy function

W = W(2~>
(i = 1' 2, 3 ), (6)
where aw CT; =A; a;: (i= 1,2,3), ( 7)

'
and p is the arbitrary hydrostatic pressure arising from the incompressibility constraint.

Use of ( 5)-( 7) leads to

and W;_, J.li denote aw;a;., oWfi72z.

Equati~ns (8) are similar in form to those given in (I) for a membrane except that a 33 # 0 here and there is now dependence on the radius r.

For the symmetric configuration considered here, the only equilibrium equation not satisfied trivially is (10) it being assumed that there are no body forces.

If a 33 is -P on r =a and zero on r = b, where P denotes the internal pressure, integration of (10) leads to

_ ~ -dr _ ~• W;. d2 P-J 2W;.-= J 2 . a r Ab (2 Az-1) ( 11)
In deriving (11) use has been made of (8 1 ) and we have adopted the notation

2 0 = ajA, 2b = bjB, which, from (3d, gives 2 R 2 2 2 2a2z-1 = A2 (2 Az-1) = t~(2b2z-1), (12)
where (13) t The subscript i is not summed here or in the sequel.

Note that equality in (13) corresponds to the membrane limit discussed in (I), while from (12) it follows that (14) with equality holding if and only if P = 0. (In this connection we recall from (I) that W;_ ~ 0 according as A 2 Az ~ I.)

The change of independent variable from r to ;, in ( 11 ) is achieved by means of ( 3 1 ) and the definition (5 2 ).

In view of ( 12) we regard P as a function of ),z and Aa• Differentiation of ( 11) with respect to A 0 then gives

ap A A Jc; 1 (A; Az-1) :;-= W;_(Au, Azl/1-a-W;_(A.h, i,z)j),b. UAa (15)
This is a straightforward generalization of a formula given by OGDEN ( 1978, equation (36)) for Az = 1.

Equation ( 15) should be compared with the corresponding equation for a membrane given in (I, equation ( 36)) to which it reduces for "small" (B-A)/A.

Clearly, a necessary condition for pressure turning-points to exist at fixed Az is that ( 16) (17) As shown in (1), the equation ( 17) characterizes the pressure turning-points in the case of a membrane. Analogous considerations in respect of a spherical shell have been discussed by [START_REF] Chadwick | t Designated (I) in the text. 1972 1977 1972[END_REF]. Experimental data indicating the occurrence of a pressure maximum for thick-walled tubes have been obtained by SKALA (1970) and[START_REF] Chadwick | t Designated (I) in the text. 1972 1977 1972[END_REF]LI (1977).

.A.-1 WA is not monotonic in A,
The resultant axial loading N on the ends of the cylinder is given by Use of ( 8) and ( 10) enables this to be expressed as ( 18)

The corresponding expression for a membrane is 2e~ where e = HjR, H being the thickness of the membrane and R its radius in the undeformed configuration. In (I) it was shown to be reasonable for the inequalities ~ A > 0, ~ic > 0 to hold, where A is the azimuthal stretch in the membrane. It is also ;e~sonable 'for these inequalities to hold for the thick-walled specimen, but in general it cannot be expected that the analogous global inequalities oNjo.A.z > 0, DNjaA.. > 0 are valid in view of the complicated form of N.

Since o-33 # 0, in general, for a thick-walled cylinder, the local form of the axial stress obtained from (8 2 ) yields

F(Aa, ~•z): ~,(~•a' Az)-A; 1 P,} F(Ab, Az)-W;,(Ab, ).:),
where F(J., Az) = ;.z-1 o-22 is the axial stress per unit undeformed area of cross-section (nominal stress). Even for these local forms it cannot in general be deduced that oF().a, Az)/oJ.z > 0 or oF(J.a, J.z)/o).a > 0. Moreover, oF().b, Az)foJ.z > 0 at fixed Aa is not a consequence of~'-> 0 since ).b depends on Az through (12).

Since N may be t'r~ated as a function of Aa and Az, one consequence of the above discussion is that inflation is not necessarily accompanied by axial shortening of the tube at fixed N. That is,

lf}~ =_oN ~~N d),a - o)
.a oJc= may be positive for certain values of ).a and A=, in contrast to the situation for a membrane (see ( 1)). Indeed, CHADWICK and HADDON (1972) have provided an example to illustrate that a tube can extend during inflation in the absence of axial loading.

If o-~2 denotes the deviatoric part of the axial stress o-22 and N' its resultant, then it is easily shown, directly from (8), that

The significance of this expression lies in the fact that for a tube with closed ends (so there is a contribution Pna 2 to the axial stress) ( 18) is replaced by N = 3N'/2. We shall return to this point later in relation to the discussion of numerical results (Section 4 ).

THE INCREMENTAL EQUATIONS

Full details of the notation and equations used here are given in [START_REF] Chadwick | t Designated (I) in the text. 1972 1977 1972[END_REF]. Here we summarize briefly the relevant information.

The incremental equilibrium equations (in the absence of body forces) can be written div s 0 = 0, (20) where s 0 is the increment in nominal stress referred to the current configuration (indicated by the subscript zero).

The incremental boundary conditions in respect of a hydrostatic pressure loading are (21) evaluated on the appropriate boundary, where superscript T denotes the transpose, n is the unit normal in the current configuration and l1 = ci 0 , ci being the increment in the deformation gradient tx.

The incremental constitutive law is written So = .::1Bl] + pt]-pfJ, (22) where .'d) is the fourth-order tensor of instantaneous moduli associated with the conjugate variables (s, tx) and fJ denotes the identity tensor. The components of ;::18 on the principal axes of the underlying deformation are given by HAUGHTON and OGDEN (1978, equations ( 84)-( 86)) and therefore not repeated here.

With the basis vectors defined in Section 2 and the choice 8, z, r of coordinates, where the non-zero values of ei. ei.k are precisely those given for ai. ai. k in (I, equation ( 23)). In the present situation, however, r varies through the thickness of the shell.

We write the increment x in the position vector x of a point in the current configuration as (24) Then, the components of '1 on the basis e 1 , e 2 , e 3 are displayed as

[ (u + Vo)/r V= Wo/r \\'= (u0 -v)/r Uz (25)
where the subscripts e, z, r denote partial derivatives.

Since the material is incompressible, we have tr('l) = u.+(u+v0 )/r+w= = 0.

(26)

BIFURCATION ANALYSIS

Prismatic bifurcations

Firstly, we assume that u, v, ware independent of z. Then (22)-( 26) together with the expressions for the components :$ ijkt from HAUGHTON and OGDEN (1978) give .:!&1212woo/r+.:i8~232rw.+.:i83232(rw

•• +w.) = 0, ( 29
)
where the prime denotes djdr.

Equation ( 29) can be solved quite generally for w at any stage of the deformation and, in particular, in the undeformed configuration. The solution has no effect on the shape of the cylinder cross-section and deserves no further consideration in the present context. We therefore set w = 0, as in the membrane case treated in (1).

The boundary conditions (21) on the lateral surfaces of the tube become

u 0 -v+rv,=0 onr=a,b (30) and (31)
In the derivation of ( 30) it has been assumed that ~3131 # 0. In fact, the inequality . .11 3131 > 0, a consequence of the Baker-Ericksen inequalities, is adopted here.

To solve equations (27) and (28) with the incompressibility condition (26) we set (

) (38) 37
Equations equivalent to (34), ( 35) and (36) were given by WANG and ERTEPINAR (1972), amongst others, for the special case of a neo-Hookean solid, and by HILL ( 1976) for a general form of strain-energy function.

For 11 = 0, the circular symmetry is maintained and the analysis is equivalent to that described in Section 2. No further discussion of this is necessary.

For n = 1, equation ( 34) is integrated to give The constant of integration has been set equal to zero to ensure that the boundary conditions ( 35) and ( 36) are compatible. Use of ( 40) in ( 36) for n = 1 shows that each of .( 35) and ( 36) reduces to rf{' + f{ = 0 on r =a, h.

{40)

Use of the formula rX = },(1-}. 2 Az) enables two further integrations to be carried out. These, consecutively, give

and ( 42)

where c 1 and c 2 are constants.

Since the differential equations and boundary conditions are homogeneous, f 1 = const. is a solution. Thus, for prismatic bifurcations to occur, a non-zero f{ is required. Substitution of ( 42) into (40) leads to the hifurcatiol1 criterion

ic - f " WvAdA . _ O -(; 2 A -1 )( W ) 2 -' icb • z A (43)
use having been made of the connection

.J& 3131 = XW;j(A 4 A;-1).

In view of (12), equation ( 43) determines the bifurcation points Aa (dependent on ;.zl. It is applicable equally to internally or externally pressurized cylinders at fixed axial extension.

The special case P = 0 corresponds to A 2 Az = 1, that is WA = 0. Separate analysis of this case directly from the differential equations and boundary conditions (noting that the underlying deformation is now homogeneous) shows that prismatic bifurcation is not possible for any mode number n except in the limit h--+ a.

From (43) a necessary condition for 11 = 1 mode prismatic bifurcation to occur is WH = 0 for some A in ().b, l.a). In the membrane limit, equation ( 43) reduces to a necessary and sufficient condition, namely Wu = 0 for some ). (as shown directly in (I)).

However, arguments presented in (I) suggest that the inequality Wv > 0 (44) holds for all A and ).=, ensuring that prismatic bifurcations are not possible for a membrane. For a thick-walled cylinder, the inequality (44) excludes n = 1 mode prismatic bifurcation but higher-order modes remain possible. We remark that although (44) is physically reasonable for a membrane it could possibly be violated within the bulk of a thick material specimen. Nevertheless, we adopt (44) here.

For mode numbers n ~ 2, WANG and ERTEPINAR (1972) have considered the problem of an externally pressurized thick-walled tube in respect of the neo-Hookean form of strain-energy function, restricting attention to il.z = 1. Under external pressure, the strains at bifurcation are "small" and the results obtainable essentially independent of the precise form of strain-energy function. We have confirmed this by calculations for a number of different strain-energy functions.

We use the class of strain-energy functions described fully in (1). This is such that (45)

In particular, in Table 1, we show our calculations in respect of the three-term form of ( 45) with material constants given by

(1(1 = 1•3, (1(2 = 5•0, (1(3 = -2•0, } (46) llf = 1•491, J.l! = 0•003, 11! = -0•023,
where 11: = J.i,/J.i and J.l is the ground-state shear modulus (see (I)).

In Table 1, values of il.b at bifurcation are given for il.z = 1, 2, 3, 4, 5 and 0•5 for a number of different values of the thickness ratio A/Bin respect of ( 45) with ( 46). For comparison we reproduce the results of ERTEPINAR (1972) for il.z = 1 and the neo-Hookean form of strain-energy function. In all cases it is found to be the n = 2 mode which occurs first. We remark that for A/B = 0•95 the critical values of il.b for il.z ~ 1 correspond to P = 0 to within an accuracy of three decimal places, this being consistent with the situation in the membrane limit.

Under internal pressure our calculations reveal that no prismatic mode of bifurcation is possible for the strain-energy functions considered.

In Fig. 1 the pressure P, as given by ( 11), is plotted as a function of il.a in respect of ( 45) with ( 46) for a number of values of il.z and for A/B = 0•85. Because the value of P increases rapidly with B/A we do not includ~ corresponding curves for smaller values of A/B. Instead, in Fig. 2, the value of A.. at the pressure turning-points is plotted against AjB. For sufficiently large wall thicknesses, pressure turning-points do not occur for all values of Az, but this is not shown in Fig. 2. So far no mention has been made of the boundary conditions on the ends of the cylinder. These can be chosen either to exclude prismatic bifurcations or to be compatible with them. They are more important in the discussion of axisymmetric and asymmetric bifurcations, which we consider next. is then substituted in the differential equation ( 53) in order to obtain fa(r) subject to (54) and ( 55) for different choices of the dimensionless geometric parameters L/B and A/B, recalling (1).

A/B 2 3 4 5 0•5 0•95 0•999 0•999 0•707 0•577 0•500 0•447 1•413 0•9 0•997 0•997 0•705 0•576 0•499 0•446 1•410 0•85 0•994 0•994 0•703 0•574 0•497 0•445 1•406 0•8 0•989 0•990 0•700 0•571 0•495 0•443 1•400 0•7 0•979 0•979 0•693 0•566 0•490 0•438 1•385 0•6 0•965 0•968 0•685 0•559 0•485 0•434 1•369 . Cl.
Results calculated on the basis of (61) show how the length of the cylinder affects the onset of axisymmetric bifurcation into then = 1 mode. Equally, and equivalently, if L is fixed and n varied then results show how a cylinder of fixed length responds to axisymmetric bifurcations in different mode numbers. This dual interpretation has been exploited in (I) and we continue it here. The spectrum of values of a from 0 to oo embraces cylinders of arbitrary length for any mode number.

The differential equation ( 53) cannot in general be solved analytically, even for very simple forms of strain-energy function, except in the case a = 0 when a simple explicit form for f 0 (r) can be written down. For finite L, a """ 0 corresponds to maintenance of the circular cylindrical shape and the solution is u = A/r and w = 0, where A is a constant. In the limit L---> oo, when ex---> 0, the solution (59) reduces to u = (A/r) sin az, w = 0, with a/ remaining finite, and these expressions replace ( 56) and ( 57) for this case.

The solution of ( 53) with a= n/A.zL, together with the boundary conditions ( 54) and ( 55), leads to a bifurcation criterion which is a relation between A.z and )"a (or A.b in view of ( 12)). Numerical calculations have been carried out for a number of values of Az and the ratio L/B. The numerical method is described in Section 4.3.

In Fig. 3(a,b,c) for L/B = 20, 10, 5 respectively, the bifurcation curves corresponding to A/B = 0•85 are shown in the (A.z, A 0 )-plane. Results for the neo-Hookean, Mooney-Rivlin and Varga forms of strain-energy function are compared with those for the three-term strain-energy function (46). Corresponding results are given in Fig. 4(a,b,c) for A/B = 0•5.

In each case, the curve of P = 0 is also shown and attention is restricted to P ~ 0.

Regions of both (axial) tension and compression are included. The curves of pressure turning-points (independent of LjB), calculated by setting oP/i3Aa = 0 in (15), are plotted in Figs 3 and 4 in respect of the three-term strain-energy function ( 46). These show that, dependent on A.., A/Band L/B, bifurcation may occur before or after the pressure maximum during inflation, just as for a membrane (see (I)).

For reasons discussed by [START_REF] Chadwick | t Designated (I) in the text. 1972 1977 1972[END_REF][START_REF] Chadwick | t Designated (I) in the text. 1972 1977 1972[END_REF] the neo-Hookean, Mooney-Rivlin and Varga forms of strain-energy function exclude certain bifurcation features. The remainder of our calculations is therefore based on the three-term strain-energy function (46). We note from Figs 3 and4, however, that there is little difference between predictions of the several strain-energy functions in the region of axial compression {Az < 1 ).

In respect of (46) the bifurcation curves corresponding to axial compression are shown on a larger scale in Fig. 5(a, b) for A/B = 0•85 and 0•5 respectively. In each case, the values LIB= 5, 10, 20 are chosen. The curve of zero resultant axial force (that is, N = 0, where N is given by ( 18)) is included in Fig. 5(a,b) to emphasize the boundary between the regimes of resultant axial tension and compression. We note that as the axial compression is increased in the absence of pressure (that is, along the curve;..;;..,= 1), the value of A.a at bifurcation increases with the value of L/B. At first sight this is at variance with one's intuition that long tubes buckle more ,. easily than shorter ones. However, it should be borne in mind that asymmetric (or bending) modes will have priority over axisymmetric modes for long tubes in compression (this is shown in Section 4.3).

For the neo-Hookean solid, WILKES (1955) considered a tube subject to end-thrust alone. In this special case, A. For each value of A.. shown, there is clearly a "cut-otr' value of the tube thickness beyond which axisymmetric bifurcations do not occur with increasing pressure. For long tubes, as illustrated in Fig. 6(b,c), the character of the results does not vary significantly with thickness except near the cut-off values. Indeed, for thicknesses up to, say, 60% of the critical value (for each A.. and L/B) the results are broadly similar to those found in (I) for a membrane and corresponding to the limit A/B = 1 here. It is an easy matter to find a value of a (if one exists in the range under consideration), and hence b from b 2 = ..t; 1 (B 2 -A 2 )+a 2 , to satisfy the bifurcation criterion as accurately as required. This numerical method has in essence been used by ERTEPINAR (1978) to investigate bifurcation of thick-walled tubes of neo-Hookean material subjected to finite twist combined with external pressure.

The above method was the basis for our calculations (corresponding to m = 0) described in Section 4.2.

Form= 1 and IX given by ( 61), bifurcation curves are plotted in the (A.z, ..1. 0 )-plane in Fig. 8(a,b) for A/B = 0•85, 0•5 respectively. The character of the results is essentially independent of A/Band, in particular, there is no "cut-off' as there is in the case of axisymmetric bifurcations. This, of course, is to be expected since solid cylinders (A/B = 0) are susceptible to bending in compression. Bending can only occur when the dilatational part of the axial loading is compressive, that is when N', given by ( 19), is non-positive. In the membrane case of (I) this reduces to the integrand in (19) being non-positive.

In the region where N' > 0, asymmetric modes do not compete with axisymmetric modes, but where N' ~ 0 whether an axisymmetric or asymmetric mode can occur first during inflation at fixed A.z depends on Az and L/B, just as for a membrane (see ( 1)), and also on the ratio A/B. This can be seen by comparing Fig. 8(a, b) with Fig. 5(a,b) respectively and also with the appropriate parts of Figs 3 and 4 respectively. Note, however, the differences in scale. We do not pursue a discussion of the comparison here since, for a membrane, this has been done in (I) and the details are broadly similar for thick-walled tubes provided A/B is not too small.

For m ~ 2, and IX again defined by (61 ), our calculations reveal no bifurcation curves in the region of the (A.., ..1. 0 )-plane considered in respect of the strain-energy function (46).

  and for (16) to hold it is necessary that .A.Wu-W;.=O for some A in the interval (.A.b, A 0 ).

  Pr = [r( :!&~ 33 3-::$'1133 + p') + :18 3 3 3 3 + .:1811 1 1 -2.!-d 1133]u./r + (.18 333 3-.18 11 3 3 )u •• +.181313(Uoo-Vo)/r 2 +.181331 v.o/r, (27) Po= [r( . .$'1331 +p')+.:16t313](uo-v)/r+(r:i8~131 +.183t3tlv. +.:183131 rv •• + (.:181331 +•:iltt33 -.!.d 1111 )u.o, (28)

v

  = g.(r) sin n8, (32) u = J.(r) cos n8,} p = h.(r) cos n8, for example. (In general, four independent functions of {) are associated with each function of ron the right side of (32).) Equation (26) gives rf~ +f.+ ng. = 0, n ~ 0. This is used to eliminate h. and g. from the other equations, leading to r{.:fl3131 r 3 f~" + (r::fl~ 131 + 3..113131 )r 2 f~' (33) + [r.JI~ 131-.fd 3131 + n 2 (2:?if 1331 + 2.16' 1133-.~ 1111-•~ 3333)]rj~}' +(n 2 -1)[r 2 .fl3 131 +r.~~131 +(n 2 -1).~3 131 +n 2 (a 1 -a 3 )]J. = 0 (34) for n ~ 1. Correspondingly, the boundary conditions (30) and (31) for n ~ 1 yield r 2 f~'+rf~ +(n 2 -1)f. = 0 on r =a, b (35) and, after some rearrangements, •~ 3131 r 3 f~" + (r.Jd~ 131 + 4.$3131 )r 2 f:' + + [r . .f.9~131-(n 2 -1 ).~ 3131 + 11 2 (2.:?if 1331 + 2~ 1133 -.~ 1111 -.~3333)]rJ;; +(n 2 -1)(r.~~131 +:14 3131 )f.=0 onr=a,b. (36) For n = 0, (31) gives (.143333-.~1133+a3)fo/r+ho = {~ where h 0 is obtained from on r =a,} on r = b, r 2 h~ = (r.Jd'1133-r.Jl~333rp' +.Jl3333 -:;J 1111 )fo, and (30) is satisfied identically.

  FIG.I. Plot of the pressure P* = P/J-l against ,t. for i.z = I, 2, 3, 4, 5 with A/B = 0•85.

4. 2

 2 Axisymmetric bifurcar;ons When equations (23) are independent of 0 the equation for v decouples from those for u and w. The solutions of this equation are of little interest in the present context and it is considered no further. The remaining equations become Pr = (r:J6'1133-.1H1111 )u/r 2 + (r.qa~333 + rp' + .16 3333)u,jr +.16 3333 u,, + .31 :nn u.z + (r-16'~233 + .::?6'nn-.Jl1122)wz/r + (:JiS' 22H + .16 3223)w"' (47)

  FIG. 3. Plot of the axisymmetric bifurcation curves in the (A=• A0 )-plane for P ~ 0 and AlB= 0•85 in respect of the Mooney-Rivlin (-• --), neo-Hookean (-•-•-•-),Varga(--) and three-term(---) strainenergy functions: (a) LIB= 20, (b) L/B = 10, (c) LIB= 5. The curve of pressure turning-points (calculated for fixed A=) is included in (a), (b) and (c) for the three-term strain--energy function.

  2 A.. = 1 so the underlying deformation is homogeneous and the equations can be solved in terms of modified Bessel functions. Wilkes (in different notation) plots A.. against et.b = mr:BfLA.; 12 and the resulting curve provides a critical value of A. •. Bifurcation into an axisymmetric mode in compression cannot occur for values of A.. greater than this critical value. Whether or not bifurcation occurs at the critical value depends on whether or not the value of nB/L determined by the corresponding critical value of et.b is an integral multiple of the ratio B/L of the considered tube. The points of the curve P = 0 in Fig.5(a, b) corresponding to bifurcation in essence show the behaviour of A.. as a function of nB/L.In Fig.6(a, b,c), for L/B = 5, 10, 20 respectively, the value of A.a at bifurcation is plotted against A/B for a number of values of A.. to illustrate the dependence of bifurcation on the tube thickness.

  FIG. 6. Plot of the value of Au at bifurcation against A/B fori.,= I, 2, 3, 4, 5 in respect of the three-term strain-energy function: (a) L/B = 5, (b) L/B = 10, (c) LIB= 20.

FIG. 8 .

 8 FIG. 8. Plot of the asymmetric bifurcation curves in the (Jc, A0 )-plane for P::;,. 0 and L/B = s. 10. 20. oo in respect of the three-term strain-energy function: (a) A/B = 0•85. (b) A/B = 0•5.

  .1& 3131 r 3 Jl" + (r.J&3131 + 31& 3131 )r 2 f{' + (rd&3131 -.;& 3131 + 21& 1331 + 2.1& 1133-.1& 1 11 1 -.J& 3333 )rf{ = 0. (39)

TABLE 1 .

 1 Values of il.b at bifurcation for various values of il.z and AjB in respect of(45) with(46)

	Az neo-Hookean
	(Ertepinar,	Three-term strain-energy function
	1972)

ACKNOWLEDGEMENT

The work of D.M.H. was supported by a Research Grant from the U.K. Science Research Council.

Pz = ;}l3232 W,, + (r:?~:3232 + .1/J 3232)w,/r+.J12222 Wzz + (:?12233 +.?13223)u,z + (r'~3223 + rp' + :?13223 + .?ltt22)uzlr, (48) and (26) reduces to u,+ u/r+ Wz = 0.

(49)

The boundary conditions (21) on the lateral surface yield (:~3333-.~2233+u3)u,+(.'?III33-.JiJ2233)u/r-p = 0 on r =a, b

and w,+uz=O onr=a,b.

(51)

The equations and boundary conditions have been made homogeneous by the omission of the incremental pressure P, this being associated with a circular cylindrical solution as discussed in Section 2. In order to solve the above system, we write

for example (~ i= 0).

Elimination of w and p between (47), ( 48) and (49) leads to r 4 (:.1l3232f"' + (r.Jl3232 + 2:Jl3232)f''/r+ (r.

To fit the boundary conditions on the ends of the cylinder we write

where A and Bare constants. The corresponding form of w, obtained from (49) and consistent with (51), is

where C 0 (r) is a function of r related to f 0 (r), while p, if required, can be obtained from (

Note that the solution corresponding to ()( = 0 is ruled out by ( 58) and we can take ()( > 0 without loss of generality. The boundary conditions (58) correspond to zero shear stress and zero axial displacement pointwise on the ends of the tube. The values of()(/ given by ( 60) are also obtained if the end conditions (58) are replaced by ones in which the (nominal) stress is held fixed pointwise and is in the axial direction. If the resultant axial loading on the end faces is held fixed, on the other hand, values of()(/ depending on .A.z and .A.a may be obtained. Such dependence arises in the membrane case, treated in (1), in which integration with respect to r through the thickness of the material is implicit from the outset, and the solution analogous to (56 1 ) involves two (in general, distinct) values of()(. Note that (56) is not the most general solution for u since (53) is quadratic in ()( 2 . For a thick-walled tube the use of "average" end conditions involves lengthy algebra, and we therefore confine attention to pointwise end-conditions which, although yielding simple and explicit values of()(/, nevertheless allow the main features of the bifurcation phenomena to be revealed. By the same token it suffices to consider only (58) and (60).

Because the length of the cylinder is at our disposal, we can arbitrarily set n = 1 in For a number of values of).., the value of ).a at bifurcation is plotted against B/L in Fig. 7, the left-hand limit corresponding to an infinitely-long tube. The ratio A/B = 0•85 has been chosen here in view of the results in Fig. 6(a,b,c) which indicate that only for short tubes (that is, for B/L larger than about 0•15) is there a significant dependence on A/B. Typical results for short tubes can be read off from Fig. 6(a).

The duality between the mode number n and the ratio BjL shows that for fixed ).. (or B/L) there is a value of B/L (or )..) greater than which axisymmetric bifurcation cannot occur in any mode number.

Asymmetric bifurcations

When dependence on both () and z is retained, equations ( 23 The incompressibility condition ( 26) is also required. The equations are solved by setting u = f(r) cos me sin rxz,} v = g(r) sin me sin rxz, w = h(r) cos me cos rxz' p = k(r) cos me sin rxz, (65) the z-dependence selected, as in ( 59), in order to satisfy end-conditions which again lead to values of rx given by (60). In order that the solutions be single-valued we take m = 0, 1, 2, ....

The expressions (65) are substituted in ( 62)-( 64) and h(r) is eliminated by means of ( 26). The resulting coupled equations for f(r), g(r) and k(r) are then put in the form dyj .

where y = (f,f' ,J", g, g'' k),

suitably non-dimensionalized. The corresponding boundary conditions, after some algebra, can be put in the form rys-Y4-myl = 0, } r 2 YJ + rYl + (rx 2 r 2 + m 2 -1 )y1 = 0, (:?11133-.?!3'2233)(y 1 + my4) + (.c?!S'3333-.?!3'2233 +a 3)ry2-ry6 = 0, each for r = a, b. the vanishing of the resulting 6 x 6 determinant of coefficients of ck.


Luca Bortolussi
email: luca.bortolussi@gmail.com

Nicolas Gast
email: nicolas.gast@inria.fr

Mean Field Approximation of Uncertain Stochastic Models

Keywords: stochastic models, population, parameter estimation, mean-field approximation, differential inclusions

We consider stochastic models in presence of uncertainty, originating from lack of knowledge of parameters or by unpredictable effects of the environment. We focus on population processes, encompassing a large class of systems, from queueing networks to epidemic spreading. We set up a formal framework for imprecise stochastic processes, where some parameters are allowed to vary in time within a given domain, but with no further constraint. We then consider the limit behaviour of these systems as the population size goes to infinity. We prove that this limit is given by a differential inclusion that can be constructed from the (imprecise) drift. We provide results both for the transient and the steady state behaviour. Finally, we discuss different approaches to compute bounds of the soobtained differential inclusions, proposing an effective controltheoretic method based on Pontryagin principle for transient bounds. This provides an efficient approach for the analysis and design of large-scale uncertain and imprecise stochastic models. The theoretical results are accompanied by an in-depth analysis of an epidemic model and a queueing network. These examples demonstrate the applicability of the numerical methods and the tightness of the approximation.

I. INTRODUCTION

The complexity of the world prevents us to describe all its aspects with full precision. Lack of knowledge and limited computational and intellectual resources force any description of a complex behaviour to be imprecise and uncertain to a certain degree. This is true also when we try to construct mathematical models of complex systems, in domains as diverse as telecommunication, molecular biology, epidemiology. Probabilistic models are one way of representing uncertainties but they depend on parameters which values are never known precisely, and can possibly be estimated from data only with a certain degree of imprecision. Furthermore, mathematical models of systems cannot provide an explicit description of the environment, which always influences the behaviour in uncontrollable ways.

As a typical example, consider a model of epidemic spreading in a population, for instance a model of disease spreading in humans [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF] or of a malware spreading in a computer network [START_REF] Stamatiou | Computer Network Epidemics: MODELS and Techniques for Invasion and Defense[END_REF]. The infection rate is a typical parameter of such models which can hardly be known exactly, and is usually estimated from available data about the initial outbreak of the epidemic [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF]. However, statistical estimation can never provide an exact value, as we have at disposal only a finite, and often insufficient, amount of data. Moreover, the infection rate itself may depend on environmental factors, which can change arbitrarily in time. For instance, in cholera spreading [START_REF] Mari | Modelling cholera epidemics: The role of waterways, human mobility and sanitation[END_REF], the level of rainfall impacts on the diffusion of the bacterium among nearby water reservoirs, leading to infection rates which can vary unpredictably in time.

Another example is provided by sociotechnical systems such as bike sharing systems. In such as system, users can pick up or return bikes at any station of their choice. The arrival rate of users in a bike station cannot be assumed fixed during the day, but it will depend on the current hour. The precise form of such a dependency is unknown, as it will be influenced by several factors like the weather, the status of public transportation, the presence of events in the city. Even if we restrict to a small time frame, like the rush hour in the morning, so that we can assume a constant arrival rate, such a rate cannot be reliably fixed to a precise value, but should rather be assumed to lie in an interval of possible values.

In this paper, we distinguish two models of uncertainties Imprecision: Some parameters ϑ can depend on features of the environment external to the model (like temperature, atmospheric weather, and so on,). ϑ may be subject to variations during the time horizon T of interest, so that considering it as a fixed parameter is an incorrect assumption that can lead to incorrect results. One way to capture such a variability, without committing to assumptions on the form of dependency of ϑ on the external / environmental factors, is to fix a set Θ of possible values for ϑ and assume that ϑ depends on time t and can take any value of Θ at any time instant, i.e. ϑ t ∈ Θ. Essentially, ϑ t plays a role similar to inputs in control theory. We call this the imprecise scenario. Uncertainty: In a simpler scenario, a parameter ϑ is assumed fixed, but its precise value is not known precisely. In this case, we assume that ϑ ∈ Θ, where Θ is the possible set of values of ϑ, as above. This will be referred to as the uncertain scenario.

In this paper, we focus our attention on models of large populations of interacting agents, formalised in terms of continuous time Markov chains (CTMC), introducing a new class of models, called imprecise population processes, which are able to capture the intrinsic uncertainty of the world. In these models, we assume that some parameters can vary in an unconstrained way within a certain range. This non-deterministic variation represents all possible ways the external environment can influence the value of these parameters, and is a model of the imprecise scenario described above. We also distinguish a simpler case, in which parameters are assumed to be fixed, i.e. not influenced by the environment, but unknown. In this case, the class of population models so obtained is considerably simpler, resulting in the so-called uncertain continuous time Markov chains, see [START_REF] Bortolussi | Smoothed Model Checking for Uncertain Continuous Time Markov Chains[END_REF].

When populations are large, as is often the case in epidemiology, biology, telecommunications, a direct analysis of the stochastic model is unfeasible, even by restricting to a statistical treatment of simulations. A viable alternative is to construct the so-called mean field approximation [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF], which provides a deterministic description of the behaviour to which stochastic trajectories converge in the limit of infinite populations.

The complexity of the analysis is further exacerbated by the presence of imprecision or uncertainty. The main contribution of this paper is to provide a characterisation of mean-field limits for imprecise population processes (and a-fortiori for uncertain CTMCs) in terms of differential inclusions (DI, [START_REF] Aubin | Differential Inclusions[END_REF]), both for transient and steady-state behaviour (when this is meaningful). We investigate the computational gains obtained in this way. With respect to classic mean field, the presence of imprecision increases the computational cost of analysis also for the mean-field limits, as we have to deal with differential inclusions. In the paper, we also present a novel method to bound the solutions of DI using control theoretic tools, namely the Pontraygin principle [7, Section 3].

These results pave the way to novel and efficient computational methods to analyse, design, and eventually control large scale systems of interacting agents, taking both uncertainty in parameter specification and imprecision due to external effects consistently into account. For instance, we can design patching (or vaccination) strategy to counteract an epidemic which is effective even if the infection rate changes in time in unpredictable ways.

The paper is organised as follows: Section II discusses Imprecise Markov Chains and the evolution of their probability mass in terms of differential inclusions. Section III introduces imprecise population models and proves their mean field behaviour. Section IV discusses computational methods for the analysis of the differential inclusion mean field limits. Section V works out in detail an example about epidemic spreading, while Section VI discusses a generalised processor sharing scenario. The final discussion is in Section VII.

Related work: Following the terminology of [START_REF] Aubin | The viability theorem for stochastic differential inclusions 2[END_REF], the models developed in our paper combine two ways for representing uncertainties: a "stochastic uncertainty", driven by a Markovian behaviour, and a "contingent uncertainty", given by possibly changing parameters. This approach is similar to the notion on stochastic differential inclusion studied in [START_REF] Aubin | The viability theorem for stochastic differential inclusions 2[END_REF], [START_REF] Kisielewicz | Stochastic differential inclusions and applications[END_REF], where the "stochastic uncertainty" is driven by a Wiener process while the "contingent uncertainty" is given by a setvalued map. Our work also builds on [START_REF] Škulj | Discrete time markov chains with interval probabilities[END_REF], that considers Markov chain with interval probabilities. The key contribution of our paper with respect to [START_REF] Škulj | Discrete time markov chains with interval probabilities[END_REF] is to extend this notions to population process and to develop a rigorous mean-field approximation of such systems. Mean field approximation of population processes has a long history, starting from the works of Kurtz [START_REF] Kurtz | Solutions of Ordinary Differential Equations as Limits of Pure Jump Markov Processes[END_REF], [START_REF] Kurtz | Markov Processes -Characterisation and Convergence[END_REF]. Mean field based analysis, for the transient and the steady state, have been applied in performance modelling and model-checking tools [START_REF] Tschaikowski | Approximate reduction of heterogenous nonlinear models with differential hulls[END_REF], [START_REF] Kolesnichenko | A logic for model-checking mean-field models[END_REF], systems biology, epidemiology [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF]. For a gentle introduction, see [START_REF] Bortolussi | Continuous approximation of collective systems behaviour: A tutorial[END_REF]. Classic mean field results require (Lipschitz) continuity of rate functions, but more general theorems can be proved for piecewise smooth rates [START_REF] Bortolussi | Hybrid Limits of Continuous Time Markov Chains[END_REF] or even general discontinuous functions by using differential inclusions [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF]. The proof of most of our convergence results are based on the construction of a proper stochastic approximation. They are similar to the ones of [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF], [START_REF] Roth | Stochastic approximations with constant step size and differential inclusions[END_REF]. For clarity and brevity, we will only give sketches of proofs.

Note that the main technical contribution of this paper is to properly define the mean-field framework for uncertainties, in a way that makes feasible the adaptation of proofs from stochastic approximation by differential inclusions [START_REF] Roth | Stochastic approximations with constant step size and differential inclusions[END_REF].

II. IMPRECISE MARKOV CHAINS

In this section we first discuss imprecise continuous-time Markov chains in general. Population processes are a subclass of this general model. We will introduce the imprecise and the uncertain models. We will then introduce briefly differential inclusions, and show how Kolmogorov equations for the probability mass generalise to differential inclusions in this setting. We also introduce the imprecise drift of an imprecise model, which will play a central role in the construction of mean field limits.

A. Imprecise and uncertain Markov chains

We consider a stochastic process X = (X t ) t≥0 that takes value in a discrete state space E ⊂ R d and that is adapted to a filtration1 F. The dynamics of the process depends on a parameter (or a vector of parameters) ϑ. We denote by Θ the set of possible parameter values of ϑ. We consider a set of transitions kernel Q ϑ on E, parametrized by ϑ ∈ Θ: For each ϑ ∈ Θ, Q ϑ is such that Q ϑ x,y ≥ 0 for x = y ∈ E and y∈E Q ϑ xy = 0. Definition 1. An imprecise continuous time Markov chain is a stochastic process X together with a F t -adapted process θ such that for all t ≥ 0:

lim h→0 1 h P(X t+h = y | F t , X t = x) = Q ϑ t xy if x = y -y =x Q ϑ t xy if x = y
The definition of an imprecise Markov chain makes no restriction on the set of processes that the varying parameter θ can take. In some cases, it can be interesting to focus on subset of processes. An example is if we assume that ϑ t is deterministic and constant in time. In that case we obtain the notion of uncertain Markov chain: Definition 2. An uncertain continuous time Markov chain is a stochastic process X such that there exists a constant parameter ϑ ∈ Θ such that for all t ≥ 0:

lim h→0 1 h P(X t+h = y | F t , X t = x) = Q ϑ xy if x = y -y =x Q ϑ xy if x = y
The definitions of imprecise and uncertain CTMC correspond to two extreme cases: in the imprecise case, the parameter function θ can be any measurable function while in the uncertain case the parameter function θ is deterministic and constant in time. It would be possible to consider cases in between, for example by restricting the set of admissible processes θ to be the set of functions that only depends on the value X(t) -such process θ are called Markovian control policies in the Markov decision processes community -or the set of time-dependent deterministic functions -which would lead to time-inhomogeneous CTMC.

Remark: Imprecise CTMC are strongly related to (continuous-time) Markov Decision Processes (MDPs) [START_REF] Baier | Efficient computation of time-bounded reachability probabilities in uniform continuous-time Markov decision processes[END_REF], the main difference being that in MDPs, the emphasis is on finding a policy that maximises some reward criteria. To to this, one usually assumes a finite or countable number of actions that can be taken non-deterministically. In ICTMC, the variability of the decision variables is usually not controllable and can be a generic random function adapted to the process. The parameter space is in general uncountable and our objective is to characterise the set of possible behaviours. Uncertainty, in the sense considered here, has been considered also in the MDS context, in the so called Uncertain (Bounded-Parameter) MDPs [START_REF] Beccuti | Markov decision petri nets with uncertainty[END_REF], [START_REF] Givan | Bounded-parameter markov decision processes[END_REF].

Example: For illustrative purposes, we consider a simple model of a bike sharing system, similar to the one of [START_REF] Fricker | Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity[END_REF]. We describe the number of bikes present in a single station, so that the state of the CTMC is given by X t ∈ {0, . . . , N }, where N is the capacity of the station, i.e. the maximum number of bikes. We assume that customers arrive at an unknown rate ϑ a , belonging to the interval [ϑ min a , ϑ max a ]. Each customer arrival brings the system from state k to k -1 (for k > 0). Similarly, we can model the return of a bike as a transition with rate ϑ r , belonging to the interval [ϑ min r , ϑ max r ], and increasing the number of available bikes from k to k + 1, provided k < N . Different choices of how the two imprecise parameters can vary give rise to different models. In the general case, we can assume that ϑ a (t) and ϑ r (t) are generic functions of time, encoding complex but unknown dependencies of the environment of customers requesting a bike and wishing to return it at the station. On the opposite spectrum, we could assume that these values are unknown but constant in time, for example if we are considering the dynamics of the station in a restricted time frame.

B. Differential inclusions

Differential inclusions (DI) are a generalisation of differential equations, and provide the natural mathematical tool to describe the transient evolution of the probability mass of an imprecise CTMC. Furthermore, as we will see in the next section, under a proper scaling, the behaviour of an imprecise Markov chain is closely related to the one of a differential inclusion corresponding to Equation [START_REF] Mari | Modelling cholera epidemics: The role of waterways, human mobility and sanitation[END_REF]. In this section, we provide a quick introduction to DI, recalling some classical definitions. See [START_REF] Aubin | Differential Inclusions[END_REF] for further details.

Let F be a set-valued function on E ⊂ R d that assigns to each x ∈ E a set of vectors F (x) ⊂ R d . A solution to the differential inclusion ẋ ∈ F (x) that starts in x is a function x : [0, ∞) → E such that there exists a measurable function

f : [0, ∞) → R d satisfying: for all t ≥ 0 f (t) ∈ F (x(t)) and x(t) = x 0 + t 0 f (s)ds.
For an initial condition x 0 , we denote by S F,x0 (or S x0 if there is no ambiguity) the set of solutions of ẋ ∈ F (x) that start in x 0 . Note that the set S F,x0 can be empty or be composed of multiple solutions, depending on the function F . When E = R d , a sufficient condition for the existence of at least one solution is that (a) for all x ∈ R d F (x) is non-empty, convex and bounded (i.e., sup x∈X,y∈F (x) y < ∞) and (b) F is upper-semi-continuous (i.e., the graph of F , {(x, y), y ∈ F (x)}, is a closed set).

Asymptotic behaviour as t goes to infinity: As the time grows, a solution x of a differential inclusion can have a chaotic behaviour and can oscillate around many points. For a given starting point x 0 , we define the limit set L F (x 0 ) as the closure of the set of points that are accumulation points of at least one solution of the differential inclusion starting in x 0 :

L F (x 0 ) := x∈S F,x 0 t≥0 closure({x(s), s ≥ t}),
where closure denotes the closure in the metric space R d .

As in [START_REF] Roth | Stochastic approximations with constant step size and differential inclusions[END_REF], we call the Birkhoff centre of the differential inclusion is the closure of the set of recurrent points of the differential inclusion:

B F = closure({x ∈ E : x ∈ L F (x)}). (1)
Intuitively, the Birkhoff centre contains all attractors, equilibria, limit cycles, and in general all points that can be visited infinitely often by one solution of the differential inclusion. As such, it provides a characterisation of the localisation in the phase space of the limit behaviour of the system. Note that, in general, the sets L F (x) and B F are not convex nor necessarily connected, even when F (x) is single valued and Lipschitz (for example the Birkhoff centre of

F (x) = {x(1 -x)} is {0, 1}).
Note that even if the case of Lipschitz ordinary differential equation, the computation of the Birkhoff centre of an ODE is complicated and even its shape can depend strongly on the parameter of the function (see the comparison between Figure 1 and Figure 2 of [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF]).

C. Imprecise Kolmogorov Equations

The time evolution of the probability mass of an imprecise Markov chain can be obtained by a generalised version of Kolmogorov equations, stated in terms of differential inclusions. For a fixed realisation ϑ t of the process ϑ, the evolution of the probability mass P (t | ϑ t ) is given by the standard (non-autonomous) Kolmogorov (forward) equations

Ṗ (t | ϑ t ) = Q ϑt •P (t | ϑ t ).
However, if we have no knowledge on the process ϑ, we can only assume that the process at time t behaves like a Q ϑ for some ϑ, implying that the transient behaviour of the probability mass will be a solution of the differential inclusion

Ṗ (t) ∈ Q • P (t), (2)
where 2) is linear, but the dimension of the differential inclusion equals the size of the state space, which is too large in most interesting scenarios, especially for population models, see Section III.

Q = ϑ∈Θ Q ϑ . Equation (

D. Generalised Drift

An important notion associated with an imprecise Markov chain is that of drift, which extends the corresponding notion for CTMCs, and describes the average increment in a point of the state space, as a function of the imprecise parameters.

Definition 3. The imprecise drift of the imprecise CTMC is the function

f : E × Θ → R d , defined by f (x, ϑ) = y∈E Q ϑ x,y (y -x). (3)
When for all x ∈ E, ϑ ∈ Θ we have y∈E Q ϑ x,y yx < ∞, the sum in Equation (3) is well defined.

III. IMPRECISE POPULATION PROCESSES

In this section we introduce the class of models of interest in this paper, namely imprecise population processes. Population processes are a very common class of models, and when relaxing the precision intrinsic in their definition, Imprecise population processes emerge naturally. After providing the definition, we shift the focus to the behaviour of imprecise population processes in the limit of an infinite population. We show that such a limit is described by a differential inclusion, which provides also information about the stationary behaviour of the imprecise process, encoded in the Birkhoff centre.

A. Definition

Let N be a scaling parameter (typically, N is the population size of the considered model). We consider a sequence of imprecise Markov processes indexed by N , denoted (X N ) N ≥0 on a sequence of subset E N ⊂ E ⊂ R d . The stochastic process X N is an imprecise process of kernel Q N,θ and X N (t) denotes the state of the process at time t. Definition 4. An imprecise (respectively uncertain) population process is a sequence of imprecise (respectively uncertain) Markov chains that satisfies the following assumptions:

(i) The chains are uniformizable: i.e., for all N , the rates are uniformly bounded:

sup x∈E N ,ϑ∈Θ |Q N,ϑ xx | < ∞ (ii)
The transitions become smaller as N grows, i.e., there exists ε > 0 such that

lim N →∞ sup x∈E N ,ϑ∈Θ y∈E N Q N,ϑ xy y -x 1+ε = 0 (iii)
The drifts are well-defined and bounded:

lim sup N →∞ sup x∈E N ,ϑ∈Θ y∈E N Q N,ϑ xy y -x < ∞
Transition classes: A direct definition of an imprecise population process by exhibiting its generator is unfeasible. A simpler definition can be obtained by specifying transition classes, similarly to [START_REF] Bortolussi | Continuous approximation of collective systems behaviour: A tutorial[END_REF], and describing the state of the system by a state vector, counting the size of each population. The idea is to specify the dynamics by a set of possible transitions or events, providing their rate, as a function of the state vector and of the imprecise parameters, and how they change the state vector of the system. Specification can be done at the level of the density (normalised population process) or on the integer-valued counting variables (non-normalised process). Usually, such specifications satisfy the conditions in the previous definition. See below and Section V for examples. Example: A simple example of a population model is provided by a slight modification of the single station bike sharing system discussed in the previous section. In this case, we can assume the total population N is given by the total amount of bike racks in the station, while the only model variable, X B (t), encodes the number of bikes available at the station in terms of the fraction of occupied bike racks. There are two transitions classes in the model: a customer arrival, taking one bike, and a biker arrival, returning one bike. The former transition class has rate N ϑ a (t), if X B (t) > 0, and changes the state of the system from X B to X B -1/N . Similarly, bike arrivals happens at rate N ϑ r (t), if X B (t) < 1, and move the system from X B to X B +1/N . The dependency of the arrival rates on N is needed to satisfy the conditions of the definition above, and intuitively describes the idea that ϑ a and ϑ r model the relative traffic volume in terms of bike density, so that the absolute traffic level is obtained by multiplication with the system size N .

B. Mean field limit (finite time-horizon)

Let f N (x, ϑ) = y∈E N Q N,ϑ
xy (yx) be the drift of the system N . We define the limit drift of the system as the convex closure of the set of the accumulation points of f N (x N , ϑ) as N goes to infinity, for all sequences x N that converge to x:

F (x) = lim N →∞ ϑ∈Θ f N (x, ϑ) ,
= y ∈ R d such that there exists a sequence φ N : lim

N →∞ φ N = ∞ ∧ lim N →∞ f φ N (x φ N , ϑ φ N ) = y . (4)
As in [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF], it can be shown that when the Q N,ϑ 's satisfy Definition 4, the differential inclusion ẋ ∈ F (x) has a unique solution as for all x ∈ R d F (x) is non-empty, convex and bounded (i.e., sup x∈X,y∈F (x) y < ∞) and F is upper-semicontinuous.

Theorem 1. Let (X N ) be an imprecise population process. Then, if X N (0) converges (in probability) to a point x, then the stochastic process X N converges (in probability) to S F,x , the set of solutions of ẋ ∈ F (x) starting in x.

Sketch of proof:

The difficult part of the proof is to show that there exists a decreasing function I(N ) such that the discrete time process X N (kI(N )) satisfies the definition of a GASP of [START_REF] Roth | Stochastic approximations with constant step size and differential inclusions[END_REF]. This can be done similarly to [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF]. Hence, the process (X N (kI(N ))) k satisfies the assumptions of [18, Theorem 3.1] which concludes the proof.

Remark: The construction of the imprecise drift generalises the definition given in [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF], where the focus was on discontinuous, but precise, rates. The proof of the theorem follows the ideas of [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF], [START_REF] Benaïm | Stochastic approximations and differential inclusions[END_REF]. In fact, when Θ is restricted to a single value, we obtain again the result of [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF].

In the case of an uncertain population process, we define the drift F ϑ (x) as a function of the state x ∈ E and the parameter ϑ ∈ Θ.

F ϑ (x) = lim N →∞ f N (x, ϑ) . (5)
It should be clear that for all x: ϑ∈Θ F ϑ (x) ⊂ F (x). Note that the quantity F ϑ (x) is non-necessarily single-valued when the dynamics is not-continuous (see [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF]).

Corollary 1. Let X N be an uncertain Markov population process, then if X N (0) converges (in probability) to a point x, then the stochastic process X N converges (in probability) to ϑ S ϑ,x , where S ϑ,x is the set of solutions of the differential inclusion ẋ ∈ F ϑ (x) starting in x.

These results are illustrated in Section V.

C. Mean field limit (stationary regime)

An uncertain Markov chain does not necessarily have a stationary behaviour. In particular, an imprecise Markov chain is not necessarily time-homogeneous (the parameter ϑ t can depend on time) and is also not necessarily a Markov chain (ϑ t can depend on all the past of the stochastic process). Nevertheless, it is possible to constraint the asymptotic regime of an imprecise Markov population process by using the asymptotic reachable set of the differential inclusion, A F , defined by:

A F = T >0 x0,t≥T,x∈S F,x 0 {x(t)} (6)
Note that the the Birkhoff centre B F is included in the set A F . The inclusion is in general strict.

The following result is a direct consequence of Theorem 1.

Theorem 2. Let (X N ) N be an imprecise population process, then there exists a sequence ε N that converges to 0 in probability and such that the distance between X N (t) and A F becomes stochastic lower than ε N as t goes to infinity. In other word:

lim N →∞ lim t→∞ d(X N (t), A F ) = 0 in probability.
When in addition of being an imprecise Markov chain, for all N , the process X N is a Markov chain and has a stationary measure µ N , we can say more. The following results characterise the stationary behaviour of the stationary measures as N grows. It shows that the sequence of stationary measures µ N concentrates on the Birkhoff centre of the differential inclusion. The Birkhoff centre, defined in Section II-B, is essentially the set of recurrent points of the differential inclusion, i.e., the set of points such that there exists a trajectory that starts at this point and comes back to this point in the future. Theorem 3. Let X be an imprecise population process such that X N is a Markov chain that has a stationary measure µ N . Let µ be a limit point of µ N (for the weak convergence). Then, the support of µ is included in the Birkhoff centre of F , defined in Equation (1): µ(B F ) = 1.

Sketch of proof:

As for the proof of Theorem 1, one can show that there exists a function I(N ) such that X N (kI(N )) is a GASP. Theorem 3 is therefore a consequence of [START_REF] Roth | Stochastic approximations with constant step size and differential inclusions[END_REF]Corollary 9].

Corollary 2. Let X be an uncertain population process such that X N is a Markov chain that has a stationary measure µ N . Let µ be a limit point of µ N (for the weak convergence). Then, the support of µ is included in the Birkhoff centre of F ϑ , defined in Equation (1).

Remark: Theorem 3 only states that the support of µ is included in the Birkhoff centre but provides no intuition on the how the probability mass is spread on this set. This results can be refined by using the notion of semi-invariant measure [START_REF] Roth | Stochastic approximations with constant step size and differential inclusions[END_REF]Definition 3.3]. However, this notion is complex and computing a semi-invariant measure of a differential inclusion seems numerically intractable, as it requires to compute a probability measure on all the possible trajectories of the differential inclusion. Hence, it the present document, we limit our exposition to the notion of Birkhoff centre, which, even if less accurate, provides a simpler characterisation that can be computed numerically.

IV. NUMERICAL METHODS

The results of the previous section imply that, for large populations, we can study the mean field differential inclusion to get insights on the transient and on the steady state of the stochastic population process. The analysis of this class of mean field models, however, is in general more challenging than ordinary differential equations (ODE). In the differential inclusion case, in fact, we usually are only able to compute bounds on the solution set of the equations. After discussing some existing approaches, we present in more detail two fast methods, one based on differential hulls and the other based on the control-theoretic Pontraygin principle.

A. Related work on numerical methods

Most recent numerical approaches dealing with imprecise deterministic processes have been developed for computing reachable sets of hybrid systems, whose continuous dynamics can be specified by (non-linear) differential equations or differential inclusions. The proposed methods in literature can be roughly divided in two classes: exact over-approximation methods and simulation-based methods [START_REF] Maler | Computing reachable sets: An introduction[END_REF].

The first class of methods manipulates directly sets of states, finitely represented, for instance, as polytopes [START_REF] Schutter | Survey of modeling, analysis, and control of hybrid systems[END_REF], ellipsoids [START_REF] Kurzhanski | On ellipsoidal techniques for reachability analysis[END_REF], or zonotopes [START_REF] Girard | Efficient computation of reachable sets of linear time-invariant systems with inputs[END_REF], or by relying on interval arithmetic [START_REF] Alefeld | Interval analysis: Theory and applications[END_REF] combined with constraint solving [START_REF] Chen | Flow*: An analyzer for non-linear hybrid systems[END_REF], [START_REF] Ramdani | Computing reachable sets for uncertain nonlinear hybrid systems using interval constraint-propagation techniques[END_REF], exploiting satisfaction modulo solvers over reals [START_REF] Gao | Dreal: An SMT solver for nonlinear theories over the reals[END_REF], [START_REF] Moura | Z3: An Efficient SMT Solver[END_REF]. The dynamics of the system is lifted at the set level, so that one computes the evolution of the reachable set under the action of the dynamics. These methods usually compute overapproximation of the real reachable set, which is formally guaranteed to contain all reachable points. Most of the methods in literature are restricted to linear systems, due to their wide diffusion in engineering applications. For this class of systems, efficient methods exist [START_REF] Girard | Efficient computation of reachable sets of linear time-invariant systems with inputs[END_REF], [START_REF] Botchkarev | Verification of hybrid systems with linear differential inclusions using ellipsoidal approximations[END_REF]. Tools like SpaceEx [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF] implement routines for linear differential inclusions capable of solving problems up to several hundred dimensions. Overapproximation methods for non-linear systems, instead, are much less developed, due to the intrinsic difficulty of the problem. Here we recall hybridization [START_REF] Dang | Computing reachable states for nonlinear biological models[END_REF], which is based on a localised linearisation of the dynamics. Another family of methods, very common for differential inclusions, is based on interval-based methods combined with constrain solving to control the over-approximation introduced by interval computations. Among them, we recall [START_REF] Chen | Flow*: An analyzer for non-linear hybrid systems[END_REF] and [START_REF] Ramdani | Computing reachable sets for uncertain nonlinear hybrid systems using interval constraint-propagation techniques[END_REF]. Other methods for integrating differential inclusions rely on error bounding [START_REF] Zgliczynski | Lohner algorithm for perturbation of odes and differential inclusions[END_REF].

An alternative approach is offered by simulation methods, which try to infer the reachable set from few simulated trajectories of the system. These methods are very effective for the limit dynamics emerging from uncertain population processes. Here we recall [START_REF] Donzé | Systematic simulation using sensitivity analysis[END_REF], which uses the sensitivity of the system with respect to initial conditions to (approximately) compute how a set of initial conditions propagates in T units of time under the action of the dynamics. The approach of [START_REF] Bortolussi | A Statistical Approach for Computing Reachability of Non-linear and Stochastic Dynamical Systems[END_REF], instead, use statistical methods borrowed from machine learning to perform inference of the reachable set with statistical error guarantees. Other methods work for more general imprecise limit models; for example, the procedure of [START_REF] Bhatia | Incremental search methods for reachability analysis of continuous and hybrid systems[END_REF]- [START_REF] Dang | Coverage-guided test generation for continuous and hybrid systems[END_REF] constructs an under-approximation of the reachable set using a Monte-Carlo sampling method.

B. Differential hull of a differential inclusion

A first method to compute the set of reachable points by the solutions of a differential inclusion is to construct a rectangular over-approximation, which is composed of two functions x and x such that x(t) ≤ x(t) ≤ x(t), for any x that is solution of the differential inclusion x. In this section, we show that we can construct two differential equations such that x and x are their solutions and such that for all t and any solution x of the differential inclusion ẋ ∈ F (x), we have x(t) ≤ x(t) ≤ x(t) (coordinate-wise). The construction of these equations is simple. This bounds are reasonably tight when the set of possible parameters θ is small. However, as we will see in Section V-D, the bounds provided by this approximation become loose when the set of possible ϑ increase.

Following the definition of [START_REF] Tschaikowski | Approximate reduction of heterogenous nonlinear models with differential hulls[END_REF], we say that a locally Lipschitz-continuous function (f , f ) is a differential hull for the differential inclusion F if for each coordinate i and each x such that x ≤ x ≤ x, we have:

x i = x i ⇒f i (x, x) ≤ inf f ∈F (x) f i x i = x i ⇒f i (x, x) ≥ sup f ∈F (x) f i .
The tightest functions f and f that are a differential hull for a differential inclusion F are:

f i (x, x) = min x∈[x,x]:xi=x i (t) min F i (x) f i (x, x) = max x∈[x,x]:xi=x i (t) max F i (x)
This leads to the following result, that is a consequence of [42, Theorem 1]. Theorem 4. Let x : [0 : T ] → R d be a solution of ẋ ∈ F (x) with initial condition x(0) = x 0 . Let (x, x) be the solution of the differential equation ẋ = f (x, x) and ẋ = f (x, x) with initial condition x(0) = x(0) = x 0 . Then, for all t, we have:

x(t) ≤ x(t) ≤ x(t)
C. Reachability as an optimal control problem: an algorithm based on Pontryagin's maximum principle In general, the bounds x and x found by the differential hull method are not tight. We can do better by exploiting Pontryagin's maximum principle, a classical method that can be used to compute numerically the exact minimal value x min (t) and maximal value x max (t) that can be reached by an imprecise fluid model at time t [START_REF] Todorov | Optimal control theory[END_REF].

Let T ≥ 0 be some fixed time and i ∈ {1 . . . d} a coordinate. Let x max i (T ) = sup x∈S F,x 0 x i (T ) be the maximal value that the ith coordinate of the solution of a differential inclusion can take at time t. The quantity x max i (T ) is the solution of the maximization problem over functions ϑ(t):

x max i (T ) := max θ x i (T ) such that for all t ∈ [0; T ]: x(t) = x + t 0 f (x(s), θ(s))ds θ(t) ∈ [ϑ min , ϑ max ]
Pontryagin's maximum principle is a set of necessary conditions that the trajectory that attains the maximum should satisfy. Following the description of [7, Section 3] these conditions are the following. If x is a trajectory that maximises x i (T ), then there exists a costate trajectory p such that p i (T ) = -1, p j (T ) = 0 for j = i and:

ẋ(t) = f (x(t), u(t)) (7) ϑ(t) ∈ arg max ϑ f (x, ϑ) T p, (8)
-

ṗ(t) = ∂ ∂x f (x, ϑ) T p (9)
where f (x, ϑ) T p denotes the scalar product between f (x, ϑ) and p.

This formulation can be used to derive a fixpoint algorithm to iteratively obtain x and ϑ. We start from an (arbitrary) initial costate p, satisfying the constraints at time T , then alternatively compute a trajectory x forward in time the x by using ( 7) and (8), holding p fixed, and computing a new p backward in time by solving [START_REF] Kisielewicz | Stochastic differential inclusions and applications[END_REF], holding x fixed. These two steps are repeted until a convergence criterion is met (typically, the sup norm of p and x varies less than a prescribed ε after an iteration).

Remark: The rectangle delimited by x min (T ) and x max (T ) provides an approximation of the reachable set of x(T ). This set is tighter than the rectangle delimited by x and x, but is not exact: in general, the set of reachable points is not a rectangle. However, the algorithm given by the iterations (7-9) can be easily extended to refine the rectangle into any convex template polyhedron by considering the minimisation problems min i i α i x i (t) for any tuple of coefficients α i .

V. ILLUSTRATIVE EXAMPLE 1: THE SIR MODEL

In this section, we apply our techniques to the well-known susceptible-infected-recovered (SIR) model. This example will serve us to illustrate the differences between the imprecise and uncertain model and to evaluate our numerical algorithms. We show in particular that the imprecise model can be solved numerically and that it provides bounds that are much tighter than the differential hull approximation.

A. The SIR model

We consider a system composed of N nodes. Each node is in one of the three states: susceptible, infected or recovered. We denote by X S (t) ∈ [0, 1] the proportion of nodes that are susceptible at time t (and by X I (t) and X R (t) the proportion of infected and recovered nodes). By definition, for each t, X S (t) + X I (t) + X R (t) = 1.

We model the dynamics of the stochastic system as follows. A susceptible node can become infected from an external source (this occurs at rate a). An infected node becomes recovered at rate b and a recovered node becomes susceptible at rate c. We consider that susceptible and infected node are moving but we do not know at which speed. We model this by considering that a susceptible node encounters an infected node and becomes infected at rate ϑX I (t).

The transitions of the system are the following: The state X S , X I , X R becomes:

• X S -1 N , I + 1 N , X R at rate N (aX S + ϑX S X I ), • X S , X I -1 N , X R + 1 N at rate N bX I and • X S + 1 N , X I , X R -1 N at rate N cX R . The sequence of processes (X N S , X N I , X N R
) N satisfies Definition 4 and is therefore an imprecise population process. Its drift is given by the triple

(f S (X S , X I , X R ; ϑ), f I (X S , X I , X R ; ϑ), f R (X S , X I , X R ; ϑ)), where: f S (X S , X I , X R ; ϑ) = -aX S -ϑX S X I + cX R f I (X S , X I , X R ; ϑ) = aX S + ϑX S X I -bX I f R (X S , X I , X R ; ϑ) = bX I -cX R (10)
As X S + X I + X R = 1, we can substitute for X R and express the drift as

f S (X S , X I ; ϑ) = c -(a + c)X S -cX I -ϑX S X I f I (X S , X I ; ϑ) = aX S + ϑX S X I -bX I (11)
By Theorem 1 and Corollary 1, as the number of object goes to infinity, the behaviour of the uncertain model (an unknown but constant ϑ) converges to the ODE ẋ = f (x, ϑ) while the imprecise model (a unknown θ that can vary in time) is included in the solution of the differential inclusion

ẋ ∈ ϑ∈[ϑmin,ϑmax] {f (x, ϑ)}
In what follows, we will compare numerically the behaviour of the two models. In particular, we show that despite the fact that the infection rate minus the recovery rate f I (X S , X I ; ϑ) = aX S + ϑX S X I -bX I is an increasing function of ϑ, the quantity x I (t) is not a monotone function of ϑ. In particular, when ϑ(t) ∈ [ϑ min, ϑ max], the proportion of infected nodes can be higher for an imprecise population process than for any uncertain population process.

In the rest of this section, we set the parameters equal to a = 0.1, b = 5, c = 1, ϑ min = 1, ϑ max = 10 and the initial conditions are X S (0) = 0.7, X I (0) = 0.3, X R (0) = 0.

B. Reachable sets in finite time

For a fixed parameter ϑ, the drift f (X S , X I , ϑ) is Lipschitz-continuous in X S and X I . Hence, the drift of the uncertain model, defined in Equation ( 5), is a single-valued function and the ODE ( Ṡ, İ) = f (X S , X I , ϑ) has a unique solution, which we denote by X ϑ S , X ϑ I . This ODE does not have a close-form solution but numerical integration is easy.

The case of the imprecise model is more complicated. By Theorem 1, whatever are the variation of the parameter θ (even if it ϑ t depends on the whole history of the process X N S , X N I ), the quantities X N S and X N I converges to a solution of the differential inclusion ( ẋS , ẋI ) ∈ {f (x S , x I , ϑ) : ϑ ∈ [ϑ min , ϑ max ]}. In particular, if S F is the set of solutions of the differential inclusions and R F (t) is the set of reachable points at some time t by the differential inclusion starting in (X S (0), X I (0)), we have:

lim N →∞ (X N S (t), X N I (t)) ∈ R F (t) := (s,i)∈S F {(s(t), i(t))}
It should be clear that the set of the possible values for (X S , X I ) of an uncertain model is included in the set of reachable points of the imprecise model:

ϑ∈[ϑmin,ϑmax] (x ϑ S (t), x ϑ I (t)) ⊂ R F (t). (12)
As we show numerically in our example (see Figure 1), the inclusion is, in general, strict.

For both models (uncertain and imprecise), we define the maximum proportion of infected nodes by

x uncertain I (t) = max ϑ I ϑ (t) x imprecise I (t) = max (X S ,X I )∈S F X I (t).
The definition of the minimum proportion of infected node is similar.

The computation of x uncertain I (t) can be done by a numerical exploration of all the parameters ϑ, or relying on the statistical method of [START_REF] Bortolussi | A Statistical Approach for Computing Reachability of Non-linear and Stochastic Dynamical Systems[END_REF]. For the computation x imprecise I (t), we formulate the problem as an optimal control problem which is to find a function θ(t) that maximises X I (t), where X I (t) satisfies ẋS , ẋI = f (x S (t), x I (t), θ(t)). We use a numerical method based on Pontryagin's maximum principle (see Section IV-C). The results are reported in Figure 1. We observe that x imprecise I (t) can be much larger than x uncertain I (t), especially for large value of t.

On Figure 2, we add to Figure 1 two examples of trajectories that attain the maximum (top) or minimum (bottom) number of infected nodes at time T = 3. These maximal and minimal trajectories are in fact solutions of bang-bang control policies. For the trajectory that maximises X I (3), the control θ(t) equals ϑ min for t < 2.25 and ϑ max for t > 2.25. For the trajectory that minimises X I (3), the control θ(t) is to start with ϑ min for t < 0.7, then use ϑ max until t = 2.2 and then use again ϑ min .

C. Steady-state regime

Computing the Birkhoff centre B F of a differential inclusion is in general not easy. A first possibility is to use Pontryagin's principle to compute the convex hull of the set of reachable points at time t. By letting t goes to infinity, this allows one to compute the convex hull of the asymptotic set A F (defined in Equation ( 6)), which is a super-set of the Birkhoff centre B F .

For the SIR model, we use a more direct approach, that is faster and more accurate for a two dimensional system. By integrating the ODE ẋ = F (x, ϑ max ), we first compute a point x 0 that is the fixed point of the uncertain model with fixed parameter ϑ max . We then compute a trajectory x 1 by integrating the ODE ẋ = F (x, ϑ min ) starting in x 0 and a trajectory x 2 by integrating the ODE ẋ = F (x, ϑ max ) starting in x 1 (∞). The two curves delimit a convex region that is included in the Birkhoff centre. We then start from any points on the surface of this region and we look for a value ϑ such that the drift is directed outside this region. If such a value exists, we then enlarge our region by computing a new trajectory with parameter ϑ. We repeat these iterations until no such point on the surface exists. When this region cannot be enlarged, we then obtain the Birkhoff centre of the differential inclusion. This is so because the drift vector never points outwards the region in any point of its boundary, meaning that no trajectory can escape from it.

In Figure 3, we compare the Birkhoff centre of the imprecise fluid model with the one of the uncertain fluid model. The steady-state of the uncertain fluid model is shown in red while the steady-state of the imprecise model corresponds to the whole convex region surrounded by the blue curve. We observe that the steady-state of the uncertain model is strictly included in the one of the imprecise model. Moreover, there are some points of the Birkhoff centre of the imprecise model for which X S is smaller and X I is larger that any stationary point of the uncertain model.

D. Comparison with the differential-hull approximation

It should be noted that the algorithms based on the Pontryagin's maximum principle provide an exact numerical method to compute the transient and steady-state behaviour of these systems. On the contrary, the differential-hull approximation, introduced in Section IV-B provides a reasonably accurate results when the range of the parameter θ is small but its accuracy is very poor when the range of parameters increases.

These facts are illustrated in Figure 5 and Figure 4. In Figure 4, we plot the solutions as a function of time. We observe that when the range of possible ϑ is small (ϑ max = 2), then the differential hull approximation is quite accurate. When ϑ max grows, the approximation is less and less accurate. For example, for ϑ max = 5, the approximation is that X I (t) ∈ [.02, 1.17]. When ϑ max = 6, the approximation is trivial for t ≥ 4, for which we have X I (t) = 0 ≤ X I (t) ≤ X I (t) = 1.

In Figure 5, we compare the possible steady-states for the imprecise model and the uncertain model with the bounds obtained by the differential-hull approach. The differential-hull approach provides a rectangular approximation for the steadystate distribution. As for the transient regime, the differentialhull approximation deteriorates non linearly in ϑ max : it is very accurate for ϑ max = 2 or 3 but very loose for ϑ max = 5 (and trivial for ϑ max ≥ 6 -not shown on the figure).

E. Comparison between the mean field approximation and stochastic simulations

Theorem 3 guarantees that when N is large, for any control policy θ, the stationary measure of the stochastic system concentrates on a subset of the Birkhoff centre as N goes to infinity. In this section, we illustrate this result by comparing a stochastic simulation of an imprecise Markov population process of the SIR model for a finite N with the limiting regime. We implemented a simulator of the SIR model with two different control functions θ(t) that are piece-wise constant:

(a) A function θ 1 that oscillates between ϑ min and ϑ max according to the following rule: if θ 2 (t) = ϑ max and X S (t) < . These functions are chosen as they induce large variations of ϑ, forcing large oscillations of the stochastic system. The control function θ 1 induces an almost-periodic behaviour of the mean field approximation. The function θ 2 is more random.

A sample path in the steady state regime of the stochastic system for N = 100, N = 1000 and N = 10000 is depicted on Figure 6, along with the Birkhoff centre. We observe that for N ≥ 1000, the stationary behaviour of the stochastic system essentially remains inside of the Birkhoff centre for both policies. The inclusion tends to be exact as N goes to infinity.

VI. ILLUSTRATIVE EXAMPLE: A GENERALISED PROCESSOR SHARING MODEL

In this section we consider another example from queueing theory, namely a closed tandem queueing network with two classes of jobs consisting of a queue which serves with GPS discipline and a delay (i.e., infinite server) station [START_REF] Iacobelli | Lumpability of fluid models with heterogeneous agent types[END_REF], [START_REF] Anselmi | Energy-aware capacity scaling in virtualized environments with performance guarantees[END_REF].

We consider an imprecise scenario for which the rate of jobs creation belongs to an interval [λ min , λ max ]. This example serves two purposes. The first is to demonstrate that our framework can be used to perform a robust optimisation of the system parameters. The second is to show in such systems, the precise distribution of the arrival process of jobs is fundamental: When jobs are created according to a Poisson process, the imprecise and uncertain case are similar. This is not true for more general arrival processes for which having a varying arrival rate can lead to a higher queue length than setting a constant λ = λ max .

A. The model

We consider a system in which N applications are treated by a single machines. There are two types of applications: There are N 1 applications of type 1 and N 2 = N -N 1 applications of type 2. Each application sends periodically a job to the machine (the jobs creation process will be described later). The jobs created by the applications of type i are called the jobs of type i. We assume that the size of a job is exponentially distributed with mean 1/µ i and that the capacity of the machine is C (i.e., on average the machines takes 1/(µ i C) unit time to serve a job of type i).

Let Q i be the number of jobs of type i ∈ {1, 2} that are waiting to be completed by the machine. The machine serves the jobs by using a GPS policy. This means it uses a fraction φ i /(φ 1 Q 1 + φ 2 Q 2 ) of its time to each job of type i, where φ i > 0 is a weight. The weight is chosen so as to give priority to a type of jobs. We refer to [START_REF] Anselmi | Energy-aware capacity scaling in virtualized environments with performance guarantees[END_REF] for a justification and an indepth analysis of the pertinence of this model in the context of cloud systems, where each application represents one virtual machine or one process running on the same physical machine.

An application that just sent a job waits until the job is completed. Once the job is received, the application becomes idle. It waits a time that we assume is randomly distributed and then sends a new job. In the numerical evaluation, we compare two possible distributions for the time between jobs:

• (Poisson) an application of type i waits an exponentially distributed time of mean 1/λ i .

• (MAP) an application of type i waits an exponentially distributed time of mean 1/a i before being active. An active application then sends a job after a time that is exponentially distributed with mean 1/λ i

We refer to the first scenario as "Poisson arrivals" and to the second one as "MAP" (Markov arrival process).

B. Population model

This models fits in our framework. For the Poisson arrival scenario, we describe this system by four variables, Q i , i = 1, 2, the GPS queue length, and D i , i = 1, 2, the number of applications that are waiting before sending a new job. We normalise the variables D i and Q i to be in [0, 1], such that D i + Q i = 1. Transitions of the system are the creation of a job of type i, bringing the system to state D i -1/N i , Q i +1/N i , at rate λ i D i , and the service of an job of type i, changing state to D i + 1/N i , Q i -1/N i . This happens at rate µ i φiNiQi φ1N1Q1+φ2N2Q2 C. The MAP scenario is modelled similarly by adding two other variables E 1 and E 2 to count the applications for which the job was completed but that are not yet active.

C. Numerical analysis

We consider now a scenario in which the arrival rate is not fixed, but rather an imprecise parameter, varying measurably in ) and λ i max = 1/(1/a i + a/λ max i ) which are such that the average time that an application waits in the Poisson and MAP scenario are equal. In each case, the initial state is Q 1 (0) = Q 2 (0) = 0.1.

We make two observations. First, for the Poisson scenario, the uncertain and the imprecise model give the same upper and lower bound on the queue length. This confirms the intuition that the higher is the arrival rate λ, the more congested is the system. A more surprising fact is that in the MAP case, the maximal queue length is significantly larger than in the imprecise scenario. This fact is counter-intuitive, yet it is a consequence of the delay introduced by the time that the application take to activate. As for the SIR example, these upper and lower bounds were computed by using Pontryagin's maximum principle.

Our imprecise approximation can also be used to robustly tune controllable parameters of our system. We consider a scenario in which we want to set the resource allocation weights φ 1 and φ 2 in order to minimise the maximum total queue length: Q(t) = max θ (Q θ 1 (t) + Q θ 2 (t)). Using the Pontryagin's principle, we can compute Q(t) and minimise it numerically with respect to the weights. As only the ratio of φ 1 /φ 2 plays a role, we fixed the parameter φ 2 = 1. In turns out that the total queue length is a convex function of φ 1 and that the minimum value of Q(t) is obtained for φ 1 = 9.0φ 2 .

VII. CONCLUSIONS

We presented imprecise population processes, which naturally capture the unavoidable uncertainty and imprecision inherent in any (stochastic) model of complex phenomena. In order to analyse efficiently these models in the large population regime, we proved mean field theorems in terms of differential inclusions. We then discussed how to numerically analyse such limit systems, and show the approximation at work in a model of epidemic spreading and of generalised processor sharing.

Future work include a deeper investigation of the Pontraygin method also for steady state computation and in combination with (non-linear) templates, to provide tight bounds to the solutions of a differential inclusion. We will also release an implementation and test the approach on larger models, to properly understand its scalability.

Fig. 1 :

 1 Fig.1: Upper and lower bounds on the number of infected nodes for an imprecise model (dashed lines) and an uncertain model (solid lines)

Fig. 2 :

 2 Fig.2: Trajectory that reaches the maximum (top) or minimum (bottom) number of infected node at time T = 3. On the left, we plot the number of infected as a function of time. On the right, we plot the number of infected as a function of the number of susceptible.

Fig. 3 :

 3 Fig.3: Steady-state regime for the imprecise and uncertain SIR models. The steady-state of the imprecise model is the convex set delimited by the blue region. The steady-state of the uncertain model is on the red line. We set ϑ max = 10ϑ min .

Fig. 4 :

 4 Fig. 4: Evolution of the proportion of susceptible (left) or infected (right) as a function of time for three values of ϑ max . We compare the differential hull approximation (in dashed line) with the imprecise model (we set ϑ min = 1).

5 Fig. 5 :

 55 Fig. 5: Comparison of the possible steady-state of the imprecise model, the uncertain and the differential hulls approximation (we set ϑ min = 1).

(a) policy θ 1 (b) policy θ 2 Fig. 6 :

 126 Fig.6: Comparison of simulations of the stochastic system where θ varies with one of the two policies (a) and (b) of Section V-E with the Birkhoff centre of the differential inclusion (in blue). As N grows, the simulation gets included in the Birkhoff centre.

1 = 1 , λ max 1 = 7 , λ min 2 = 2 , λ max 2 = 3 ,

 11172223 We computed the maximal values of Q 1 (t) and Q 2 (t) as a function of time for the imprecise and the uncertain model. The results are reported in Figure7. The parameters of the simulations are µ 1 = 5, µ 2 = 1, φ 1 = φ 2 = 1, λ min a 1 = 1 and a 2 = 2. We set λ i min = 1/(1/a i + a/λ min i

Fig. 7 :

 7 Fig. 7: GPS model: Maximal queue length as a function of time for the uncertain and the imprecise model. This Figure is to be compared with Figure 1.

A filtration is a set of σ-algebra (Ft) t≥0 such that Fs ⊂ Ft ⊂ F for each t ≥ s ≥ 0. X is adapted to F means that Xt is Ft measurable for every t.

ACKNOWLEDGMENT

This work is supported by the EU project QUANTICOL, 600708.


The Interpretative Frame: the link between the Institution and its Artifact Carrier Artifacts can convey institutions [START_REF] Blanc | Against the digital revolution? Institutional maintenance and artefacts within the French recorded music industry[END_REF]. Or stated differently, actors can transform and manipulate material objects so they reflect and shape "cultural-cognitive, normative and regulative elements that provide stability and meaning to social life" (Scott, 2008: 222). Because it operates by representing the culture, values and symbols associated with a particular institution and/or by being infused with new institutional content, an artifact is named instantiation [START_REF] Hilpinen | Artifact[END_REF].

Through a case study implemented in Denmark and France, the paper studies how actors implement intervention works, such as contemporary adjustments -renovation or extension -, in listed buildings without disrupting their embodied Heritage. A listed building is a protected monument that highlights national pride or memory. In the study, six listed buildings are analysed: three in Denmark and three in France. The listed building's legitimacy relies on its authenticity whose respect by actors is essential to maintain the institutional protection, as it is the material representation, or instantiation of the Listed-Buildings Institution. However, intervention works to change such an artifact lead to various debates among actors, as the majority of current listed buildings were not originally constructed to last, i.e. to be transmitted to future generations [START_REF] Choay | L'allégorie du Patrimoine (4 th Edition)[END_REF]. One debate during intervention works tackles the issue of what needs to be or not to be considered in terms of Heritage.

Through a constructivist grounded-theory methodology [START_REF] Charmaz | Constructing Grounded Theory (2 nd Edition)[END_REF] and a Scandinavian Institutionalist lens (Czarniawska & Joerges, 1996), the research allows the understanding of how actors succeed in symbolically constructing the instantiation, so they can latter materially work on the artifact in a way it keeps conveying the institution [START_REF] Monteiro | Recovering materiality in institutional work: prizes as an assemblage of human of material entities[END_REF]. Indeed, within neo-institutional theory, current studies only focuses on how the modification of the instantiation can impact or change the institution [START_REF] Jones | Rebels with a cause: formation, contestation and expansion of the de Novo category "Modern Architecture" (1870-1975[END_REF]Raviola & Norbäck, 2013). What is still unknown is nevertheless how an instantiation can be changed while keeping conveying the institution that circumscribed it. In the study, I therefore emphasise that actors need first to design, with the help of the three institutional pillars given by Scott (2013), a shared interpretative frame to select relevant building materials. By doing so, they could thus modify an existing building in regards to what building's authenticity deserves to be respected. This frame acts in fact as a preliminary step to implement the materialisation of the intangible ideas.

The paper is thus focused on the interpretative frame that symbolically links the institution and its instantiation and on how actors design it. More specifically, I explore what constitutes such a frame and what its role is regarding the interrelations that exist between an institution, its artifact carrier and the actors who work on it. Indeed, if the Scandinavian Institutionalism literature already explains how actors can translate into practice an intangible idea, the question of the components of such an institutional frame remains overlooked [START_REF] Cornelissen | Putting framing in perspective: a review of framing and frame analysis across the management and organizational literature[END_REF]. That is why the dissertation contributes to the neo-institutional literature by arguing that such frame is built by means of the three institutional pillars, which are the components actors play with to know to what extent they can unfold action.

Throughout the study, I analyse the practice implemented by a collective of actors who have to tangibly modify an artifact, here a listed building, while keeping its instantiational character coming from the Listed-Buildings Institution. Consequently, to enable such one and only shared material practice, the paper underlines the importance of such an interpretative frame so the actors can share and intertwine their interpretations of the building's authenticity, i.e. the main leitmotiv on which institution of Listed-Buildings relies and takes its legitimacy from, in order to work towards the same goal. The aim for actors is thus to use the interpretative frame as a way to collectively interpret one specific but essential institutional feature in order to collectively do a practice that fits with it.

De facto, I argue that this a posteriori construction of the interpretative frame facilitates collective decision-making, as it acts as a shared and stabilised knowledge resource among actors. And by extension, I demonstrate and picture how the translation of the interpretative

Extended

  Abstract for OAP 2016 "Materiality & Institutions" Sylvain Colombero Post-Doc at Mines ParisTech -PSL Research University 60 boulevard Saint-Michel, 75272 Paris Cedex 06 sylvain.colombero@mines-paristech.fr

frame into an artifact reinforces the legitimacy of the institution and its taken-forgrantedness.