Nicolas Gast
email: nicolas.gast@inria.fr

Benny Van Houdt
email: benny.vanhoudt@uantwerpen.be.

A Refined Mean Field Approximation

Mean eld models are a popular means to approximate large and complex stochastic models that can be represented as N interacting objects. Recently it was shown that under very general conditions the steady-state expectation of any performance functional converges at rate O(1/N ) to its mean eld approximation. In this paper we establish a result that expresses the constant associated with this 1/N term. This constant can be computed easily as it is expressed in terms of the Jacobian and Hessian of the drift in the xed point and the solution of a single Lyapunov equation. This allows us to propose a rened mean eld approximation. By considering a variety of applications, that include coupon collector, load balancing and bin packing problems, we illustrate that the proposed rened mean eld approximation is signicantly more accurate that the classic mean eld approximation for small and moderate values of N : relative errors are often below 1% for systems with N = 10.

INTRODUCTION

Stochastic models have been used to assess the performance of computer (and many other) systems for many decades. As a direct analysis of large and complex stochastic models is often prohibitive, approximations methods to study their behavior have been devised. One very popular approximation method relies on mean eld theory. Its widespread use can be explained by the relative ease involved to dene and solve a mean eld model in combination with its high accuracy for large systems. More specically, many mean eld models are specied by a simple set of ordinary dierential or dierence equations and in many cases convergence of the stochastic system towards the solution of the mean eld model (over both nite and innite time scales) has been established as the number of components N in the system tends to innity, see [START_REF] Benaim | A class of mean eld interaction models for computer and communication systems[END_REF][START_REF] Bortolussi | Continuous approximation of collective system behaviour: A tutorial[END_REF][START_REF] Gast | Markov chains with discontinuous drifts have dierential inclusion limits[END_REF][START_REF] Thomas | Solutions of Ordinary Dierential Equations as Limits of Pure Jump Markov Processes[END_REF][START_REF] Thomas | Strong approximation theorems for density dependent Markov chains[END_REF][START_REF] Boudec | A generic mean eld convergence result for systems of interacting objects[END_REF]. Moreover these mean eld approximations may give rise to explicit expressions for the performance measures of interest (e.g., [START_REF] Minnebo | A Fair Comparison of Pull and Push Strategies in Large Distributed Networks[END_REF][START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF][START_REF] John | On the power of (even a little) centralization in distributed processing[END_REF][START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF]) or can be solved numerically without much eort (e.g., [START_REF] Gast | A Mean Field Model of Work Stealing in Large-scale Systems[END_REF][START_REF] Massoulié | Coupon Replication Systems[END_REF][START_REF] Van Houdt | A mean eld model for a class of garbage collection algorithms in ash-based solid state drives[END_REF]).

Although mean eld theory provides no guarantees with respect to the accuracy of a nite system of size N , there are many examples where the accuracy was demonstrated (using simulation) to be quite high even for systems of moderate size, e.g., for N ⇡ 50 or less. Nevertheless this accuracy very much depends on the exact parameter settings, for instance mean eld models for load balancing systems are known to be quite inaccurate under high loads even for moderate sized systems.

The main contribution of this paper is to establish a result that applies to a broad class of mean eld models, including density dependent population processes [START_REF] Thomas | Solutions of Ordinary Dierential Equations as Limits of Pure Jump Markov Processes[END_REF] and the class of discrete-time models of [START_REF] Benaim | A class of mean eld interaction models for computer and communication systems[END_REF], and that provides a signicantly more accurate approximation for nite N than the Coupon Supermarket Pull/push Simulation (N = 10)

1.530 2.804 2.304 Rened mean eld (N = 10) 1.517 2.751 2.295 Mean eld (N = 1)

1.250 2.353 1.636

Table 1. Illustration of the refined approximation accuracy. The columns correspond to the mean number of users with coupon i (K = 5, Table 2), the average queue length for the supermarket model ( = 0.9, Table 4) and for the pull/push model ( = 0.9, Table 7).

classical mean eld approximation. When considering a mean eld model X (N ) described by the density-dependent population process of Kurtz, where X (N ) i (t ) is the fraction of objects in state i at time t, the classical mean eld approximation shows that if the corresponding ODE model has a unique attractor , then

lim N !1 E (N ) f X (N ) g = 0.
This shows in particular that for a continuous function h, we have:

lim N !1 E (N ) f h(X (N ) ) g = h( ).
We establish conditions that show that for a function h, there exists a constant V h such that

E (N ) f h(X (N ) ) g = h( ) + V h N + o ✓ 1 N ◆ .
The constant V h is expressed as a function of the Jacobian and the Hessian matrices of the drift (in ) and the solution of a single Lyapunov equation. As such it can be computed eciently even when the number of possible states i of an object is large (less than one second for a 500dimensional model). We also establish that under an additional smoothness condition, the o(1/N ) can be expressed as V (2) h /N 2 + • • • + V (k ) h /N k + o(1/N k ), but the computation of these constants is computationally intensive.

The rened mean eld model that we propose consists in approximating E (N ) [h(X (N ) )] by h( ) + V h /N and maintains many of the attractive features of the classic mean eld approximation that resulted in its widespread use, but signicantly improves upon its accuracy for nite N . We demonstrate that this rened approximation signicantly improves the accuracy of the mean eld approximation by using a variety of mean eld models that include the coupon replication model [START_REF] Lin | Stochastic analysis of le-swarming systems[END_REF][START_REF] Massoulié | Coupon Replication Systems[END_REF] and a number of load balancing schemes [START_REF] Minnebo | A Fair Comparison of Pull and Push Strategies in Large Distributed Networks[END_REF][START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF][START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF][START_REF] Xie | Power of D Choices for Large-Scale Bin Packing: A Loss Model[END_REF]. In some special cases the constant V h can be expressed in closed form, while for remaining cases the computation time for a numerical evaluation of V h is negligible.

Our contributions are twofold:

• From a theoretical point of view, the main contribution is to prove the existence of the corrected term V h and to provide an algorithmic way to compute it (Theorem 3.1). This proves that it is possible to dene and compute the rened approximation h( ) + V h /N . • From a more practical perspective, we use a variety of examples that demonstrate that the rened approximation h( ) + V h /N is often highly accurate. To give a avor of the accuracy of the rened approximation, we provide some examples in Table 1 picked from Sections 4, 5 and 6. This table illustrates that the rened approximation is very accurate, even for N = 10. Moreover, it is much more accurate than the classic mean eld approximation. Our goal when presenting many examples is to show that the rened approximation is accurate in many scenarios. Also, each of the examples developed in Sections 4-7 serves a particular purpose.

Section 4 shows that it is sometimes possible to obtain a closed form formula for the rened term V h . Section 5 presents the supermarket model, which is one of the most widely studied mean eld models. We also use this example to show that our rened approximation is accurate enough to quantify 1/N -eects such as the impact of choosing with or without replacement among N servers. Section 6 shows that there are some cases for which the rened approximation is less than 1% from the actual value although the original mean eld approximation underestimate the response time by a factor 2. Section 7 presents an example in which the quantity of interest is a blocking probability. All the previous examples are continuous-time examples that fall into the scope of density-dependent population process. Appendix A shows that it is possible to apply our results to discrete-time models as well.

Our results are obtained by using Stein's method (see [START_REF] Braverman | Stein's method for steady-state diusion approximations: an introduction through the Erlang-A and Erlang-C models[END_REF] for a good introduction to Stein's method). The use of Stein's method to compute bounds for stationary distributions has been recently popularized in [START_REF] Braverman | Stein's method for steady-state diusion approximations of M /Ph/n + M systems[END_REF][START_REF] Gurvich | Diusion models and steady-state approximations for exponentially ergodic Markovian queues[END_REF]. This methodology has then been used in [START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF] and further extended in [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF][START_REF] Ying | Stein's Method for Mean Field Approximations in Light and Heavy Trac Regimes[END_REF] to establish the rate of convergence of stochastic processes to their mean eld approximation, in light or heavy trac. Our paper strongly relies on the methodology developed in those papers. In particular, it is shown in [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF] that the mean eld approximation is 1/N -accurate. We improve this result as we are able, for a large class of models, to express the constant V h by a simple formula that can be easily evaluated numerically. Note that the idea of using h( ) + V h /N as an improved approximation was already proposed in [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF] for the two-choice model, where the constant V h was estimated by simulation. The idea of improving classical approximations via Stein's method was also presented in [START_REF] Braverman | Stein's method for steady-state diusion approximations of M /Ph/n + M systems[END_REF] but, according to the authors, their computations seem hard to generalize to high dimensional settings. On the contrary, our method scales as the cube of the number of dimensions of the model which makes it applicable to models with hundreds of dimensions.

Roadmap. The rest of the paper is organized as follows. In Section 2 we introduce the model. We develop the theory in Section 3. We present several models in Sections 4 to 7 for which we compare the rened approximation with simulation. We conclude in Section 8.

MODEL AND NOTATIONS

Our result applies to many mean eld models, including density dependent population processes [START_REF] Thomas | Solutions of Ordinary Dierential Equations as Limits of Pure Jump Markov Processes[END_REF] or the classical model of [START_REF] Benaim | A class of mean eld interaction models for computer and communication systems[END_REF].

Density-Dependent Population Processes

In this paper, we consider a model that generalizes the notion of density-dependent population processes introduced by Kurtz in [START_REF] Thomas | Solutions of Ordinary Dierential Equations as Limits of Pure Jump Markov Processes[END_REF]. Consider a sequence of continuous-time Markov chains (CTMC) (X (N ) (t )) N . For each N , X (N ) (t ) evolves on a bounded subset E of R d . There exists a (nite or countable) set of vectors L 2 R d and a set of functions (N ) `: E ! R + such that X (N ) (t ) jumps from x to x + `/N at rate N (N ) `(x ) for each `2 L. We will refer to the chain X (N ) (t ) as the system of size N . For such a system, we dene the drift f (N ) as

f (N ) (x ) = X `2L ` (N ) `(x ).
The drift corresponds to the expected variation of a chain X (N ) (t ):

f (N ) (x ) = lim dt !0 1 dt E f X (N ) (t + dt ) X (N ) (t )|X (N ) (t ) = x g .
The classical notion of a density-dependent population process corresponds to the case where the transitions rates (N ) `(x ) do not depend on N .

Scaling Assumptions

In the following, we will prove our result under the assumptions (A0)-(A4) below. Note that Assumption (A2) ensures that from any initial condition x 2 E, the ODE ẋ = f (x ) has a unique solution. We will denote by t x the value of this solution at time t:

t x = x + Z t 0 f ( s x )ds.
The operator is called the semi-group of the dierential equation.

(A0) sup x 2E, N 1 P `2L (N ) `(x )k`k 2 < 1. (A1) There exist two functions f : E ! E and f : E ! E such that, uniformly in x: f (N ) (x ) = f (x ) + 1 N f (x ) + O ( 1 N 2
). (A2) f is twice-dierentiable and its second derivative is uniformly continuous and bounded. (A3) The ODE ẋ = f (x ) has a unique attractor and this attractor is exponentially stable, i.e.

there exists a, b > 0 such that for all x 2 E:

k t x k  ae bt .
(A4) For each N , the stochastic process has a unique stationary distribution. Note that for a function h : E ! R, we denote by E (N ) [h(X (N ) )] the expected value of h(X (N ) ) in steady-state.

The rst assumption (A0) is mostly technical and is needed so that the variance of the stochastic process does not explode in nite time. The second assumption is always true if the rate functions `do not depend on N (which is the case for the classical model of Kurtz [START_REF] Thomas | Solutions of Ordinary Dierential Equations as Limits of Pure Jump Markov Processes[END_REF]). It is needed to have a O (1/N ) convergence term and will be commented later.

Discrete-Time Model of [2]

In this paper, we mainly focus on mean eld interacting models that can be modeled as densitydependent population processes. Another type of model that is popular in the performance evaluation community is the mean eld interaction model of [START_REF] Benaim | A class of mean eld interaction models for computer and communication systems[END_REF]. In this section, we briey sketch how to adapt our result in this case.

The authors of [START_REF] Benaim | A class of mean eld interaction models for computer and communication systems[END_REF] consider a discrete-time model with N objects, where S N n (t ) 2 {1 . . . d } denotes the state of the object n at time slot t. The objects are assumed to be interchangeable, meaning that the empirical measure X (N ) is assumed to be a discrete-time Markov chain, where X (N ) i (t ) is the proportion of objects in state i at time t, dened as

X (N ) i (t ) = 1 N N X n=1 1 {S N n (t )=i } .
The main assumption of [START_REF] Benaim | A class of mean eld interaction models for computer and communication systems[END_REF] is that there exists a function (N ), called the intensity of the stochastic process, that vanishes as N goes to innity (i.e., lim N !1 (N ) = 0) and such that lim

N !1 (N )E f X (N ) (t + 1) X (N ) (t ) | X (N ) (t ) = x g = f (x ),
and such that the number of objects that change state during one time slot is bounded by (N ). It is then shown in [2, Theorem 1] that, the stochastic process X (N ) (t/ (N )) converges in probability to a deterministic quantity t x at rate p (N ), where t 7 ! t x is the solution of the dierential equation starting in x at time 0.

These discrete-time Markov chains can be transformed into continuous time Markov chains by assuming that each time-slot lasts for an exponential amount of time with mean 1/ (N ). This new chain has the same stationary distribution as the original discrete-time chain. Hence, as noted in [START_REF] Bortolussi | Continuous approximation of collective system behaviour: A tutorial[END_REF]Section 6] it is possible to transform most of the discrete-time models of this type into density-dependent population processes for which we can apply Theorem 3.1 to show that

E (N ) f h(X (N ) ) g = h( ) + (N )V h + o( (N )).
We present an example in Appendix A.

MAIN RESULTS

Computation of the Constant

We consider a mean eld model that has a unique xed point that is a global attractor. Let Q be the positive denite symmetric matrix whose (n, m)-th entry is equal to:

Q n,m = X `2L
`n`m `( ).

We denote by A the rst derivative of f in (i.e. the Jacobian in ) and by B the second derivative of f , that is, for i, j, k 1 , k 2 2 {1, . . . ,d }:

A i, j = @ f i @x j ( ), (B j ) k 1 ,k 2 = @ 2 f j @x k 1 @x k 2 ( ).
As indicated by Assumption (A3), we assume that the ODE system ẋ = f (x ) is exponentially stable, which implies that the matrix A is a Hurwitz matrix (meaning the real part of any eigenvalue of A is negative) [START_REF] Khalil | Nonlinear Systems[END_REF]. Therefore, by [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF], the Lyapunov equation AX + XA T + Q = 0 has a unique solution. We denote by W the unique symmetric matrix (as Q is symmetric) that satises

AW + W A T + Q = 0. (1)
Last, we dene by C the vector :

C = f ( ) = lim N !1 N ( f N ( ) f ( )).
The vector C is equal to zero if the transition rates do not depend on N (this is the case in all of our examples except for the coupon collector application in Section 4.1 and the supermarket without replacement of Section 5.3).

We are now ready to state the main theorem of the paper, whose proof is postponed to Section 3.4.

T 3.1. Assume that the model satises (A0-A4). Let h : E ! R be a twice-dierentiable function that has a uniformly continuous second derivative. Then,

lim N !1 N ⇣ E (N ) f h(X (N ) ) g h( ) ⌘ = X i @h @x i ( )V i + 1 2 X i, j @ 2 h @x i @x j ( )W i j , (2)
where the matrices A, C and W are dened above and V i is equal to:

V i = X j (A 1 ) i, j 2 6 6 6 6 6 4 C j + 1 2 X k 1 ,k 2 (B j ) k 1 ,k 2 W k 1 ,k 2 3 7 7 7 7 7 5 . (3)
In the following, we will mainly use this result to develop a rened approximation of the expectation of X (N ) i , which is asymptotically equal to i + V i /N + o(1/N ). We also note that this result can also be used to compute the co-variance matrix of the vector p N (X (N ) ). Both of these results are a direct consequence of Theorem 3.1 and are summarized in the following corollary.

C 1. Assume that the model satises (A0-A4). Then, using the same notations as in Theorem 3.1, we have:

(i) For any coordinate i, lim N !1 N ⇣ E (N ) f X (N ) i g i ⌘ = V i ;
(ii) For any coordinates i, j, the covariance satises:

lim N !1 N cov ⇣ X (N ) i , X (N ) j ⌘ = W i j .
P. The proof of (i) is a direct application of Theorem 3.1 with the function h(x ) = x i . For (ii), as i is deterministic, we have

N cov(X i , X j ) = N cov(X i i , X j j ) = N E (N ) f (X i i )(X j j ) g N (E (N ) [X i ] i )(E (N ) [X j ] j ) applying Theorem 3.1 with h(x ) = (x i i )(x j j )
shows that the rst term converges to W i j and the second is of order O (1/N ). ⇤

The above result shows that

E[X (N ) i ] is of order i + V i /N plus a o(1/N )-term.
In fact, a natural extension of Theorem 3.1 shows that, if the drift f is dierentiable 3 times, then this higher-order term is a O (1/N 2 )-term, and in that case we have :

E (N ) f X (N ) i g = i + 1 N V i + O ✓ 1 N 2 ◆ .
This result can be readily extended to higher-order dierentiable drifts.

T 3.2. Assume that the model satises Assumptions (A0-A4). If in addition the drift is continuously dierentiable k times, then for each coordinate i, there exists a series of constants V (2) i . . .V (k 1) i such that

E (N ) f X (N ) i g = i + 1 N V i + 1 N 2 V (2) i • • • + 1 N k 1 V (k 1) i + o ✓ 1 N k 1 ◆ ,
where V i is dened as in Theorem 3.1.

The result is an (almost direct) consequence of Lemma 3.5 that can be used to show that similarly to V i /N , the constant V (j ) i can be expressed as a function of some integral of the ow t and its derivative in up to the jth one. However, apart from particular cases, the computational cost of a numerical evaluation of these constants grows quickly with k. Hence, we believe that from a practical point of view, studying beyond k = 2 is not feasible. We provide some details about it in Appendix B.

Algorithmic Considerations

In some cases an explicit expression for the matrix W and the vector V in Theorem 3.1 can be obtained (e.g., see Section 4). However, in general the constant V i must be computed by a numerical procedure that consists of the following steps:

(1) Computation of the xed point and denition of Q.

(2) Computation of the rst and second derivative A and B.

(3) Computation of W by solving the Lyapunov Equation ( 1). (4) Application of (3) to compute the vector V . This procedure is a direct application of Theorem 3.1 and can be easily automated 1 . We further note that most of the conditions needed to apply the theorem can be checked automatically, at the exception of Assumption (A3): proving that a dynamical system has a unique attractor is undecidable [START_REF] Moore | Unpredictability and undecidability in dynamical systems[END_REF].

The time-complexity of the numerical procedure is essentially O (d 3 ), where d is the dimension of the mean eld model. The main step in the computation is the determination of the solution W to the Lyapunov equation AX + XA T + Q = 0. This solution can be computed in O (d 3 ) time using the Bartels-Stewart algorithm [START_REF] Bartels | Solution of the Matrix Equation AX + XB = C [F4][END_REF] (which is implemented by the MATLAB function lyap and by the function scipy.linalg.solve_lyapunov of the scipy package of python [START_REF] Jones | SciPy: Open source scientic tools for Python[END_REF]). The other steps of the algorithm are less costly. Because we assume that there exists a unique xed point that is exponentially stable, the xed point can be easily computed by numerically integrating the ODE and/or using a Newton-Raphson iteration close to . Once the xed point has been computed, the matrices A and B can be computed by using a symbolic toolbox. The computation of V i for all i can also be done in O (d 3 ).

In our numerical examples, we tested up to d = 500 for which the computation time was less than one second. We further note that in many cases the xed point can be obtained in closed form. Also, in some particular cases the entries of ( A) 1 can also be expressed in explicit form (e.g., see Section 5). Even in the cases when W and A 1 have to be computed numerically, the numerical computation of these constants is several orders of magnitude faster than simulation.

Extensions to Infinite-Dimensional Models and Unbounded State-Spaces

To simplify the exposition of the theory, we mostly restrict our attention to population processes that evolve in a subset of R d . This excludes innite-dimensional population models that arise naturally in queuing networks when the queues are unbounded like in the examples of Section 5 and Section 6 for which the state-space is a subset of [0, 1] Z + . The main reason to work with nite-dimensional models is that in many cases, we evaluate numerically the constants V and W by using numerical methods to invert A and to solve the Lyapunov equation ( 1). This requires the use of nite matrices. In our numerical examples, the only innite-dimensional models are the supermarket model of Section 5 and the push-pull model of Section 6. To obtain our numerical results, we truncated the models to models of dimension d = 40 for the supermarket model and d = 500 for the push-pull model. Higher values of d lead to the same rened constants V i up to the oating point errors.

From a theoretical point of view, Theorem 3.1 can be easily generalized to a Banach state space. One key issue, however, is that Assumption (A3) will in general not hold if the state-space E is not bounded, which is in general the case when E is innite-dimensional (in an unbounded space with bounded drift, the attractor cannot attract the trajectories in nite time). A possible way to overcome this diculty is to replace Assumption (A3) by the weaker Assumption (A5): which is to assume that (A3) holds on a bounded subset B ⇢ E and that the sequence of stationary measures of the stochastic processes concentrates on B at rate 1/N 3 : (A5) There exists a bounded set B ⇢ E such that:

(i) The ODE ẋ = f (x ) has a unique attractor and there exists a, b > 0 such that for all x 2 B:

k t x k  ae bt .
(ii) There exists a constant B > 0 such that in steady-state: P[X (N ) < B]  B/N 3 . Note that the only dierence between Assumption (A3) and Assumption (A5)-(i) is that the latter holds for x 2 B (and not for all x 2 E as in Assumption (A3)).

As for (A3), establishing (A5) is model-specic. For the supermarket and the pull-push model, Assumption (A5)-(i) was established in [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF] and [START_REF] Minnebo | A Fair Comparison of Pull and Push Strategies in Large Distributed Networks[END_REF]. Assumption (A5)-(ii) can be proved by using a coupling argument that shows that the processes is bounded by a system of N independent M/M/1 queues, see [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF]Section 6.4] for more details. This coupling argument seems to be easily generalizable to many models of resource allocation strategies.

We then obtain the following theorem:

T 3.3. Assume that the state-space E is a Banach space and that assumptions (A0,A1,A2,A4,A5) hold. Then, for any twice dierentiable bounded function h, Theorem 3.1 continues to hold and the constants V and W satisfy the same equations.

S . Let B be the bounded set of Assumption (A5). We proceed as for the proof of [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF]Theorem 3.2]. First, it should be noted that one may assume without loss of generality that h equals h( ) outside a bounded set. Then, as shown in [9, Equation ( 10)], we have:

N (E f h(X (N ) ) g h( )) =N E (N ) f ( L (N ) )Gh(X (N ) )1 {X (N ) <B } g + N E (N ) f ( L (N ) )Gh(X (N ) )1 {X (N ) 2B } g .
If h is bounded, then by assumption (A5)-(ii), the right-hand side of the rst line is of order O (1/N 2 ) (this can be shown by adapting the proof of [9, Lemma 6.4]). Hence, the result follows if one can show that the second term converges to Equation (2).

To show that, one needs to generalize Lemma 3.4, 3.5 and 3.6 to a Banach space. The essential dierence lies in the fact that the h-th derivative of function x 7 ! R 1 0 t xdt is bounded on any bounded set, but not necessarily on the whole state-space E (see [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF]Lemma 6.3]). This is why we need (A5). Apart from that, all the other arguments continue to hold -for instance, the uniqueness of the solution of the Lyapunov equation comes from the link between the exponential stability of the ODE and the Lyapunov equation in a Banach space that is established in [START_REF] Datko | Uniform asymptotic stability of evolutionary processes in a Banach space[END_REF][START_REF] Vu | On the Lyapunov equation in Banach spaces and applications to control problems[END_REF]. ⇤

Proof of Theorem 3.1

In what follows, for a function h : E ! R, we will denote by Dh the derivative of h (if it exists) and by D 2 h its second derivative. The notations Dh(x ) • `and D 2 h(x ) • (`, `) stand respectively for P i (@h(x )/@x i )`i and P i, j (@ 2 h(x )/(@x i @x j ))`i `j . For a multivariate function : E ! E, we also use the notations D and D 2 for the rst and second derivative. (D ) i j denotes @ i /@x j and (D 2 ) i,k 1 ,k 2 denotes @ 2 i /(@x k 1 @x k 2 ). P T 3.1. We denote by L (N ) the generator of the chain X (N ) (t) and by the generator of the ODE. They associate to a dierentiable function h : E ! R the functions

L (N ) h : E ! R and h : E ! R dened by (L (N ) h)(x ) = X `2L N (N ) `(x )(h(x + Ǹ ) h(x )), ( h)(x ) = Dh(x ) • f (x ) = X `2L `(x )Dh(x ) • `.
The proof can be divided in three steps. The rst step is to use Stein's method and the ideas developed in [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF][START_REF] Ying | On the Approximation Error of Mean-Field Models[END_REF] to show that, in steady-state, we have

E (N ) f h(X (N ) ) g h( ) = E (N ) " (L (N ) ) Z 1 0 (h( s X (N ) ) h( ))ds # . (4)
This is a direct consequence of Equation ( 10) of [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF].

The second step is to compare the generators L (N ) and to relate the above equation and the dierential of the function G :

x 7 ! R 1 0 ( s x
)ds, which is well-dened by (A3). This is done in Lemma 3.4.

Using perturbation theory, we then show in Lemma 3.5 that the function [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF], satises the assumptions of Lemma 3.4. This shows that:

G h : x 7 ! R 1 0 (h( s x ) h( ))ds, which is the solution of the Poisson equation ( G h )(x ) = h( ) h(x )
lim N !1 N E (N ) " (L (N ) ) Z 1 0 (h( s X (N ) ) h( ))ds # = DG h ( ) • C + 1 2 X `2L `( )D 2 G h ( ) • (`, `).
The proof is completed by Lemma 3.6 in which we derive the general formula presented in Equation [START_REF] Bortolussi | Continuous approximation of collective system behaviour: A tutorial[END_REF]. ⇤ L 3.4. Let be a twice dierentiable function that has a uniformly continuous second derivative that is bounded. Then, for x 2 E:

lim N !1 N (L (N ) ) x = D (x ) • f (x ) + X `2L `(x )D 2 (x ) • (`, `).
Moreover the convergence is uniform in x.

P. By denitions of the generators, we have

(L (N ) ) x = X `2L N (N ) `(x ) (x + Ǹ ) (x ) ! D (x ) • f (x ) = X `2L N ⇣ (N ) `(x ) `(x ) ⌘ (x + Ǹ ) (x ) ! + X `2L `(x ) N (x + Ǹ ) N (x ) D (x ) • `! . ( 5
)
We rst bound the rst summation in Equation ( 5). As has a bounded second derivative, there exists a constant c > 0 such that k

(x + Ǹ ) (x ) D (x ) • `/N k  c k`k 2 /N 2
. This shows that (uniformly in x), the rst summation in Equation ( 5) can be simplied to:

X `2L N ⇣ (N ) `(x ) `(x ) ⌘ (x + Ǹ ) (x ) ! = X `2L N ⇣ (N ) `(x ) `(x ) ⌘ D (x ) • Ǹ + O k`k 2 N 2 !! = (D (x ) + O ✓ 1 N ◆ ) • ( f N (x ) f (x )).
which is equal to

D (x ) • f (x )/N + O (1/N 2 ) by Assumption (A1).
For the second summation in Equation ( 5), as D 2 is uniformly continuous in x, the quantity

N 2 (x + Ǹ ) N 2 (x ) N D (x ) • `converges (uniformly in x) to 1 2 D 2 (x ) • (`, `
). As a consequence, the second summation in Equation ( 5) is equal to

1 2 P `2L `(x )D 2 (x ) • (`, `)/N + o(1/N ). ⇤ L 3.5. Let ( t ) : x 7 ! t
x be the ow associated to the drift f : E ! R and k 0. Assume that f and are k-times dierentiable with uniformly continuous derivatives. Assume that the ODE ẋ = f (x ) has a unique exponentially stable attractor (Assumption (A3)). Let h : E ! R be a k-times dierentiable function that has a uniformly continuous kth derivative. Then, the function

G h : x 7 ! R 1 0 (h( s x ) h( ))ds is k-times continuously dierentiable. Its kth derivative is bounded and equal to R 1 0 D k (h s )(x )ds.
P. According to [9, Lemma 6.3], there exists two constants b > 0 and > 0 such that for all t: kD t k  be t . According to [START_REF] Eldering | Normally Hyperbolic Invariant Manifolds-the Noncompact Case[END_REF]Lemma C.1], this implies that for any i  k + 1, the ith dierentiable of t exists and there exists a constant M i such that for all t: kD i t k  M i e t . Let G : x 7 ! R 1 0 ( s x )ds. As D i t converges exponentially fast to 0, the function G is i-times dierentiable and its ith dierentiable D i G is equal to

D i G (x ) = Z 1 0 D i t (x )ds.
By using the chain rule, the kth derivative of G h can be expressed as a function of the derivatives of order 1 to k of G and h (this is a multivariate generalization of the Faà di Bruno formula). ⇤ L 3.6. The following holds:

(i) The matrix equation AX + XA T + Q = 0 has a unique solution, that we denote by W . Moreover, for k 1 , k 2 , we have :

X `2L `( ) X n,m `n`m Z 1 s=0 (e As ) k 1 ,n (e As ) k 2 ,m ds = W k 1 ,k 2 . (6)
(ii) The derivative of G satises:

(DG ( )) i, j = ( A) 1 i, j , (7)
(iii) The second derivative of G h at satises:

X `2L `( )D 2 G h ( ) • (`, `) = X i, j @ 2 h @x i @x j ( )W i, j + X i @h @x i ( ) X j ( A) 1 i, j X k 1 ,k 2 (B j ) k 1 ,k 2 W k 1 ,k 2 . (8)
P. (i) -By denition of Q n,m , Q is a positive denite matrix. As A is Hurwitz, there exists a unique positive denite matrix W that satises AW + W A T + Q = 0 [START_REF] Simoncini | Computational methods for linear matrix equations[END_REF]. Moreover, this matrix W satises:

X n,m Z 1 s=0 (e As ) k 1 ,n Q n,m (e A T s ) m,k 2 ds = W k 1 ,k 2 . (9)
By denition, Q n,m = P `2L `( )`n`m. This implies that Equation ( 6) is equivalent to (9). (ii) -By the chain rule the derivative with respect to t of the Jacobian D ( t ) of the map t : x ! t x is given by

d dt (D( t )(x )) i, j = @ @x j ( f i t )(x ) = X u (D f ( t x )) i,u (D( t )(x )) u, j . As t = , we have D f ( t ) i,u = A i,u and therefore d dt (D ( t )( )) = AD ( t )( ). This implies that (D( t )( )) i, j = (e At ) i, j , (10)
which yields (7) as R 1 t =0 (e At ) i, j dt = ( A) 1 i, j due to assumption (A3). (iii) -By the chain rule, the derivative of the Hessian D 2 ( t ) of the map t is given by

d dt (D 2 ( t )(x )) j,nm = X u (D f ( t x )) j,u (D 2 ( t )(x )) u,nm + X k 1 ,k 2 (D 2 f ( t x )) j,k 1 k 2 (D( t )(x )) k 1 ,n (D( t )(x )) k 2 ,m .
Applying the above equality at implies that

d dt (D 2 ( t )( )) j,nm = X u A j,u (D 2 ( t )( )) u,nm + X k 1 ,k 2 B j,k 1 k 2 (e At ) k 1 ,n (e At ) k 2 ,m .
The variation of constants method shows that

(D 2 ( t )( )) i,nm = X j Z t s=0 (e A(t s ) ) i, j X k 1 ,k 2 B j,k 1 k 2 (e At ) k 1 ,n (e At ) k 2 ,m ds.
By switching the order of integration this implies that

Z 1 0 (D 2 t ( )) i,nm = X j ( A 1 ) i, j X k 1 ,k 2 B j,k 1 k 2 Z 1 s=0 (e As ) k 1 ,n (e As ) k 2 ,m ds. ( 11
)
By denition of G h , we have:

(D 2 (G h )( )) nm = X i @h @x i ( ) Z 1 0 (D 2 ( t )( )) i,nm dt + X i, j @ 2 h @x i @x j ( ) Z 1 0 (D( t )( )) i,n (D( t )( )) j,m dt .
The equality in (8) now follows by combining Equation ( 6), ( 10) and [START_REF] Gast | Markov chains with discontinuous drifts have dierential inclusion limits[END_REF] with the above equation. ⇤

APPLICATION 1 : COUPON REPLICATIONS

We consider the coupon replication model introduced in [START_REF] Massoulié | Coupon Replication Systems[END_REF] and [START_REF] Lin | Stochastic analysis of le-swarming systems[END_REF]. This application demonstrates that there are (exceptional) cases where an explicit expression can be obtained for the correction term V i of the rened mean eld model. In such a system there are K dierent coupons and users arrive with coupon i at rate i . Users are characterized by their current collection of coupons, and leave the system as soon as they collected the K coupons. Users regularly sample the population of users, either by sampling a user that currently holds the same number of coupons (i.e., the layered model) or by sampling any user (i.e., the at model). If the sampled user holds at least one coupon that is not owned by the sampling user, the latter adds such a coupon to his collection.

Layered Model

In this subsection we focus on a single layer of the layered model (see [START_REF] Massoulié | Coupon Replication Systems[END_REF]Section 3]). More specically we discuss the dynamics of the rst layer as the ODE of the other layers has exactly the same structure. Let K > 2 be the number of coupons and i the rate of users arriving in layer 1 with coupon i. Let x i (t ) be the rescaled number of users in the rst layer at time t that hold coupon i. The transitions of the system are (for i = 1, . . . , K):

x 7 ! x + 1 N e i at rate e i (x ) = i , x 7 ! x 1 N e i at rate e i (x ) = x i * , 1 x i P K k =1 x k + - .
The rst line corresponds to the arrival of a user with a coupon i and the second line to a user leaving to the next layer because he found a second coupon. The set of ODEs is therefore given by

d dt x i (t ) = f i (x ) = i x i (t ) * , 1 x i (t ) P K k=1 x k (t ) + - .
Let be the unique xed point of f . We have

f i (x + ) = i ( i + x i ) * , 1 ( i + x i ) P K k =1 ( k + x k ) + - for 1  i  K. Let = P K i=1 i , then we nd A i, j = @ f i @x j ( ) = ( ( i ) 2 / 2 i = j, 2 i / 2 i , j, (12)
and

(B j ) k 1 ,k 2 = @ f j @x k 1 @x k 2 ( ) = 8 > > > > > < > > > > > : 2 2 j /( 3 ) j = k 1 = k 2 , 2( j ) 2 /( 3 ) j , k 1 , k 2 , 2 j ( j )/( 3 ) j = k 1 , k 2 , 2 j ( j )/( 3 ) j = k 2 , k 1 .
Symmetric case: From now on we focus on the symmetric case i = for all i 2 {1, . . . , K } (which was proven to minimize the mean sojourn times in [START_REF] Massoulié | Coupon Replication Systems[END_REF]). In this case i = K/(K 1), meaning the number of users in the system with coupon i is estimated by the classic mean eld limit as N K/(K 1). Due to [START_REF] Gurvich | Diusion models and steady-state approximations for exponentially ergodic Markovian queues[END_REF], A is given by

A = 1 K 2 ⇣ K (K 2)I + ee T ⌘ ,
where I is the identity matrix and e a column vector of ones. As A is the sum of a diagonal matrix and a rank one matrix, the Sherman-Morrison-Woodbury formula implies

(A 1 ) i, j = 1 (K 1)(K 2) 1 {i=j } K K 2 ,
where 1 {Z } = 1 if Z is true and 1 {Z } = 0 otherwise. The matrix Q is a diagonal matrix with all its entries equal to 2 . As A is symmetric the unique solution W of the Lyapunov equation AX + XA T + Q = 0 is given by A 1 . Some basic algebra shows that X

k 1 ,k 2 (B j ) k 1 ,k 2 W k 1 ,k 2 = 2(K 1) 2 (K 2)K 2 ,
and 1 2

X j (A 1 ) i, j X k 1 ,k 2 (B j ) k 1 ,k 2 W k 1 ,k 2 = (K 1) K (K 2) .
This implies that the mean eld correction term given by (3) has the following simple explicit form

V i = (K 1) K (K 2) ,
and the rened mean eld model suggest that the mean number of users that hold coupon i in the symmetric system with i = N is approximately In Table 2 we present a comparison of the classic and rened mean eld result with simulation. For this simulation users select a user uniformly at random after an exponential amount of time with mean one and we allow users to select themselves. We observe that the rened model is very close to the result obtained by simulation. For the system where we demand that the selected user is another user, the rate of the transitions with `= e i becomes dependent on N and one nds lim

N K K 1 + K 1 K (K 2) . K N i N i N + V i simulation
N !1 N ( (N ) e i ( ) e i ( )) = i P k ,i k ⇣ P K k =1 k ⌘ 2 ,
which equals (K 1)/K2 in the symmetric case. As the row sum of A 1 equals K/(K 1) we get the following rened mean eld approximation for the number of users in the system holding coupon i in case a user always selects another user

N K K 1 + K 1 K (K 2) 1 K = N K K 1 + 1 K (K 2)
, meaning the correction term with respect to the classic mean eld reduces by a factor K 1.

Flat Model

We now consider the at model considered in [19, Section 3.1] with m = 1, which is a reformulation of the model presented in [START_REF] Massoulié | Coupon Replication Systems[END_REF]Section 4]. Let K > 2 be the number of coupons that a user needs to collect and the rate of users arriving with a single coupon. Let i (t ) be the fraction of the total number of users in the system at time t that hold i coupons. Dene

q i ( ) = K 1 X j=1 p i, j j P K 1 u=1 u
, where p i, j = 1 if j > i and p i, j = 1 C i j /C K j , for j  i, with C b a the number of ways to select a dierent items from a set of b dierent items.

The density dependent Markov model contains the following transitions: e 1 ( ) = , e i +1 e i ( ) = i q i ( ), for i = 1, . . . , K 1, and e K 1 ( ) = K 1 q K 1 ( ). Hence, the set of ODEs is given by

d dt i (t ) = i 1 (t )q i 1 ( (t )) i (t )q i ( (t )),
for i = 2, . . . , K 1 and d dt 1 (t ) = 1 (t )q 1 ( (t )). We compute a xed point 2 of this set of ODEs using [START_REF] Massoulié | Coupon Replication Systems[END_REF]Proposition 3] and note that for any xed point we have = i q i ( ), for

N 1 N 1 N + V 1 simulation 9 N 9 N + V 9 simulation
A i, j = 1 {i >1} p i 1, j i 1 p i, j i + 1 {i=1} 1 {i=j } i + 1 {i 1=j } i 1
, and

(B j ) k,s =(p j,k + p j,s ) j 2 (p j 1,k + p j 1,s ) j 1 2 1 {j >1} 1 {j=1} 2 2 + 1 {k =j } j p j,s ! + 1 {s=j } j p j,k ! 1 {k =j 1} j 1 p j 1,s ! 1 {s=j 1} j 1 p j 1,k ! .
Note that for j > i > 1 we have A i, j = ( i 1 i )/ and for k, s > j > 1 we have (B j ) k,s = 2( j j 1 )/ 2 . Further, the entries of A are independent of and the ones of B are inversely proportional to . As Q is a tridiagonal matrix with Q n,n = 2 and Q n,n+1 = Q n,n 1 = , we note that the entries of W are proportional to and therefore the entries V i are independent of as in the layered model.

Table 3 illustrates the accuracy of the rened mean eld model in case K = 10 for i = 1 and i = K 1 = 9. Similar to the layered model we observe that the rened model is much more accurate that the classic mean eld limit especially for small N .

APPLICATION 2 : THE SUPERMARKET MODEL

In this section, we consider the supermarket model introduced in [START_REF] David | The Power of Two Random Choices in Randomized Load Balancing[END_REF][START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF] where each job is allocated a server according to the "power-of-two-choices" algorithm. This model -or variants of this model -have been heavily studied in the literature and often serve as a basic example of a mean eld interaction system. In this section, we show how our theory can be used to rene the analysis of this model. We also show that our rened approximation is accurate enough to dierentiate between choices with replacement and without replacement.

The supermarket model

We consider a system composed of N identical servers. Jobs arrive at the system according to a Poisson process of rate N . The service time of each job is exponentially distributed of mean 1. For each incoming job, two servers are picked at random. The jobs is allocated to the server with the smallest queue size. Let X (N ) i (t ) denote the fraction of servers with queue size at least i at time t. As the role of each server is symmetrical, X (N ) (t ) is a Markov chain whose transitions are as follows. A departure from a server with i 1 jobs modies X into X N 1 e i and occurs at rate N (X i X i+1 ). An arrival at a queue with i jobs modies X into X + N 1 e i . The least loaded among the two has i 1 jobs if both servers have more than i 1 jobs but not both have i jobs. This probability depends on whether the two servers are picked with replacement or without replacement:

• If one picks two servers with replacement, this occurs with probability X 2 i 1

X 2 i .
• If one picks two servers without replacement, this occurs with probability X i 1 (NX i 1 1)/(N 1) X i (NX i 1)/(N 1)

This leads to two (slightly) dierent models. For the classical two-choice model (with replacement), the transitions are, for i 1:

X ! X 1 N e i at rate N (X i X i+1 ) X ! X + 1 N e i at rate N (X 2 i 1 X 2 i ) (13)
The transition of the model without replacement are:

X ! X 1 N e i at rate N (X i X i+1 ) X ! X + 1 N e i at rate N (X i 1 N X i 1 1 N 1 X i N X i 1 N 1 ) (14)
The drift of the classical model (Equation ( 13)) does not depend on N and its ith coordinate is:

f i (x ) = (x 2 i 1 x 2 i ) + (x i+1 x i ). ( 15
)
This drift is also the limit of the model without replacement (Equation ( 14)) which shows that the mean eld approximation of both models coincides as N goes to innity. This explains why in the literature, there is no distinction made between the two models : the classical convergence result cannot dierentiate the two models. In fact, the error between Equation ( 13) and Equation ( 14) is of order O (1/N ), which is of the same order as our rened term of Theorem 3.1. We will explore this in more detail in Section 5.3 to show that our results allows one to study the impact of choosing with or without replacement.

Refined Model for the Classical Two-Choice

In this section, we rst study the rened model for the classical two-choice model (i.e., with replacement Equation ( 13)). Its drift is given by Equation [START_REF] Thomas | Solutions of Ordinary Dierential Equations as Limits of Pure Jump Markov Processes[END_REF]. It has a unique xed point = (1, 1 , 2 , . . . ) given by i = 2 i 1 .

5.2.1

Computations of A, B and W . To apply Theorem 3.1, the rst step is to compute the Jacobian and Hessian of f . To compute the Jacobian and Hessian of f , let us denote by i = x i i . Developing f ( + ), when ! 0, we obtain for i 1: (note: some terms cancel because f ( ) = 0)

f i ( + ) = 2 ( i 1 i 1 i i ) + i+1 i + ( 2 i 1 2 i ) = 2 i 1 i 1 (2 i + 1) i + i+1 + ( 2 i 1 2 i ).
The Jacobian and Hessian in are (for i, j, k 1 , k 2 1):

A i, j = @ f i @x j (0) = 8 > > > > > < > > > > > : 2 i 1 = 2 2 i 1 if j = i 1 2 i 1 = 2 2 i 1 if j = i 1 if j = i + 1 0 otherwise (16)
(B j ) k 1 ,k 2 = @ 2 f j @x k 1 @x k 2 (0) = 8 > > < > > : 2 if k 1 = k 2 = j 2 if k 1 = k 2 = j 1 0 otherwise (17)
The matrix Q of Equation ( 1) is:

Q = P `2L Q `with (Q `)n,m = `n`m `( ).
For the 2-choice model, this matrix is a diagonal matrix whose (i, i)-entry is

Q i,i = ( e i ( ) + e i ( )) = 2( 2 i 1 2 i +1 1 )
Let W be the (unique) solution of the Lyapunov equation AX + XA T + Q = 0. The constant V i of Theorem 3.1 is (for i 1):

V i = 1 2 X j (A 1 ) i, j X k 1 ,k 2 (B j ) k 1 ,k 1 W k 1 ,k 2 , = 1 2 X j A 1 i, j 2 ( W j, j + W j 1, j 1 ) = X j (A 1 i, j A 1 i, j+1 )W j, j (18)
It is not hard to verify that the general term of

A 1 i, j A 1 i, j+1 is A 1 i, j A 1 i, j+1 = ( 2 i j 1 2 i 2 j+1 if j < i 0 otherwise.
It is unclear whether the matrix W can be computed in closed form.

Average queue length.

In Table 4, we display the average queue length, computed by simulation or computed by using the mean eld approximation plus

P i V i /N
, where V i is the constant given by Equation [START_REF] Boudec | A generic mean eld convergence result for systems of interacting objects[END_REF]. Again, we observe that the rened mean eld approximation greatly improves the quality of the approximation. Even for N = 10, the rened model has a relative error below 5%. Moreover, as the error of rened-mean eld decreases as O (1/N 2 ), the rened model quickly becomes very accurate. We conrm the observation made in [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF] that the error made by the classical mean eld approximation seems to grow as 1/( 1) as the load approaches 1. In [START_REF] Gast | Expected values estimated via mean-eld approximation are 1/N-accurate[END_REF], it was conjectured that the constant P i V i /N was close (but not equal to) 2 /(2 2 ). Our numerical evaluation conrms this tendency.

N = 10 N = 20 N = 30 N = 50 N = 100 N = 1 Simulation ( = 0.7)
The results in Table 4 also suggest that the dierence between the mean eld approximation and the values obtained by simulation again grows quickly as approaches 1. By Theorem 3.2, it is possible to develop a second order approximation of the form + V /N + V (2) /N 2 . By using a numerical tting of the parameter, we propose to approximate the average queue length by the quantity a(N , ), dened by

a(N , ) = 1 X i=1 i ( ) + 2 2N (1 ) + 1 20N 2 (1 ) 2 . ( 19
)
The above formula suggests that error of the rened approximation grows as 1/(1 ) 2 as approaches 1.

In Table 5, we provide some numerical results showing the accuracy of this second order approximation against simulation result. This new approximation appears to be very accurate for all tested values. N = 10 N = 20 N = 30 N = 50 N = 100 Simulation ( = 0.7)

Distribution of queue size.

In Figure 1, we plot the marginal law of the number of jobs at a server. We compare the mean eld approximation (dashed line), the simulation results and the rened mean eld approximation, given by i + V i /N , where V i is given by Equation [START_REF] Boudec | A generic mean eld convergence result for systems of interacting objects[END_REF]. We restrict our observation to N = 50 and = 0.95, but other values are similar.

In Figure 1(a), we observe that the rened approximation predicts the shape of the distribution with great accuracy : for N = 50, the curves are almost indistinguishable. The situation of the tail is dierent. In Figure 1(b), we plot in log-scale the quantity E[X (N ) i ] = P[queue length i] as a function of i. We observe that while the rened mean eld approximation greatly improves the accuracy for the average response time or the shape of the distribution, it does not accurately predict the tail behavior: it is an improvement compared to the mean eld approximation, but remains quite inaccurate. This can be explained as follows. Theorem 3.1 guarantees that, uniformly in i,

i + V i N P [queue length i] = O ✓ 1 N 2 ◆ .
However, the relative error i + V i /N P[queue length i] /P[queue length i] is not guaranteed by this result to remain small as i grows because the quantity P[queue length i] converges quickly to 0 as i grows.

Impact of Choosing Without Replacement

We now focus on the comparison of the two-choice model with replacement (Equation ( 13)) and without replacement (Equation ( 14)). The only dierence between the two models is the rate of arrivals. Rewriting Equation [START_REF] Jones | SciPy: Open source scientic tools for Python[END_REF] and Equation ( 14), the rate of a transition X 7 ! X + e i /N (for i 1) is one of the following:

i (x ) = (x 2 i 1
x 2 i ) (with replacement, Eq. ( 13))

(N ) i (x ) = (x i 1 x i )(x i 1 + x i 1 N ) N N 1
(without replacement, Eq. ( 14))

For each i, the function (N ) i (x ) converges to the function i (x ) as N goes to innity. Moreover, we have

lim N !1 N ( (N ) i (x ) i (x )) = (x i 1 x i )(x i 1 + x i 1).
(a) queue-length distribution (b) tail behavior (log-scale). Theorem 3.1 shows that, for a coordinate i, the expected number of servers with a queue size i or more, E[X i ], can be expressed as

E f X (N ) i g = i + V i N + E i N + O ✓ 1 N 2 ◆ , (20)
where V i is the rened constant of the classical two-choice (with replacement), dened in Equation [START_REF] Boudec | A generic mean eld convergence result for systems of interacting objects[END_REF] and E i is dened by :

E i = 1 X j=1 A 1 i, j ( j 1 j )( j 1 + j 1) = 1 X j=1 ( 2 j j )(A 1 i, j+1 A 1 i, j ) = i 1 X j=1 (1 j )2 i j 1 2 i 2 j 5.
3.1 Average queue length. For the average queue length, the dierence is the sum over all i of the above expression. This shows that for large N , the average queue length of the two-choice model without replacement will be E/N smaller than the one of the two-choice model with replacement, where

E = 1 X i=1 E i = 1 X i=1 i 1 X j=1 (1 j )2 i j 1 2 i 2 j = 1 X j=1 1 X i=1 ( 2 i +j 2 j 2 i +j 1 )2 i 1
As an illustration, we compare in Table 6 the rened approximation and values obtained by simulation for various values of and N . We observe that in most cases, the error E/N is a good predictor of the impact of choosing with or without replacement. can be expressed as

P i i + V /N + E/N + O (1/N 2 )
where P i i is the mean eld approximation and P i i + V /N is the rened approximation for the classical model with replacement. In Figure 2, we compare the constant V and E by plotting them as a function of . We plot, in solid blue the rened constant V and in dashed orange the term E that corresponds to choosing with replacement and without replacement. We observe that the constant V is always larger than E. Moreover, when is close to 1, the error due to with or without replacement grows, but grows much slower than the correction term V . Numerically, we observe that E seems to grow as log(1/( 1)) as goes to one (similar to the average queue length) whereas V seems to grow as 1/( 1).

Fig. 2. Two-choice : Comparison of the constants V (error of the mean field) and E (dierence between choosing with or without replacement).

APPLICATION 3 : PUSH-PULL STRATEGIES

In this section we illustrate the rened mean eld approximation on a model used to study a class of pull and push strategies in a distributed network consisting of N servers. The main reason for considering this application is that it illustrates that the rened mean eld model can drastically improve the accuracy of the classic mean eld model where for N = 10 errors close to 100% reduce to less than 1% (see Table 7).

The push-pull model is as follows. Jobs arrive at each server according to a Poisson process. Each server has a single server with an innite waiting room and serves jobs in FCFS order. Under a pull strategy idle servers sample a server uniformly at random at rate r and if the sampled server holds at least 2 jobs, a job is exchanged. In case of a push strategy the servers with pending jobs sample servers in order to nd an idle server. As shown in [START_REF] Minnebo | A Fair Comparison of Pull and Push Strategies in Large Distributed Networks[END_REF], the performance of both the pull and push strategy can be captured using the mean eld model considered in this section.

Let x k be the fraction of the servers with k or more jobs. There are three types of transitions (arrivals, departures or job transfers):

e k (x ) = (x k 1 x k ), e k (x ) = (x k x k +1 ), e 1 e k (x ) = r (1 x 1 )(x k x k+1 ).
One easily veries that (as is a xed point)

f 1 (x + ) = (1 + + r 2 )x 1 + (1 + r (1 )r )x 2 rx 1 x 2 , f i (x + ) = x i 1 ( + 1 + (1 )r )x i + (1 + (1 )r )x i+1 + ( i i+1 )rx 1 + rx 1 (x i x i+1 ), for i > 1. Dene = /(1 + (1 )r ), such that i = i 1 for i 1, then A i, j = 8 > > > > > > > > > > > > < > > > > > > > > > > > > : 2 r / i = j = 1, + (1 ) 2 r i = 2, j = 1, (1 ) i r i> 2, j = 1, 1 j = i + 1, ( 1
+ ) j = i > 1, j = i 1 > 0, 0 otherwise, and
(B j ) k 1 ,k 2 = 8 > > < > > : r j 1, (k 1 , k 2 ) = (1, j + 1) or (j + 1, 1) r j > 1, (k 1 , k 2 ) = (1, j) or (j, 1) 0 otherwise.
Note that the transposed matrix A T is the subgenerator matrix of a phase-type distribution as its diagonal entries are negative, its o-diagonal entries are nonnegative and its row sums are all equal to zero, except for the rst row which sums to minus one [START_REF] Latouche | Introduction to Matrix Analytic Methods and stochastic modeling[END_REF]. Hence, the eigenvalues of A all have negative real parts and the Lyapunov equation AX + XA T + Q = 0 has a unique solution. 7. Mean queue length under pull/push with r = 1/(1 ): simulation vs refined mean field approximation Fig. 3. eue length distribution under pull/push with r = 1/(1 ): simulation vs refined mf. approximation.

APPLICATION 4 : BIN PACKING AND BLOCKING PROBABILITIES

In this subsection we consider the bin packing model of [START_REF] Xie | Power of D Choices for Large-Scale Bin Packing: A Loss Model[END_REF], where the total job arrival rate is N and each server is a loss queue with capacity K. Arriving jobs pick 2 servers at random and join the queue with the largest spare capacity. If both selected queues have no spare capacity the incoming job is blocked and considered lost. While the focus was mainly on the mean response time in case of the previous applications, we now look at the accuracy of the rened mean eld model when considering small blocking probabilities.

There are two types of transitions (arrivals and departures):

e k (x ) = (x 2 k 1 x 2 k ), e k (x ) = k (x k x k +1 ),
for k = 1, . . . , K, where x 0 = 1 and x K +1 = 0. eld model is signicantly closer. Note that, as for the supermarket and the pull/push model, the rened mean eld is an improvement compared to the classical mean eld but does not fully predict the tail behavior: this explains why the relative error is larger when K = 40 than when K = 10. The blocking probability p blocked of the rened mean eld model was computed by noting that the rate at which jobs are processed should match the rate at which jobs are accepted into the system, i.e.,

(1 p blocked ) = K X i=1 ( i + V i /N ).
For the classic mean eld approximation this is equivalent to stating p blocked = 2 K as the queue lengths are asymptotically independent, but this is not the case for nite N . This is further illustrated by Figure 4, which compares the server occupancy distribution of the rened mean eld approximation with simulation for = 9 and K = 10. It shows that while the probability of having K jobs tends to decrease as we decrease N , the blocking probability nevertheless increases (see Table 8) due to the correlation in the joint queue lengths.

CONCLUSION

In this paper, we proposed a rened mean eld approximation that is shown to be signicantly more accurate than the classic mean eld limit for systems composed of a limited number of objects. This renement exists in approximating E (N ) [h(X (N )] by h( ) + V h /N , where is the xed point (i.e., h( ) is the classic mean eld limit) and V h is the limit of N (E (N ) [h(X (N )] h( )) as N tends to innity. We showed that the constant V h can be expressed via the Jacobian and Hessian of the drift f in and the solution of a single Lyapunov equation, which can be evaluated in less than a second for models with as many as 500 dimensions. This gives the following drift equations, where K = d 1:

f 0 (x ) = K 1 X k =1 q k x k (1 (x )) q 0 x 0 (x ) + q K x K = K X k =0
q k x k (1 (x )) q 0 x 0 + q K x K (x ),

f k (x ) = (x )q k 1 x k 1 q k x k , for k 1.
which is the ODE presented in [START_REF] Cho | On the asymptotic validity of the decoupling assumption for analyzing 802[END_REF]]. Its xed point is given by the famous Bianchi's formula:

x k = k q k x 0 = k q k P K j=0 j /q j ,
where = 1 e P k q k x k . To determine the rened mean eld approximation, the rst and the second derivative of the drift can be easily obtained using a symbolic computation toolbox like SymPy. Also, we note that this model has a correction term C due to the limit lim N !1 (N ) k = .

B HIGHER ORDER TERMS

In this section, we sketch how the constants D (2) . . . D (k ) can be computed. By using Lemma 3.5, the proof of Lemma 3.4 can be adapted to show that, for a function h that is k + 1 times dierentiable, we have: As for Equation [START_REF] Gast | Markov chains with discontinuous drifts have dierential inclusion limits[END_REF], by using the chain rule, the expression R 1 0 D k t ( )dt can be expressed as an integral of functions that involve up to the kth derivatives of the drift f evaluated in as well as some factors of the form Z 1 0 exp(At ) i 1 , j 1 exp(At ) i 2 , j 2 . . . exp(At ) i k , j k dt .

(L (N ) h)(x ) = ( h)(x ) + 1 2N X `2L `(x )D 2 h(x ) • (`, `) + • • • + 1 k!N k 1
B.1 Is a numerical evaluation possible?

To evaluate the expression R 1 0 D k t ( )dt, one suers two drawbacks. First, the complexity of this expression grows quickly with k, and second it requires to generalize the Lyapunov equation (1) in order to evaluate numerically the integral [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF].

In what follows, we show how to evaluate numerically this integral by using the Kronecker product . Let M be the (d k ) ⇥ (d k ) matrix dened by :

M = A A • • • A.
(k times)

In terms of indices, for k = 3, this means that M i 1 i 2 i 3 , j 1 j 2 j 3 = A i 1 j 1 i 2 j 2 i 3 j 3 + i 1 j 1 A i 2 j 2 i 3 j 3 + i 1 j 1 i 2 j 2 A i 3 j 3 .

We have:

exp(Mt ) = exp(At ) ⌦ exp(At ) ⌦ • • • ⌦ exp(At ). (k times)
Again, in terms of indices, for k = 3, this means that: exp(Mt ) i 1 i 2 i 3 , j 1 j 2 j 3 = exp(At ) i 1 , j 1 exp(At ) i 2 , j 2 exp(At ) i 3 , j 3 .

As A is Hurwitz, so is M. This implies that R 1 0 exp(Mt )dt = M 1 . In particular, this shows that Z 1 0 exp(At ) i 1 , j 1 exp(At ) i 2 , j 2 exp(At ) i 3 , j 3 dt = (M 1 ) i 1 i 2 i 3 , j 1 j 2 j 3 .

As a result, the terms of Equation ( 22) can be computed by inverting the d k ⇥ d k -matrix M which can be done in O (d 3k ). Note that for k = 2, this method is less ecient than solving the Lyapunov equation (which can be done in O (d 3 )).

ACKNOWLEDGMENT

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

Fig. 1 .

 1 Fig.1. Two-choice model : Distribution of the queue length estimated by simulation, the mean field approximation and its refined version.

5. 3 . 2

 32 Error of the mean field approximation v.s. error of the model with/without replacement. As indicated by Equation (20), the average queue length for the two-choice model without replacement N = 10 N = 20 N = 50 = 0.

Fig. 4 .

 4 Fig. 4. Bin packing: Server occupancy distribution or = 9 and K = 10: simulation vs refined mf. approximation.

  )D k h(x ) • (`, . . . , `)

Table 2 .

 2 Mean number of users with coupon i in the symmetric layered model with arrival rate i N = N /10

	5 10 1.250	1.517	1.530
	20 2.500	2.767	2.772
	50 6.250	6.517	6.522
	8 10 1.143	1.289	1.291
	20 2.286	2.432	2.435
	50 5.714	5.860	5.861

Table 3 .

 3 Mean number of users with i = 1 and i = 9 coupons in the flat model with arrival rate = N /10 and K = 10 i = 1, . . . , K 1. Using this equality and denoting P i i = one nds (for any xed point)

	10 1.009	1.018	1.019	1.839	1.752	1.752
	20 2.019	2.028	2.024	3.678	3.591	3.588
	50 5.047	5.056	5.062	9.194	9.107	9.100

Table 4 .

 4 Two-choice model : Average queue length for various values of and N . We compare simulation with the refined mean field approximation

	1.2194 1.1735 1.1584 1.1471 1.1384	-
	Rened mean eld 1.2150 1.1726 1.1584 1.1471 1.1386 1.1301
	Simulation ( = 0.9) 2.8040 2.5665 2.4907 2.4344 2.3931	-
	Rened mean eld 2.7513 2.5520 2.4855 2.4324 2.3925 2.3527
	Simulation ( = 0.95) 4.2952 3.7160 3.5348 3.4002 3.3047	-
	Rened mean eld 4.1017 3.6578 3.5098 3.3915 3.3027 3.2139

Table 5 .

 5 Two-choice model : Average queue length for various values of and N . We compare simulation with a second order approximation defined by Equation[START_REF] Lin | Stochastic analysis of le-swarming systems[END_REF].

		1.2194 1.1735 1.1584 1.1471 1.1384
	a(N , )	1.2173 1.1723 1.1580 1.1467 1.1383
	Simulation ( = 0.9) 2.8040 2.5665 2.4907 2.4344 2.3931
	a(N , )	2.8077 2.5677 2.4932 2.4357 2.3937
	Simulation ( = 0.95) 4.2952 3.716 3.5348 3.4002 3.3047
	a(N , )	4.3164 3.7151 3.5369 3.4024 3.3061

Table 6 .

 6 Two-choice : comparison with and without replacement. The values in parenthesis (rmf=•) correspond to the refined mean field approximation, the others to the values obtained by simulation. The rows "withwithout" are the dierence between "with replacement" and "without replacement".

Table

	0.80	N Simulation	10 1.5569 1.4438 1.3761 1.3545 20 50 100	1 -
		Rened mean eld 1.5473 1.4403 1.3761 1.3547 1.3333
	0.90	Simulation	2.3043 1.9700 1.7681 1.7023	-
		Rened mean eld 2.2945 1.9654 1.7680 1.7022 1.6364
	0.95	Simulation	3.4288 2.6151 2.1330 1.9720	-
		Rened mean eld 3.4369 2.6232 2.1350 1.9723 1.8095

Proceedings of the ACM on Measurement and Analysis of Computing Systems, Vol. 1, No. 2, Article 33. Publication date: December 2017.

A prototype of a tool that implements this algorithm is available at https://github.com/ngast/rmf_tool/.Proceedings of the ACM on Measurement and Analysis of Computing Systems, Vol. 1, No.

2, Article 33. Publication date: December 2017.

Global attraction of this xed point is still an open issue. Proceedings of the ACM on Measurement and Analysis of Computing Systems, Vol. 1, No. 2, Article

Publication date: December 2017.

The sum formula (3) therefore simplies to

(A 1 ) i, j W 1, j + / -= r * . ,

as W is symmetric. We now argue that (A 1 ) i, j (A 1 ) i, j+1 = ( i j / i > j,

A T is the subgenerator matrix of a phase-type distribution, as such entry (j, i) of ( A T ) 1 is equal to the expected time spent in phase i given that we start in state j until absorption. When i  j it is clear that the expected time spent in phase i remains the same whether we start from state j or j + 1 (as any path from phase j + 1 to phase i must visit phase j rst). Therefore (A 1 ) i, j (A 1 ) i, j+1 = 0 for i  j. When i > j we nd that (A 1 ) i, j (A 1 ) i, j+1 = (( A T ) 1 ) j+1,i (( A T ) 1 ) j,i which is the expected time spent in phase i starting from phase j + 1 under taboo of phase j. This expectation is given by i (j+1) /( / ) due to the birth-death structure [START_REF] Latouche | Introduction to Matrix Analytic Methods and stochastic modeling[END_REF].

Due to [START_REF] Minnebo | A Fair Comparison of Pull and Push Strategies in Large Distributed Networks[END_REF] we obtain

with W the unique solution of the Lyapunov equation

A comparison between the mean queue length of the rened mean eld and simulation is presented in Table 7 when r = 1/(1 ) for various N and values. We remark that for = 0.95, the mean eld approximation has relative errors close to 90%, while the rened model has relative errors below 0.5%.

Figure 3 illustrates the accuracy of the rened mean eld approximation on the queue length distribution when = 0.90. Simulation results are based on a run of length 10 8 . As for the supermarket model (see Figure 1(b)), we observe again that although the rened approximation greatly improves the accuracy of the estimation of the average response time compared to the classical mean eld approximation, it is less accurate for the tail behavior. The explanation is similar to that of §5.2.3: in Figure 3, we observe that the relative error grows with i because it is of the order of O (1/N )/P[queue length i], but the absolute error is bounded uniformly in i. It is shown in [START_REF] Xie | Power of D Choices for Large-Scale Bin Packing: A Loss Model[END_REF] that the corresponding system of ODE has a unique xed point . Moreover, around this point we have

The sum formula (3) can therefore be written as

with W the unique solution to AX + XA T + Q = 0 where Q is a diagonal matrix with

), for i = 1, . . . , K, with 0 = 1 and K +1 = 0.

Table 8 illustrates the accuracy of the rened mean eld model for various parameter settings for and K for the blocking probability. While the blocking probability of the classic mean eld approximation can be several orders of magnitude below the simulation result, the rened mean One of the key assumptions to establish our main result is that the drift f is twice-dierentiable. For systems with a discontinuous drift (e.g., [START_REF] John | On the power of (even a little) centralization in distributed processing[END_REF]) the decay of E (N ) [h(X (N )] h( ) may be slower than 1/N if the xed point is located at a point in which the drift is discontinuous. When the drift is up to k times dierentiable, the rened approximation can in principle be further improved to obtain a o(1/N k ) error. However, at this stage it is not clear whether an ecient computational method can be devised to compute the constants of the 1/N i terms for i = 2, . . . , k.

APPENDIX A DISCRETE-TIME MODEL : A 802.11 EXAMPLE

In this section, we show how we can adapt our results to the model of the 802.11 MAC protocol presented in [START_REF] Cho | On the asymptotic validity of the decoupling assumption for analyzing 802[END_REF]. This models falls into the class of mean eld interacting models of [START_REF] Benaim | A class of mean eld interaction models for computer and communication systems[END_REF] and the purpose of this section is to show that one can easily apply our results to such models.

The authors of [START_REF] Cho | On the asymptotic validity of the decoupling assumption for analyzing 802[END_REF] presents a discrete-time model of the 802.11 MAC protocol. There are N transmitters and each transmitter can be in a state {0, . . . ,d 1}. At each time slot, a transmitter that is in a state k attempts to transmit with probability q k /N . If it is the only one transmitting, then the transmission is successful and its state becomes 0 at the next time slot. If there are two or more nodes transmitting, then the state becomes (k + 1) mod d.

Let X (N ) i (t ) denote the fraction of transmitters in state i at time t. If n is a node in state k, then the probability that at least one other node transmits at the same time is k (X (N ) (t )), where:

where Nx s is the number of nodes in state s. This model falls into the class of models of [START_REF] Cho | On the asymptotic validity of the decoupling assumption for analyzing 802[END_REF] with an intensity (N ) = 1/N . By considering that the time steps last for a random time that is exponentially distributed with mean 1/N , we obtain the following Markov process with transitions given by (for k 2 {0, . . . ,d 1}) The rst line corresponds to a successful transmission, the second to a collision. The last one is because any attempted transmission from state d 1 goes back to state 0 (successful or not). Note that this model has also some transitions that aect k 2 objects but their rate is O (N 1 k ) and they are therefore negligible. Note that, as N goes to innity, we have

which shows that model falls into our density-dependent population process introduced in Section 2.1.