Luca Bortolussi

Nicolas Gast

Mean-Field Limits Beyond Ordinary Differential Equations

Keywords: Population models, Markov chain, Mean-field limits, Differential inclusions, Hybrid systems

We study the limiting behaviour of stochastic models of populations of interacting agents, as the number of agents goes to infinity. Classical mean-field results have established that this limiting behaviour is described by an ordinary differential equation (ODE) under two conditions: (1) that the dynamics is smooth; and (2) that the population is composed of a finite number of homogeneous sub-populations, each containing a large number of agents. This paper reviews recent work showing what happens if these conditions do not hold. In these cases, it is still possible to exhibit a limiting regime at the price of replacing the ODE by a more complex dynamical system. In the case of non-smooth or uncertain dynamics, the limiting regime is given by a differential inclusion. In the case of multiple population scales, the ODE is replaced by a stochastic hybrid automaton.

Introduction

Many systems can be effectively described by stochastic population models, for instance biological systems [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF], epidemic spreading [START_REF] Andersson | Stochastic Epidemic Models and Their Statistical Analysis[END_REF], queuing networks [START_REF] Boudec | Performance Evaluation of Computer and Communication Systems[END_REF]. These systems are composed of a set of objects, agents, or entities interacting together. Each individual agent is typically described in a simple way, as finite state machines with few states. An agent changes state spontaneously or by interacting with other agents in the system. All transitions happen probabilistically and take a random time to be completed. By choosing exponentially distributed times, the resulting stochastic process is a continuous-time Markov chain with a finite state space. Many numerical techniques exist to compute probabilities of such chains [START_REF] Baier | Model-checking algorithms for continuous-time Markov chains[END_REF], and they are part of state-of-the art stochastic model checking tools like PRISM [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF] or MRMC [START_REF] Katoen | A Markov reward model checker[END_REF].

These techniques, however, are limited in their applicability, as they suffer from the state-space explosion: the state-space grows exponentially with the number of agents and even simple agents, when present in large quantities, can generate a huge state space which is far beyond the capabilities of current tools.

This results in the need for approximation techniques to estimate the probabilities and the behaviours of the system. A classic way is to resort to stochastic simulation, which scales better but is still a computationally intensive process for large populations. Precisely in this regime of large populations, mean field analysis offers a viable, often accurate, and much more efficient alternative. The basic idea of mean field is that, when counting the number of agents that are in a given state, the fluctuations due to stochasticity become negligible as the number of agents N grows. For large N , the system becomes essentially deterministic.

A series of results, e.g., [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Bortolussi | Continuous Approximation of Collective Systems Behaviour: a Tutorial[END_REF][START_REF] Kurtz | Solutions of Ordinary Differential Equations as Limits of Pure Jump Markov Processes[END_REF], have established that when the state space of each agent is finite and the dynamics is sufficiently smooth, the system's behaviour converges as N goes to infinity to a limiting behaviour that is described by system of ordinary differential equations (ODE). The dimension of this system of ODE is equal to the number of states of the individual agents, but independent of the population size N . The dimension of the differential equation is typically small, hence the numerical integration of these equations is extremely fast. These results show that the intensity of the fluctuations goes to zero as 1/ √ N . This approach is used in many domains, including computer-based systems [START_REF] Gast | Mean field limit of non-smooth systems and differential inclusions[END_REF][START_REF] Gast | Transient and Steady-state Regime of a Family of List-based Cache Replacement Algorithms[END_REF][START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF][START_REF] Van Houdt | A Mean Field Model for a Class of Garbage Collection Algorithms in Flash-based Solid State Drives[END_REF], epidemic or rumour propagation [START_REF] Chaintreau | The age of gossip: spatial mean field regime[END_REF][START_REF] Hu | Optimal channel choice for collaborative ad-hoc dissemination[END_REF] or bike-sharing systems [START_REF] Fricker | Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity[END_REF]. Is it also used to construct approximate solutions of stochastic model checking problems [START_REF] Bortolussi | Model checking single agent behaviours by fluid approximation[END_REF][START_REF] Bortolussi | Model Checking Markov Population Models by Central Limit Approximation[END_REF][START_REF] Bortolussi | Stochastic approximation of global reachability probabilities of Markov population models[END_REF][START_REF] Bortolussi | Fluid Model Checking of Timed Properties[END_REF].

However, these limiting results have two main shortcomings. First, these models cannot deal with discontinuities on the rates of interaction between agents, or uncertainty in model parameters in an obvious way. Second, being able to approximate the number of agents by a continuous variable requires all populations to be large. These limitations are essentially due to restricting the attention to a limiting regime that can be expressed in terms of smooth ODEs.

In this document, we show that by enlarging the set of possible limiting regimes, it is possible to extend the classical framework in multiple directions. We first begin in Section 2 by a concise introduction to classical mean field models and their ODE limits. This section requires basic knowledge of CTMCs and ODEs. We then show in Section 3 how discontinuities and uncertainties can be treated uniformly and consistently considering mean field limits in terms of differential inclusions. We then tackle the presence of multiple population scales in Section 4. We show that when the number of agents in some populations go to infinity while others remain finite, the mean field limit is naturally expressed as a stochastic hybrid automaton, where continuous-deterministic and discretestochastic dynamics coexist and modulate each other. Last, we mention other related work and extensions of this framework, for instance to cooperative games, in Section 5.

The classical mean field framework

In this section, we will introduce the fundamental mean field approximation. We assume the reader familiar with basic concepts of Markov Chains in Continu-ous Time, see e.g. [START_REF] Norris | Markov chains[END_REF][START_REF] Durrett | Essentials of Stochastic Processes[END_REF] for an introduction, and with ordinary differential equations.

We start in Section 2.1 by introducing a framework to describe the class of systems amenable of mean field analysis, namely Markov population processes (see also the chapter on spatial representations [START_REF] Galpin | Spatial representations and analysis techniques[END_REF]). We illustrate these concepts in Section 2.2 by means of a classic epidemic spreading model. In Section 2.3, we describe the basic mean-field theorems.

Population Continuous-Time Markov Chains

Population continuous-time Markov chains (PCTMCs) describe a set of interacting agents, which can have different internal states. Interactions involve a small number of agents, and can happen randomly in time, according to an exponential distribution with system-dependent rate. We describe these systems in terms of counting variables and transition classes, following the conventions of [START_REF] Bortolussi | Continuous Approximation of Collective Systems Behaviour: a Tutorial[END_REF][START_REF] Henzinger | Formalisms for Specifying Markovian Population Models[END_REF].

More specifically, a PCTMC M model is a tuple (X, T , x 0 , N ), where -X = (X s ) s∈S ∈ R |S| is the population vector. The state space of an agent is S and X s ∈ N counts the number of agents in state s ∈ S. The state space of the model is a finite or countable subset of R |S| . -T is the set of transition classes, each of the form η = (φ η (X), v η , f η (X)), where • φ η (X) ∈ {0, 1} is a guard predicate, representing a subset of S in which the transition is active; • v ∈ R n is an update vector, encoding the relative change of X induced by the firing of the transition η: the new state will be X + v η ; • f η (X) is the rate function, giving the rate at which an η transition is fired as a function of the state space of the system. Typically, f η (X) is a (locally) Lipschitz continuous function of the population variables. -x 0 ∈ R |S| is the initial state of the system. -N is the population size.

Each PCTMC model M defines a CTMC X(t) on the state space S. This chain is characterised by the inifinitesimal generator matrix Q [START_REF] Durrett | Essentials of Stochastic Processes[END_REF], whose offdiagonal entries are given by

Q x1,x2 = η∈T s.t. x2=x1+vη φ η (x 1 )f η (x 1 ).
An important concept related to population models is the system size, N . Typically, system size is the total (initial) population. In some domains, though, like biochemical networks or ecological models, N may represent another measure of size, like the volume or the area in which the dynamics described by a PCTMC happens. We refer to [START_REF] Bortolussi | Continuous Approximation of Collective Systems Behaviour: a Tutorial[END_REF] for a deeper discussion of this.

Example: SIR epidemic spreading

As a simple and illustrative example, we consider the spreading of a disease in a community of N agents (which can be humans, animals, computers). The state space of an agent is S = {S, I, R}. This model is one of the classical examples of a Markov population process and is referred to as the SIR model.

The contagion happens when a susceptible agent (X S ) enters in contact with another agent who turns out to be infected (at rate k si X S X I /N ) or enters in contact with an external source of the disease (at rate k i S). Infected individuals spontaneously recover at rate k r , and become Recovered (X R ) and immune from the disease. This immunity, however, can be lost with rate k s .

Formally, the model can be described as a PCTMC with three variables (X S , X I , X R ), each taking values in the integers {0, . . . , N }, as no birth or death events are considered. The model has four transition classes, all having a guard predicate evaluating to true (= 1) in all states :

-Internal infection: (true, e I -e S , k si X S X I /N ); -External infection: (true, e I -e S , k i X S ); -Recovery: (true, e R -e I , k r X I ); -Immunity loss: (true, e S -e R , k s X R );

Classic mean field equations

mean field theory answers the following question about population models: what happens when the population is very large? More specifically, it can be shown that, for a large class of models, the dynamics of the system greatly simplifies as the system size goes to infinity. The classic theorem, dating back to the 1970s [START_REF] Kurtz | Solutions of Ordinary Differential Equations as Limits of Pure Jump Markov Processes[END_REF], shows that trajectories of suitably rescaled processes for large populations look deterministic, and in fact converge to the solution of an ordinary differential equation (ODE).

An important operation in the path to mean field is to normalise population processes, dividing variables by the system size, and updating accordingly the transitions. This allows one to compare different models, as they will now have the same scale, intuitively these are population densities. Roughly, the deterministic behaviour appears because fluctuations around the mean of a population process grow as √ N , hence while normalising, i.e. dividing by N variables, fluctuations will be of magnitude 1/ √ N , and will thus go to zero. More formally, consider a population model M N , where we make explicit the dependence on the system size, and define its normalised version MN as follows:

-Population variables (and initial conditions) are rescaled by N : XN = X/N ; -Transition rates and guard predicates are expressed in the normalised variables, by substituting N XN for X in the functions: f N ( XN ) = f (N X) and φN ( XN ) = φ(N X); -Update vectors are rescaled by N , too: vN = v/N ;

The CTMC associated with the normalised model will be denoted by XN (t).

Example. Consider the SIR model of Section 2.2. Its normalised version, for a population of N agents, has the following four transition classes:

-Internal infection: vsi N = (e I -e S )/N , f N si ( X) = N k si XS XI ; -External infection: vi N = (e I -e S )/N , f N i ( X) = N k i XS ; -Recovery: vr N = (e R -e I )/N , f N r ( X) = N k r XI ; -Immunity loss: vs N = (e S -e R )/N , f N s ( X) = N k s XR . As we can see, all transition rates depend linearly on system size. When this happens, rates are called density dependent [START_REF] Bortolussi | Continuous Approximation of Collective Systems Behaviour: a Tutorial[END_REF], a condition that usually guarantees the applicability of the mean field results.

Drift. The main quantity required to define mean field equations is the drift. The drift is the average direction of change of the population model, conditional on being in a certain state at some time t. The drift of the normalised model is

F (x) = η∈T φη (x) fη (x)v η , (1)
Usually, mean field is defined under some additional restrictions on the population model:

(C1) Guards are true for any x ∈ S, hence the indicator function can be safely removed from the drift:

F N (x) = η∈T f N η (x)v N η . (C2) F is a Lipschitz continuous function.
Note that by definition of the rescaled model, the drift F (x) does not depend on N , because the update vectors are rescaled by 1/N while the transition rates are rescaled by N . When the drift F N (x) does depend on N , condition (C2) can be replaced in all theorems by a condition (C2'): F N (x) converges uniformly as N → ∞ to a Lipschitz continuous function F (x).

Note that in the SIR model, Conditions (C1) and (C2) are satisfied. For the second one, in particular, we can see that by multiplying a normalised update vector, e.g. vi N , by the corresponding rate, e.g. f N i ( X) = N k i XS XI , the dependency on system size cancels out, so that the drift F N (x) is independent of N . Lipschitz continuity3 is easily proved. See [START_REF] Bortolussi | Continuous Approximation of Collective Systems Behaviour: a Tutorial[END_REF] for a deeper discussion of these conditions. In the following, we will discuss how to weaken these assumptions.

The following theorem can then be proved (see [START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF][START_REF] Bortolussi | Continuous Approximation of Collective Systems Behaviour: a Tutorial[END_REF][START_REF] Darling | Differential equation approximations for Markov chains[END_REF]):

Theorem 1. Assume the above conditions C1 and C2 hold and that X0 converges to x 0 almost surely (resp. in probability) as N goes to infinity. Let x be the solution of the ODE:

d dt x(t) = F (x(t)) x(0) = x 0 . (2)
Then, for any T > 0, lim

N →∞ sup t≤T XN (t) -x(t) = 0 almost surely (resp. in probability).
The theorem essentially states that trajectories of the PCTMC, for large N , will be indistinguishable from the solution of the mean field ODE restricting to any finite time horizon T > 0. This can be seen as a functional version of the law of large numbers. An example of the theorem at work, for the SIR model, can be seen in Figure 1. For the SIR model, the mean field approximation is given by the following system of ODEs:

ẋS = -k i x S -k si x S x I + k s x R ẋI = k i x S + k si x S x I -k r x I (3) ẋR = k r x I -k s x R
We simulate the model for population size of N = 10, N = 100 and N = 1000 agents and report the evolution of the numbers of susceptible or infected agents as a function of time. The parameters are

k i = k si = k s = k r = 1 and the initial conditions are X I (0) = X S (0) = 2X R (0) = 2N/5
. Each plot contains three curves: a sample path of one simulation, the mean field (ODE) approximation and an average over 104 simulations. This figure illustrates how large the population size N has to be for the approximation to be accurate. We observe that the mean field approximation describes correctly the overall dynamics of the PCTMC for N = 100 and N = 1000. In fact, it can be shown that the rate of convergence in Theorem 1 is of the order of 1/ √ N [START_REF] Darling | Differential equation approximations for Markov chains[END_REF][START_REF] Benaim | A class of mean field interaction models for computer and communication systems[END_REF] but with bounds that often an underestimation of the real convergence speed. In practice, we often observe that the convergence is quicker than this bound. This is particularly true when one considers the average stochastic value: E [X(t)]. To illustrate this fact, we simulated the SIR model 10 4 times to compute the values E [X(t)]. We report the evolution of E [X S (t)] and E [X I (t)] with time in Figure 1. We observe than, already for N = 10, the ODE x(t) is extremely close to the value of E [x(t)] computed by simulation. For N = 100 and N = 1000, the curves are indistinguishable.

Notice that in Theorem 1, the restriction to finite time horizons is fundamental, as convergence at steady state does not necessarily hold. An example is given by the SIR model when

k i = 0, k si = 3, k s = k r = 1 and x R (0) = x S (0) = x I (0) = 1/3.
This initial condition is a fixed point of the mean field ODE which therefore predicts an endemic equilibrium. In the PCTMC model, however, for any N , the epidemic always extinguishes (i.e. eventually X I (t) = 0) because there is no external infection in this situation. For a deeper discussion of this issue, see the next section as well as [START_REF] Bortolussi | Continuous Approximation of Collective Systems Behaviour: a Tutorial[END_REF].

In fact, steady state results can be obtained at the price of adding two additional conditions:

(C3) For any N , the PCTMC has a steady-state distribution π N . The sequence of distributions π N is tight. 4 (C4) The ODE (2) has a unique fixed point x * to which all trajectories converge. [XI] (ave. over 10000 simu) xI (mean-field approximation) Fig. 1. Simulation of the SIR model: comparison between ODE and simulation for various values of N . We observe that the simulation converges to the ODE as N goes to infinity. Moreover, even for N = 10, the average simulation is very close to the ODE.

The condition C3 is a natural condition that is in general true for PCTMC models. For example, if the population CTMC of size N has a finite number of states, it has a steady-state distribution. Moreover, the tightness of the measure is true if X is almost surely bounded.

However, condition C4 is a condition that is often difficult to check for a given set of ODEs. Proving it requires to exhibit a Lyapunov function witnessing global attractiveness of the unique equilibrium point.

Theorem 2. Assume that the above conditions C1, C2, C3 and C4 hold. Then, π N converges weakly to the Dirac measure x * .

Similarly to Theorem 1, it can be shown that under mild additional conditions, the speed of convergence of the steady-state distribution π N to x * is also 1/ √ N . For example, it is shown in [START_REF] Ying | On the Rate of Convergence of Mean-Field Models: Stein's Method Meets the Perturbation Theory[END_REF] that this holds when x * is exponentially stable, in which case we have

E X N (∞) -x * 2 = O(1/ √ N ),
where X N (∞) denotes a random point distributed according to the stationary measure π N .

Non-continuous dynamics and uncertainties

The classical mean field models make the assumption that the drift F , given by ( 1), is Lipschitz-continuous. Yet, this is not the case in many practical problems. For example, this occurs when a transition η has a guard predicate φ η that is true only in a sub-part of the domain. This causes a discontinuity in the drift between the two regions where the predicate is true or false. In this case, the ODE ẋ = F (x) is often not well-defined. A classical way of overcoming this difficulty is to replace the ODE by a differential inclusion ẋ ∈ F (x).

In this section, we give some general discussion of differential inclusions. We then show how the classical mean field results can be generalised to differential inclusion dynamics in the case of discontinuous dynamics. Last, we show in Section 3.3 how this framework can be used to deal with uncertainties in the parameters. This section collects results from [START_REF] Bortolussi | Mean Field Approximation of Imprecise Population Processes[END_REF][START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF].

The differential inclusion limit

Let M = (X, T , x 0 , N ) be a population model as defined in Section 2.1. Each transition class has the form η = (φ(X), v η , f η (X)). Recall that the drift of the stochastic system, defined in Equation (1), is

F (x) = η∈T φη (x) fη (x)v η , (4)
The classical mean field results presented in the previous section (Theorems 1 and 2) apply when the guard predicates φ η always evaluate to true. This condition guarantees that, if the functions f η are Lipschitz continuous, the function F is also Lipschitz continuous. This ensures that the ODE ẋ = F (x) is welldefined: from any initial condition, it has a unique solution. However, this no longer holds when guard predicates can take the two values true and false. In this case, the drift F is not continuous and Theorems 1 and 2 no longer apply. One of the reasons for the inapplicability of those theorems is that when F is not a Lipschitz-continuous function, the ODE ẋ = F (x) does not necessarily have a solution. For example, let F : R → R be defined by

F (x) = -1 if x ≥ 0 and F (x) = 1 if x < 0. The ODE ẋ = F (x)
starting in 0 has no solution. A natural way to overcome this limitation is to use differential inclusions.

Let G be a multivalued map, that associates to each x ∈ R S a set G(x) ⊂ R S . A trajectory x is said to be a solution of the differential inclusion ẋ ∈ G(x) starting in x 0 if:

x(t) = x 0 + t 0 g(s)
ds, where for all s: g(s) ∈ G(x(s)).

The sufficient condition for the existence of at least one solution x : [0, ∞) → R d of a differential inclusion from any initial condition x 0 is the following (see [START_REF] Aubin | Differential Inclusions[END_REF]):

(C5) For all x ∈ R d , G(x) is closed, convex and non-empty; sup x∈R d |G(x)| < ∞
and G is upper-hemicontinuous. 5Theorem 3 ( [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF]). Let M = (X N , T , x 0 , N ) be a population model with drift F such that there exists a function G satisfying (C5) such that for all x ∈ R d : F (x) ∈ G(x). Let S be the set of solutions of the differential inclusion ẋ ∈ G(x) starting in x 0 . Then, for all T , almost surely

lim N →∞ inf x∈Sx 0 sup t∈[0,T ] X N (t) -x(t) = 0.
In other words, as N grows, the distance between the stochastic process X and the set of solutions of the differential inclusion goes to 0. If this set has a unique solution, then X N converges to this solution. Theorem 3 is a generalisation of Theorem 1 that relaxes the condition (C1) by using a larger drift function G. We will see in the next section a natural way to choose G when the drift is not continuous. The price to be paid by this generalisation is composed of two drawbacks. First, differential inclusions can have multiple solutions. Theorem 3 implies that X N gets closer to the solutions of the differential inclusion but does not indicate towards which solutions the process will converge. Second, the speed of convergence of X N to S is unknown, apart in the special case of one-sided Lipschitz drift, for which the distance between X N and S decays in 1/ √ N (see [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF]). In general, the convergence appears to be slower in the case of non-continuous dynamics (see Figure 2).

The steady-state of a non-continuous PCTMC can also be approximated by using the same approximation. In fact, the results of Theorem 2 can also be directly generalised to the case of non-continuous dynamics: if the differential inclusion has a unique point x * to which all trajectories converge, then the steady-state distribution of X N concentrates on x * as N goes to infinity (see [START_REF] Bortolussi | Mean Field Approximation of Imprecise Population Processes[END_REF][START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF]).

Application to discontinuous dynamics

A natural way to define a multivalued map that satisfies (C5) is to consider the multivalued map F , defined by

F (x) = ε>0 convex hull   x : x-x ≤ε F (x )   . (5)
It is shown in [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF] that if F is bounded, then F satisfies (C5). By Theorem 3, this implies that, regardless of the properties of the original drift F , the trajectories of the stochastic system X N converge to the solution of the differential inclusion ẋ ∈ F (x).

When F is continuous at a point x ∈ R d , the F (x) = F (x). When F is not continuous in x ∈ R d , F (x) is multivalued. To give a concrete example, let us consider the SIR model of Section (2.2) in which we add an additional transition corresponding to the treatment of some infected people. This treatment is applied when the proportion of infected people is greater than 0.3 and changes an infected person into a susceptible individual at rate k t . This corresponds to a transition class:

-Treatment: η = (1 X I ≥0.3 , e S -e I , k t X I )
Adding this transition to the original ODE (3), the drift is given by

F (x) =   -k i x S -k si x S x I + k s x R + k t 1 x I ≥0.3 x I k i x S + k si x S x I -k r x I -k t 1 x I ≥0.3 x I k r x I -k s x R   , (6)
where the guard predicate leads to the term 1 x I ≥0.3 . This drift is not continuous in x. In fact, it can be shown that because of this discontinuity, the corresponding ODE has no solution on [0, ∞) starting from x 0 = (.4, .4, .2). The corresponding F defined by Equation ( 5) is then given by

F (x) =   -k i x S -k si x S x I + k s x R + k t 1 x I >0.3 x I + k t [0, x I ]1 x I =0.3 k i x S + k si x S x I -k r x I -k t 1 x I >0.3 x I -k t 1 x I =0.3 [0, x I ] k r x I -k s x R   ,
where the notation a + [b, c] denotes the set [a + b, a + c].

It can be shown that the differential inclusion ẋ ∈ F (x) has a unique solution, x. Hence, Theorem 3 applies to show that X N converges to x as N goes to infinity. To illustrate this fact, we simulated the modified SIR model with the treatment policy and report the results in Figure 2. We observe that, as stated by Theorem 3, X N converges to x as N goes to infinity. In this case, the convergence appears to be slower than for the Lipschitz-continuous case. This is especially visible when looking at the average values E [X S ] and E [X I ]: in Figure 1, we observe that for the Lipschitz-continuous case, E [X S ] is almost equal to x s already for N = 10. In the non-continuous case, Figure 2 indicates that E [X I ] does converge to x I but at a much slower rate. [XI] (ave. over 10000 simu) xI (mean-field approximation) Fig. 2. Simulation of the SIR model: comparison between differential inclusions and simulation for various values of N . We observe that the simulation converges to the solution of the differential inclusion as N goes to infinity. Moreover, even for N = 10, the average simulation is very close to the solution of the differential inclusion.

Imprecise and uncertain models

Stochastic models are one way of representing uncertainties in a system but they depend on parameters whose precise values are not always known. The differential inclusion framework is also well adapted to study models with imprecise or unknown parameters. Following [START_REF] Bortolussi | Mean Field Approximation of Imprecise Population Processes[END_REF], we distinguish two ways to model uncertainties in models of complex systems:

-Imprecise scenario: Some parameters ϑ can depend on features of the environment external to the model. We 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. that ϑ t ∈ Θ. -Uncertain scenario: In a simpler scenario, a parameter ϑ is assumed fixed, but its precise value not known precisely. In this case, we just assume that ϑ ∈ Θ, where Θ is the possible set of values of ϑ, as above.

An imprecise or uncertain PCTMC model is a tuple (X, T , x 0 , Θ, N ), where Θ is a set of parameters. The difference with classical PCTMC is that the rate function f η (X, θ) of each transition class η ∈ T depends on a parameter θ ∈ Θ.

The differential inclusion framework can be used to study the limits of imprecise and uncertain PCTMC. For the uncertain scenario, there is a differential inclusion ẋ ∈ F (x, ϑ) associated with each parameter ϑ. Denoting by S ϑ the set of solutions of this differential inclusion, Theorem 3 shows that, as N goes to infinity, any sequence of uncertain trajectories X N converges to the set S uncertain = ϑ∈Θ S ϑ . The differential inclusion corresponding to the imprecise scenario is ẋ ∈ ϑ∈Θ F (x, ϑ). Denoting by S imprecise the set of solutions of this differential inclusion, Theorem 3 shows that, as N goes to infinity, any imprecise trajectory X N converges to S imprecise . Some numerical methods are developed in [START_REF] Bortolussi | Mean Field Approximation of Imprecise Population Processes[END_REF] to compute or approximate the set of solutions of differential inclusions corresponding to the imprecise and the uncertain model. In particular, we obtain the most precise results by describing the set of reachable values of x at time t as a maximisation problem. Then Pontryagin's maximum principle [START_REF] Todorov | Optimal control theory[END_REF] can be used to numerically compute the solution.

Hybrid mean field limits

In the previous sections, we considered scenarios where all populations of the model are large and grow with the system size. This allows one to prove that their density has vanishing fluctuations around the mean, given by the solution of the mean field equation. However, in many practical cases, there may be multiple population scales in a model, typically in the form of some entities being present in small numbers, independent of the total population size. Examples can be found in genetic regulatory networks, where genes are present in a fixed quantity, typically one or few copies, or more generally in the presence of a centralised form of control [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF]. This suggest that in these scenarios we need to consider mean field models in which the continuous and deterministic limit dynamics of parts of the system coexists with the discrete and stochastic dynamics of other parts. Mathematically, this behaviour is captured by stochastic hybrid systems (SHS) [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF][START_REF] Davis | Markov Models and Optimization[END_REF], which will be introduced in the next subsection.

Stochastic Hybrid Systems

We introduce a model of SHS, essentially borrowing from the treatment of [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF][START_REF] Bortolussi | Hybrid) automata and (stochastic) programs: The hybrid automata lattice of a stochastic program[END_REF] of a class of stochastic hybrid processes known as Piecewise-Deterministic Markov Processes (PDMP) [START_REF] Davis | Markov Models and Optimization[END_REF].

We consider two sets of variables, the discrete variables Z = Z 1 , . . . , Z k and the continuous variables Y = Y 1 , . . . , Y m . The former describes populations that remain discrete also in the mean-field limit, while the second describes populations that will be approximated as continuous. We call E = E d × E c the hybrid state space of the SHS, with E d ⊂ N k a countable set of possible values for Z, and with E c ⊂ R m the continuous state space in which variables Y can take values. Each possible value that the vector Z can take is called a discrete mode, and can be identified with a node in a graph describing the transitions of the discrete states. This graph-based point of view is taken in the definition of stochastic hybrid automata, see e.g. [START_REF] Bortolussi | Hybrid) automata and (stochastic) programs: The hybrid automata lattice of a stochastic program[END_REF].

The evolution of the continuous state is governed by an m-dimensional vector field F (Z, Y), depending on the continuous and the discrete variables. Hence, the continuous variables will evolve following the solution of the differential equation defined by F , which can be different in each discrete mode z. Such a mode-specific continuous dynamics is one of the characteristic features of SHS.

The dynamics of the discrete state is governed by a stochastic Markovian dynamics, specified by two quantities: a rate function λ(Z, Y), depending both on discrete and continuous variables, and a jump or reset kernel R(Z, Y, •), specifying for each Z, Y a distribution on E, giving the state in which the system will find itself after a jump of the discrete transition. For the purpose of this chapter, we can restrict ourselves to finitely supported reset kernels, defined by a finite set of pairs of update vectors {(v d j , v c j ) | j = 1, . . . , h} and associated probability functions p j (Z, Y), giving the likelihood of jumping from state Z, Y to state Z + v d j , Y + v c j , if a stochastic event fires when the system is in state Z, Y. This, in turn, happens after an exponentially distributed delay with rate λ(Z, Y).

Discrete and continuous dynamics in a SHS are intertwined. The system starts in a given state z 0 , y 0 , and its continuous state evolves following the solution of the initial value problem d dt Y(t) = F (z 0 , Y(t)), Y(0) = y 0 . This continuous flow will go on until a discrete event will happen, at a random time governed by an exponential distribution with rate λ(Z, Y(t)). Note that, as Y will change value in time following the flow of the vector field, the rate of a discrete jump is also time-dependent. When a discrete transition happens, say in state z, y, then the system will jump to the state z + v d j , y + v c j with probability p j (Z, Y). Note that both the discrete mode and the value of continuous variables can change. From this new state, the system continues to evolve following the dynamics given by the vector field in the new discrete mode. The overall dynamics is given by an alternation of periods of continuous evolution interleaved by discrete jumps. For a proper mathematical formalisation of this process, we refer the interested reader to [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF][START_REF] Davis | Markov Models and Optimization[END_REF].

From PCTMC to SHS

In this section we will show how to construct a SHS from a PCTMC (X, T , x 0 , N ), and how to guarantee the asymptotic correctness of the method. The starting point is a partition of the variables X of the PCTMC into two distinct classes: discrete and continuous. We will denote discrete variables with Z and continuous ones with Y, so that X = Z, Y. Transitions T also have to be separated in two classes: discrete T d and continuous T c . Intuitively, continuous transitions and variables will define the continuous dynamics, and discrete transitions and variables the discrete one. The only request is that continuous transitions do not affect discrete variables, i.e. for each η ∈ T c , v η [Z] = 0, where v η [Z] denotes the vector v η restricted to the components of Z.

Remark 1. The choice of how to partition variables and transitions into discrete and continuous is not obvious, and depends on the model under consideration. Often, this is easily deduced from model structure, e.g. due to the presence of conservation laws with a small number of conserved agents. A further help comes from the request to make explicit in the rates and updates the dependency on system size N . An alternative is to define rules to automatically switch between a discrete or a continuous representation of variables and transitions, depending on the current state of the model. We refer the interested reader to the discussion in [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF] for further details.

Normalisation of continuous variables. To properly formalise hybrid mean-field limits, we need to perform a normalisation operation on the continuous variables, taking system size into account (hence we will use the superscript N from now on). This can be obtained as in Section 2.3, by introducing the normalised variables ŶN = Y/N , and expressing rates, guards, and update vectors with respect to these normalised variables. Note that normalised update vectors vN η are divided by N only in the continuous components, as the discrete variables are not rescaled. Transitions, after normalisation, have to satisfy some scaling conditions:

(Tc) Continuous transitions η ∈ T c are such that f N η (Z, Ŷ)/N → f η (Z, Ŷ), as N → ∞, uniformly on Ŷ for each Z. The limit function is required to be (locally) Lipschitz continuous. Furthermore, their non-normalised update is independent of N . Guards can depend only on discrete variables: φ η = φ η (Z). (Td) Discrete transitions η ∈ T d are such that their rate function f N η (Z, Ŷ) converges (uniformly in Ŷ for each Z) to a continuous function f η (Z, Ŷ). Their normalised jump vector vN η has also to converge to a vector vη as N diverges. Guards can depend only on discrete variables: φ η = φ η (Z).

Note that, for discrete transitions, we consider the change in the normalised continuous variables, and we admit that the update vectors can depend on N . In particular, the update for continuous variables can be linear in N , thus resulting in a non-vanishing jump in the density, in the large N limit. This means that limit discrete transitions may also induce jumps on continuous variables.

Construction of the limit SHS. Given a family of PCTMC models (X, T , x 0 , N ), indexed by N , with a partition of variables into Z, Y and transitions into T c , T d we can normalise continuous variables and formally define the SHS associated with it:

-The vector field defining the continuous dynamics of the SHS is given by the following drift:

F (Z, Ŷ) = η∈Tc v η I{φ η (Z)}f η (Z, Ŷ)
-The jump rate of the SHS is given by

λ(Z, Ŷ) = η∈T d I{φ η (Z)}f η (Z, Ŷ)
-The reset kernel is specified by the pair of update vectors (v η [Z], vη [ Ŷ]) and by the probability

p η (Z, Ŷ) = I{φ η (Z)}f η (Z, Ŷ) λ(Z, Ŷ) ,
for each transition η ∈ T d . -The initial state is z 0 , ŷ0 .

SIR model with vaccination. We consider now an extension of the SIR model of Section 2.2, with the possibility of starting a vaccination campaign of susceptible individuals. The model has an additional variable, X V ∈ {0, 1}, which is going to be the only discrete variable of the system and encodes if the vaccination is in force or not. We further have two additional transitions:

-Vaccination of susceptible: (true, e R -e S , k v X S X V ); -Activation of the vaccination policy: (X V = 0, e V , k a X I /N );

The first transition, which will be a continuous transition, models the effect of vaccination, moving agents from susceptible to recovered state. Note that the rate depends on X V , hence the transition is in force only if X V = 1. We could alternatively specify the same behaviour by introducing a guard in the transition, depending only on discrete variables, which would result in a rate active only in a subset of discrete modes. The second transition will be kept discrete and model the activation of the vaccination policy. Its rate depends on the density of infected individuals (the higher the infected, the higher the activation rate). The guard on X V allows the activation transition to be in force only when the vaccination is inactive.

In Figure 3 (left), we show a trajectory of the system for N = 100, and compare it with a trajectory of the limit SHS. Parameters are

k si = 0, k i = 1, k r = 0.1, k s = 0.01, k v = 2, k a = 2, X S (0) = 0.95N , X I (0) = 0.05N , X R (0) = 0.
As we can see, around time t = 2.5 there is a sudden drop on the number of susceptible individuals, caused by the beginning of vaccination. The similarity between the PCTMC and the SHS trajectories is a clear hint on the existence of an underlying convergence result.

Mean-field convergence results. Consider a family of PCTMC model (X, T , x 0 , N ), and denote by (Z N (t), ŶN (t)) the normalised CTMC associated with it, for size N , where we made explicit the partition of variables into discrete and continuous. Denote by (z(t), y(t)) the limit SHS, constructed according to the recipe of this section. We can then prove [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF] the following theorem: Theorem 4. Assume (z N 0 , ŷN 0 ) → (z 0 , y 0 ) and that the transitions of the PCTMC model satisfy conditions (Tc) and (Td). Then (Z N (t), ŶN (t)) ⇒ (z(t), y(t)), for all times t ≥ 0, where ⇒ denotes weak convergence. 6This theorem states that the distribution of the PCTMC will look like the distribution of the SHS for large N . In particular, if E d is finite (i.e. there is a finite number of discrete modes) and E c is compact, then all conditional and unconditional moments of the distribution converge.

In Figure 3 (right), we see the theorem at work in the epidemic with vaccination example. The figure compares the empirical cumulative distribution of the density of infected individuals at time t = 10. The curves look quite similar already for N = 100, and are almost identical for N = 1000.

Extensions of the hybrid limit framework

The hybrid mean field limit presented in the previous section can be extended in many ways, as discussed in [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF]. Here we will sketch them briefly, referring the interested reader to [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF] for further details. One possible direction to enrich the framework is to consider forced transitions. In the context of SHS, these are discrete jumps happening as soon as a condition on the system variables becomes true. Typically, they are introduced by constraining the continuous state space E c (in a mode dependent way), and forcing a jump to happen as soon as the trajectory of the continuous variables hits the boundary ∂E c [START_REF] Davis | Markov Models and Optimization[END_REF]. Then a jump is done according to the reset kernel R, whose definition has to be extended on the boundary ∂E c . Hence, discrete jumps may happen at stochastic times, or when the condition for a forced jump is met.

In the PCTMC setting, introducing forced transitions requires us to allow transitions with an infinite rate and with a non-trivial guard, firing as soon as their guard becomes true. Their guards, then, can be used to constrain the continuous state space. Hence, E c will be defined in each mode as the interior of the complement of the region obtained by taking the union of the satisfaction sets of all the guards of forced transitions. The reset kernel in a point of the boundary ∂E c will then be defined by the active immediate transitions at that point, choosing uniformly among the active transitions. 7As an example, consider again the SIR model with vaccination, but assume its activation and deactivation is threshold-based: when the density of infected becomes greater than a threshold I high , the vaccination is started, while if it falls below I low , the vaccination is stopped. In the PCTMC model, we would have the following two transitions in place of the stochastic one discussed previously:

-Activation of the vaccination policy: (X I ≥ N I high , e V , ∞); -Deactivation of the vaccination policy: (X I ≤ N I low , -e V , ∞).

Theorem 4 is readily extended to the of instantaneous transitions, under some additional technical conditions on the vector fields (called transversal crossing), see [START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF] for details. The validity of the result for the SIR model with vaccination is illustrated in Figure 4, where we compare simulations of the PCTMC model and the limit SHS, for thresholds I high = 0.3 and I low = 0.2. Other extensions of the hybrid mean field include dealing with guards in deterministic and stochastic transitions. These introduce discontinuities in the model, and require further technical conditions for the limit theorems to hold. As for guards in the continuous transitions, this in fact requires one to introduce in the hybrid context the mean field techniques based on differential inclusions of Section 3.

Related work and examples

Load balancing and discontinuous dynamics. The use of mean field approximation is popular for studying load-balancing policies in server farms. In such a system, an object represents a server and its state is typically the number of jobs that are waiting in its queue. A popular randomised load balancing policy is the two-choice policy: for each incoming packet two servers are picked at random and the job is allocated to the least-loaded of the two. This policy has been successfully analysed by classical mean field techniques in [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF] where it is shown that it leads to an important gain of performance compared to a purely random allocation. The classical approach then fails when one considers a centralised load balancing policy such as join the shortest queue because these policies lead to discontinuous dynamics. As demonstrated in [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF][START_REF] Tsitsiklis | On the power of (even a little) resource pooling[END_REF][START_REF] Fricker | Incentives and redistribution in homogeneous bike-sharing systems with stations of finite capacity[END_REF], these problems can be modelled and resolved by using differential inclusions.

Heterogeneous systems and uncertainties. Another problem where differential inclusions can help is the case of heterogeneous systems. In such cases, there is a large number of objects each having distinct parameters. One example is the caching model of [START_REF] Gast | Transient and Steady-state Regime of a Family of List-based Cache Replacement Algorithms[END_REF] in which an object i has a popularity p i . One possibility to solve the problem is to consider a system of N ×S ODEs, where N is the number of objects and S the dimension of the state space. This method scales linearly in the number of objects but might still be problematic for large populations. An alternative is to consider upper and lower bounds on the dynamics [START_REF] Tschaikowski | Approximate Reduction of Heterogenous Nonlinear Models with Differential Hulls[END_REF] or to study a PDE approximation when the number of objects is large, see for example [START_REF] Fricker | Mean field analysis for inhomogeneous bike sharing systems[END_REF].

Hybrid mean field limits. The use of hybrid approximation of population models is quite common in the area of systems and synthetic biology, where genes and often mRNA molecules are present in such low numbers (genes usually in one copy) that classic mean field assumptions are not correct and can lead to models failing to capture important features of the system like bursting protein expression [START_REF] Krn | Stochasticity in gene expression: from theories to phenotypes[END_REF]. In addition to the investigation of hybrid limits, carried out independently for general population models [START_REF] Bortolussi | Limit Behavior of the Hybrid Approximation of Stochastic Process Algebras[END_REF][START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF] and more specifically for gene networks [START_REF] Crudu | Convergence of stochastic gene networks to hybrid piecewise deterministic processes[END_REF], considerable work has been done in hybrid simulation [START_REF] Pahle | Biochemical simulations: stochastic, approximate stochastic and hybrid approaches[END_REF] and in developing moment closure techniques for hybrid limits [START_REF] Hasenauer | Method of conditional moments (MCM) for the Chemical Master Equation: A unified framework for the method of moments and hybrid stochastic-deterministic models[END_REF].

Mean field games. Game theory studies the decisions taken by competing rational agents. Recently, the notion of mean field games has been introduced in [START_REF] Lasry | Mean field games[END_REF][START_REF] Huang | Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle[END_REF] to model decisions in systems composed of a large number of agents. In a mean field game, each agent tries to minimise an objective function that depends on the average behaviour of the population of agents but not on the action of a precise agent. This simplifies the analysis of Nash equilibria that are replaced by mean field equilibria. This theory is used for modelling purposes but also to solve optimisation problems in a decentralised way [START_REF] Yang | A mean-field control-oriented approach to particle filtering[END_REF][START_REF] Gast | Impact of demand-response on the efficiency and prices in real-time electricity markets[END_REF]. It can be shown that in certain cases, mean field equilibria are the limit of a sub-class of Nash equilibria [START_REF] Doncel | Mean-Field Games with Explicit Interactions[END_REF].

Modelling languages supporting a mean field semantics. In the last ten years, there has been a considerable interest in extending stochastic modelling languages, in particular stochastic process algebras, in order to derive automatically mean field equations. Examples are e.g. [START_REF] Tribastone | Scalable differential analysis of process algebra models[END_REF][START_REF] Ciocchetta | Bio-PEPA: A framework for the modelling and analysis of biological systems[END_REF][START_REF] Bortolussi | Dynamical Systems and Stochastic Programming: To Ordinary Differential Equations and Back[END_REF]. This work has also been extended to generate hybrid semantics in [START_REF] Bortolussi | Hybrid Dynamics of Stochastic Programs[END_REF][START_REF] Bortolussi | Hybrid) automata and (stochastic) programs: The hybrid automata lattice of a stochastic program[END_REF][START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF]. See also the chapter of this book on CARMA [START_REF] Loreti | Modeling and analysis of collective adaptive systems with CARMA and its tools[END_REF].

Conclusion

In this document, we reviewed the notion of mean field limits of a stochastic population process and presented two extensions. The classical mean field results show that, under some conditions, a stochastic population process converges to a deterministic system of ODEs as the number of objects of the population grows. We have shown that, by replacing the system of ODEs by either a differential inclusion or a hybrid system, it is possible to enlarge the set models for which a mean field limit exists. We illustrated these notions by using a classical SIR example. The last section gives a few pointers to papers in which these frameworks can be applied or generalised.

  simulation)[XS] (ave. over 10000 simu) xs (mean-field approximation) )[XI] (ave. over 10000 simu) xI (mean-field approximation) simulation)[XS] (ave. over 10000 simu) xs (mean-field approximation) )[XI] (ave. over 10000 simu) xI (mean-field approximation) simulation)[XS] (ave. over 10000 simu) xs (mean-field approximation)

  (ave. over 10000 simu) xI (mean-field approximation) (ave. over 10000 simu) xI (mean-field approximation) simulation) [XS] (ave. over 10000 simu) xs (mean-field approximation)

Fig. 3 .

 3 Fig.3. Left: comparison of a simulated trajectory of the PCTMC model, for N = 100, with a simulated of the limit SHS. Right: comparison of the empirical cumulative distribution of XI at time t = 10 of the limit SHS and the PCTMC models for N = 100 and N = 1000.

Fig. 4 .

 4 Fig. 4. Left: comparison of a simulated trajectory of the PCTMC model with instantaneous transitions (normalised variables), for N = 100 (left) and N = 1000 (right), with a simulated trajectory of the limit SHS.

In fact, Lipschitz continuity is satisfied only locally, but this enough for mean field convergence to work.

A sequence of distributions is π N tight if their probability does not escape to infinity, i.e. for each there is a compact set K such that π N (K) ≥ 1 -for each N .

G is upper-hemicontinuous if for all x, y ∈ R d , xn ∈ R d , yn ∈ F (xn), limn→∞ xn = xand limn→∞ yn = y, then y ∈ F (x).

In fact, weak convergence holds for (z, y) as processes in the Skorokhod space of cadlag functions, see[START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF]. For a definition of weak convergence, see[START_REF] Billingsley | Probability and measure[END_REF].

In[START_REF] Bortolussi | Hybrid behaviour of Markov population models[END_REF], weights are introduced to solve non-determinism between instantaneous transitions. Furthermore, the possibility of seeing a chain of instantaneous events firing is taken into account. Termination of this chain is discussed in[START_REF] Bortolussi | Hybrid) automata and (stochastic) programs: The hybrid automata lattice of a stochastic program[END_REF] (where it is proved undecidable for countable state spaces), and in[START_REF] Galpin | HYPE: Hybrid modelling by composition of flows[END_REF], where sufficient and testable conditions for termination are given.


Mathieu Petitcuenot

Bernard Anselmetti

Keywords: Tolerancing, Defect analysis, Digital straightening, Complex surface

 for flexible parts, which enables us to define a stressed state in order to measure the part by straightening it to simulate its position in the mechanism. The originality of this approach is that the straightening operation is performed numerically by calculating the displacement of a cloud of points. The results lead to a quantification of the global defects through various simple models and enable us to extract local defects. The outcome is an acceptable tolerance solution. The procedure is first developed for the simple example of a steel bar with a rectangular cross section, then applied to an industrial case involving a complex 3D surface of a turbine blade. The specification is described through ISO standards both in the free state and in the straightened state.

Fig. 1. Examples of parts and defects

The detection of local defects requires the surface to be measured using a large number of closelyspaced points. The resolution must be adjusted to the size of the defect being sought. For example, using scanning or optical measurement techniques, the resolution can be as small as 0.3 mm. However, the objective is not to measure the roughness or waviness of the part.

Traditional form specification is unsuitable because acceptable global defects can be much larger than the maximum acceptable local defects.

Metrology techniques are used to compare the measured part to the nominal CAD model and display the deviations at various points. The specification is satisfied if all the points are inside the tolerance zone defined by the ISO 1101 standard [START_REF]Geometrical product specifications (GPS) -Geometrical tolerancing -Tolerances of form, orientation, location and run-out[END_REF].

In our new approach, an intermediate, distorted model (e.g. through twisting or bending) is defined. The residual deviation between the distorted model and the measured part is calculated. In addition, the parameters of the deformation can be used to assess the quality of the actual part.

Bibliographical review

This work discusses two applications in the fields of metrology, defect analysis, calculation and tolerancing. The measurements were carried out using a coordinate measuring machine for the first application and an optical measuring device for the turbine blade.

All measuring devices have their own characteristics and particularities. Li [START_REF] Li | Free-form surface inspection techniques state of the art review[END_REF] presented a general state-of-the-art of various types of measuring machines. Even though the emphasis is primarily on CMM measurement, the characteristics, advantages and drawbacks of contact-free techniques for free-form inspection are also discussed.

Similarly, Savio [START_REF] Savio | Metrology of freedom shaped parts[END_REF] presented and compared contact and non-contact inspection tools along with some applications. His study focuses more particularly on the shape of a turbine blade and the problem of the distribution of points along the aerodynamic profile. The difficulty is that the curvature of an airfoil becomes very high near the thin edges. The choice of the density of the points on an airfoil is a common concern. In our work, the points were sampled on the pressure side and the suction side, excluding the edges in order to avoid problems in fitting the actual profile to the CAD model.

Several studies were aimed at the identification of manufacturing defects in parts. Radvar-Esfahlan [START_REF] Radvar-Esfahlan | Nonrigid geometric metrology using generalized numerical inspection fixtures[END_REF][START_REF] Radvar-Esfahlan | Robust generalized numerical inspection fixture for the metrology of compliant mechanical parts[END_REF] proposed an inspection methodology based on a knowledge of predictable defects, such as the effect of gravity or residual strains on large slender parts. Assuming conservation of distance, he applied a finite element simulation to a CAD model, then performed rigid registration with the sampled scanned workpiece, and finally used non-rigid registration to eliminate the springback effect of the manufacturing process.

The registration of point clouds with the CAD model is an important step in the analysis process. Many works have used Bourdet's concept of small displacement torsor [START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF], which lends itself particularly well to small strains and small displacements. The next references used that approach to build mathematical models of defects.

Zhu [START_REF] Zhu | Efficient registration for precision inspection of free-form surfaces[END_REF] conducted an exhaustive study on the efficient registration of free-form surfaces. Along with Savio [START_REF] Savio | Metrology of freedom shaped parts[END_REF], he raised the issue of the point density of the sample. He used the root-mean-square error criterion, which depends on the entire sample of points rather than on the end points or boundary points alone. The least-squares method we used in this paper for registration and calculation led us to choose a homogeneous point distribution. Zhu noted that after he performed the registration the turbine blade became "slightly twisted": our work shows how such a defect can be identified and quantified.

Various types of models for the representation of parts with defects can be found in the literature along with the association method used to fit the model to the measured element.

The modal method developed by Samper [START_REF] Samper | Form Defects Tolerancing by natural Modes Analysis[END_REF] uses the identification of real surfaces from a defect base composed of a series of eigenmodes of the surfaces. The first eigenmodes are associated with a surface's size and position, then with its orientation. The subsequent modes present ripples of decreasing periods which can be used to represent sheer, concavity and all the smaller ripples which can occur on a surface. A Fourier model is used for the decomposition. This work can be extended to roughness.

In [START_REF] Adragna | Analysis of shape deviations of measured geometries with a modal basis[END_REF], Adragna presented a state-of-the-art review of the modal method, beginning with the techniques based on vibration studies, such as Formosa's work [START_REF] Formosa | Modal expression of form defects[END_REF]. In addition, Adragna introduced shape defects called "technological modes". He proposed to create a modal basis composed of defects which may occur during the manufacturing process, but cannot be defined directly by decomposition into vibration modes. One of the drawbacks of the decomposition into eigenmodes is that the mode shapes are continuous, whereas a technological mode can include a peak, a step or an abrupt change of slope in the geometry. Adragna recommended that technological modes be put at the beginning of the modal basis because they are likely to occur frequently.

A similar approach was followed by Favrelières [START_REF] Favreliere | 3D statistical analysis and representation of form error by a modal approach[END_REF], who characterized the actual geometry using a decomposition into an orthogonal basis of vibration modes.

Other methods of the same type have been developed, e.g., in the case of plane parts, the decomposition of shape defects into Fourier series by Capello [START_REF] Capello | Harmonic fitting approach for plane geometry measurements[END_REF] and the identification of shape defects by discrete cosine transform by Huang [START_REF] Huang | Mode-based decomposition of part form error by discrete-cosinetransform with implementation to assembly stamping system with compliant parts[END_REF].

The analysis of shape defects in cylindrical surfaces was addressed in several works. For example, Summerhays [START_REF] Summerhays | Optimizing discrete point sample patterns and measurement date analysis on internal cylindrical surfaces with systematic form deviations[END_REF] and Henke [START_REF] Henke | Methods for evaluation of systematic geometric deviations in machined parts and their relationships to process variables[END_REF] recreated shape errors using an association of Fourier series and Chebyshev polynomials completed by a basis made of measured eigenmodes. Cho [START_REF] Cho | Roundness modeling of machined parts for tolerence analysis[END_REF] carried out a similar study and compared his results with those of Damir [START_REF] Damir | Approximate harmonic models for roundness profiles[END_REF].

These types of analyzes are applicable when the shape or the cause of the defect is unknown. Otherwise, the origin of the deviations can help define a mathematical model to represent the defect.

In the aviation industry, the measurement of long tubes is carried out using a part holder which consists of fixed supports. For such a case, Mounaud [START_REF] Mounaud | Integrating the flexibility of components in the assembly of aeronautics hydraulic systems[END_REF] applied an upward force to each section in order to cancel the effect of gravity. Radvar-Esfahlan [START_REF] Radvar-Esfahlan | Nonrigid geometric metrology using generalized numerical inspection fixtures[END_REF][START_REF] Radvar-Esfahlan | Robust generalized numerical inspection fixture for the metrology of compliant mechanical parts[END_REF] did the same using a finite element calculation.

The scientific issue

The main problem consists in identifying a local geometric defect (e.g. a hollow or a bump) in a globally deformed part and in proposing a functional tolerancing methodology according to ISO standards.

First, we will introduce the approach through the example of a flat part which is identified using a cloud of points obtained with a Coordinate Measuring Machine (CMM). Then, we will extend the method to the case of a complex surface, such as a blade, which is measured using a large cloud of points.

Identification methodology for the two types of defects

Digital straightening

Notations

The following notations will be used throughout the paper: Pi a point on the nominal (CAD) surface Mi a measured point Ei the initial deviation between the measured point Mi and the point on the nominal surface 𝒛 vector of the beam's axis 𝒅𝑷 𝒊 (dxi,dyi,dzi) vector displacement of point Pi dues to digital straightening 𝒏 𝒊 (nxi,nyi,nzi)

vector normal to the nominal surface at Pi ei the residual deviation after straightening: ei = Ei -𝒅𝑷 𝒊 𝒕

• 𝒏 𝒊

Part identification

Initially, the part being studied was deliberately distorted by tightening it in a vise and applying a torque using a wrench.

Then, the distorted part was measured using a CMM. A coordinate system was defined based on 14 points in the built-in zone (the green zone on the left), and 154 points were measured in the distorted zone (the red zone on the right of Figure 2). Considering the distortion applied to the test part, the model (Fig. 3) consists of a plane zone on the left and, on the right beyond line D, the superposition of a parabolic model representing bending and a helicoid representing torsion.

Thus, the deformed shape is characterized by the deflection f and the angle a at the end of the part. The designer can choose the equations which express the deflection and the rotation.

𝑑𝑦 𝑖 = 𝑓.
𝑧 𝑖 ² 𝐿²

The angular deflection of section zi is: 𝛼 𝑖 = 𝑎.

𝑧 𝑖 𝐿

The displacement in y direction of a point Pi(xi,yi,zi) due to torsion is :𝑑𝑦 𝑖 = 𝑎.

𝑥 𝑖 .𝑧 𝑖 𝐿

The general expression of the deformed shape is: 𝑑𝑦 𝑖 = 𝑓.

𝑧 𝑖 ² 𝐿² + 𝑎.

𝑥 𝑖 .𝑧 𝑖 𝐿

Under the small displacement torsor assumption [START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF], the displacement 𝑑𝑦 𝑖 is a linear function of parameters f and a.

Thus, the objective is to find the parameters f and a which best represent the part.

The association principle

We chose the least-squares method to carry out the association between the actual part and the model.

The case of Fig. 2 was implemented as an Excel worksheet using the deformation model defined in 2.1.3. The measured point coordinates were used to calculate the straightened coordinates:

𝑦 𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙 = 𝑦 𝑚𝑒𝑎𝑠𝑢𝑟𝑒𝑑 -𝑓. (𝑧 𝑚𝑒𝑎𝑠𝑢𝑟𝑒𝑑 ) 2 𝐿 2 -𝑎.
(𝑥 𝑚𝑒𝑎𝑠𝑢𝑟𝑒𝑑 ). (𝑧 𝑚𝑒𝑎𝑠𝑢𝑟𝑒𝑑 ) 𝐿 The objective of the least-squares method is to minimize the sum of the squared residual deviations with respect to the nominal plane (0, x, z). In order to do that, we introduced a cell containing the sum of the squared y residual into the Excel file, then minimized that cell by varying the parameters a and f thanks to the Excel solver:

𝐺𝑜𝑎𝑙 = ∑(𝑦 𝑟𝑒𝑠𝑖𝑑𝑢𝑎𝑙 ) 2

Results of the digital straightening op eration

The results are shown in Fig. 4. 4a shows that the part can be viewed as a twisted beam with maximum deviation Ei max = -0.616 mm.

Fig. 4b was obtained after digital straightening with an angle a = -0.0192 rad ≈ -1.1 ° and a deflection f = -0.037 mm. The residual deviations after digital straightening correspond to local defects. In that configuration, the maximum deviation was ei max = 0.192 mm. The fact that the residual deviation ei max is smaller than Ei max proves that the proposed model represents the part's actual shape quite well.

The deflection thus applied enables one to identify the bending defect. The negative value of f along the y-axis shows that in the actual shape one end points downwards (Fig. 4a). The negative angle represents a negative twisting about the z-axis.

Evolution of the sections

The residual deviations at each point after digital straightening characterize the part. The principle of the method consists in associating each section with a nominal profile and using the least-squares method to minimize the deviation between the nominal profile and the measured points (Fig. 5). Each section is characterized by two parameters, a translation w of the center of the section and a rotation Fig. 5. Deviations of a section after straightening

The evolution of the two parameters w and δ of several sections of the part after digital straightening (Fig. 6) gives a good indication of the defects along the z-axis of the part. In that particular case, the evolution is smooth, but this figure could also have shown a wave or a jump, which could have been unacceptable. Fig. 6. Evolution of the position and orientation parameters of the part's sections Today, there is no specification enabling one to tolerance this continuity of the evolution of sections. A general definition of the "straightness of the part's axis" could be introduced at the center of each section, which would require the construction of the axis of a part of any shape.

Tolerancing of the part

The ISO 10579 standard [START_REF]Geometrical product specifications (GPS) -Dimensioning and tolerancing -Non-rigid-parts[END_REF] enables one to define a protocol for straightening a part. This straightening is usually achieved by applying loads, e.g. by holding the part with clamps. The originality of our approach is that straightening is not performed physically, but digitally. Therefore, tolerancing must be expressed by an equation which characterizes this digital straightening (Fig. 7). Consequently, this makes it possible to specify the shape of a side, of a section or of the entire part.

The tolerances marked "F" (for "free") must be satisfied in the free state, i.e. without straightening, with deviations Ei. The unmarked tolerances must be satisfied in the stressed state after digital straightening, i.e. with the deviations ei (Fig. 8). Fig. 8. Specifications relative to the free state and after digital straightening

Other examples of distorted models

Regardless of the part's configuration, its global defects must be defined using a parametric model. For a cylindrical part, it is easy to define a model using 4 parameters: -the cone angle -the base radius -the ovality -the radial deflection at the center of the generatrix. For large beam-like parts with constant cross sections, many bending models can be imagined depending on the loading cases and boundary conditions applied to the beam. This work proposes two deformation modes for this type of part:

Mode 1: (Fig. 9) Built-in condition at A and deflection at B with a readjustment of the beam's direction and the base coordinate system through two parameters, decY(bottom) and decY(top). By similarity with the bending of a beam with constant cross section, the deformed shape is a polynomial of degree 3.

Mode 2: (Fig. 9) Pinned connections (i.e. no translations) at ends A and B, and a deflection f at the center of the beam. Twisting of the part around the z-axis with angles a at B and 0 at A (Fig. 9 and10).

Fig. 9. Examples of deformation modes

For the second mode, a twist angle is applied at B. The resulting shape is a helicoid (Fig. 10). If the angle is too large, it is necessary to determine the displacements of some points along the x and y axes because the small displacement assumption [START_REF] Bourdet | The concept of the small displacement torsor in metrology[END_REF] is not satisfied.

Fig. 10. A twisted model

For each model, the equations giving the displacement at each point must be defined using the following parameters:

-a the rotation angle -fx, fy the deflections in x and y direction (Fig. 9)

Mode 1:

The displacement of point Pi is decomposed in x and y directions: 𝑥 𝑖 -𝑥𝑐 𝐿 Other models can be created, e.g. for parts with varying sections. The finite element method can be used to calculate the deflection at each point under a nominal force or moment [START_REF] Radvar-Esfahlan | Nonrigid geometric metrology using generalized numerical inspection fixtures[END_REF], [START_REF] Radvar-Esfahlan | Robust generalized numerical inspection fixture for the metrology of compliant mechanical parts[END_REF]. Due to the linearity assumption, one obtains the deformed shape as a function of the force and moment.

𝑑𝑥 𝑖 = 𝐷𝑒𝑐𝑋 (𝑏𝑜𝑡𝑡𝑜𝑚)
By comparing several calculation criterions, one can select parameter values which minimize the sum of the deviation absolute values

Application to complex industrial surface

Application to a turbine blade

In the context of the manufacturing of aeronautical parts, Snecma wishes to improve its understanding of the deviations measured on turbine blades (Figure 11). Indeed, the measured deviations can be due to various types of defects on different scales. For example, a 0.05 mm distortion along the trailing edge has a greater impact on performance than a 0.3 mm global defect due to twisting of the part. However, the former generates turbulences, whereas the latter only deflects the air flow.

Our approach was applied to a blade designed by Snecma using a cloud of points, spanning the entire part (about 600,000 points), obtained by means of an optical measuring tool [20].

The main objective was to quantify the deviations Ei between the actual part and the nominal geometry, then the deviations ei between the actual part and the distorted geometry.

Identification of the global form defects

Two studies were carried out using the deformation modes defined previously: -clamping of the blade at its base in the turbine wheel, with no twisting (Fig. 11 The measurements were performed by Snecma using an optical system developed by GOM [20], which consists in a stereovision device projecting an array of blue fringes (to avoid interaction with the natural light [START_REF] Savio | Metrology of freedom shaped parts[END_REF]).

The points were measured in a measurement coordinate system. The first fitting of the point cloud with the CAD model was achieved by considering several surfaces and matching 6 target points with the corresponding CAD surfaces (Fig. 12). These surfaces were not machined and provided a constant coordinate system throughout the manufacturing process. This is the raw reference frame. The time required to perform these operations (acquisition of the point cloud and registration) was about 10 minutes for an accurate view of the whole part with triangulation, but the initial setting is time-consuming (about 4 hours), so profitability is very dependent on the number of parts measured.

Point extraction for straightening

The objective is to determine the straightening parameters, particularly deflection and twisting, to be applied to the point cloud in order to minimize the sum of the squared deviations. The Excel solver was used so the operation would remain as generic as possible.

In order to avoid numerical problems and excessive computation time because of too many points, the calculation was performed on only a few sections (Fig. 13).

Fig. 13. The measured cloud and the sample extracted

The least-squares calculation requires that the real point Mi coordinates and the normal 𝒏𝒊 components of the surface are known.

The construction of the sample is illustrated in Fig. 14. The optical measurements lead to a cloud of real points and a triangular mesh connecting these points. The expected sections are generated in the CAD system by extracting a small number of nominal points Pi along the section line (38 points per section in this particular experiment). At each point 𝑃𝑖 of a nominal section, the deviation Ei between the mesh and the CAD model is defined along the normal 𝒏𝒊. The lines (Pi, 𝒏𝒊) and the mesh intersect at points Mi considered as real points:

E i = 𝑷 𝒊 𝑴 𝒊 𝒕 • 𝒏 𝒊
The points along the trailing edge are excluded from the calculation because the local radius is too small and the normal vectors are not constant with regard to deviation dues to twist.

Thus, each section is characterized by a series of nominal points Pi from the CAD model, the normal 𝒏𝒊 at each point and the deviation Ei with respect to the actual surface materialized by the mesh of the measured points.

This extraction was automated through a quick development in the CATIA CAD software. As written before and in order to avoid numerical problems due to the curvature, the sample distribution was created without taking into account the leading and trailing edges. At each point Pi, the deformation models lead to 𝑑𝑥 𝑖 and 𝑑𝑦 𝑖 according to chosen mode. Then, the Excel solver changes the values of the parameters in order to minimize the sum of the squared deviations. The computation time is only a few seconds, depending on the size of the sample.

Another calculation could minimize the sum of the shape defects in the 9 sections by minimizing:

𝐷 = 2 ∑ 𝑀𝑎𝑥(|𝑒 𝑖 |) 9 𝑗=1
One advantage of that method would be that the calculation would be unaffected by the density of the point cloud.

Digital straightening of the actual part

With the values obtained, a displacement field was applied to the whole set of the measured points.

Then, the new cloud of points obtained after digital straightening was imported into CATIA for comparison with the nominal model. Finally, the residual deviations (local defects) were obtained through a deviation analysis and plotted as a color map (Fig. 15).

The values obtained (hidden for confidentiality) show the global shape defects of the blade (deflection and twisting).

The color map of the residual deviations characterizes the local defects of the surface. In our example, the main defects, before digital straightening, are a lack of material on the suction side of the leading edge and an excess of material on the pressure side. After digital straightening we note that some deviations was due to global defects. As explained in Section 2.4 in the case of a simple example, this geometric analysis method can be specified with tools of the ISO 10579 standard [START_REF]Geometrical product specifications (GPS) -Dimensioning and tolerancing -Non-rigid-parts[END_REF], leading to the definition of a stressed state.

Fig. 16 shows the blade case with deformation Mode 1. The coordinate system is defined directly on the CAD model created from the measurement cloud.

As mentioned before, the specifications can be controlled at two different states (stressed or free), depending on the presence of modifier F. So there are two dedicate check procedures for both case.

The first one, with F modifier, consists to verify if the measurement cloud is included in the tolerance zone, defined by the ISO specification, immediately after the measurement clouds registration on the datum reference frame of the CAD model. In that case no digital straightening is realized, this is the "free" state.

In the second case, without F modifier, the measurement cloud is fitted on the CAD model next digitally straighten in a "stressed" state defined by defect modes, then straighten points cloud is compared to the tolerance zone.

The last two figures correspond to the industrial requirement to analyze the defects and validate the conformity of a blade.

Conclusion

The digital straightening method presented in this work can be applied to slender parts with constant or variable sections, or even to thin parts such as sheet metal or crankcases, once the deformation model has been chosen.

The ISO tolerancing proposed is based on a dual specification in the free state and in the straightened state through a digital process, leading to the deformation equation and the objective function to be minimized.

In the case of the simple part discussed in section 2, the measurements were carried out by means of a measuring machine. For the blade, the measurement of a large point cloud was performed using an optical device.

In both cases, the model's parameters characterized the general shape of the part. The residual deviations after straightening led to the representation of the local defects.

Several types of straightening models can be tested in order to identify the global shape deviations and the smaller residual deviations of the part. In addition, the smoothness of the variation in the directions of the part's sections can be described.

The residual defects can be easily visualized as color maps. The corresponding computation time is small (only a few seconds), but the time required for the measurements remains too long (4 hours setup time + 10 minutes per part for the measurements and triangulation). Today, this method can be used only to qualify a manufacturing process or to perform a sample measurement during production.

New portable optical measurement systems should make it possible to measure parts directly in the workshop and to check raw parts and machined workpieces in only a few seconds.

However, additional effort is required in order to be able to adjust the density of the point clouds for determining the straightening parameters, especially at the leading and trailing edges, using the least-squares criterion and comparing the results with those of the minimization criterion of the sum of the sectional shape defects.

The demonstrator which has been achieved meets the needs of Snecma Company. The company decided to integrate this application into the CATIA environment to use it in the central metrology services, especially to validate the production process of raw or finished blades.

Acknowledgement

Fig. 2 .

 2 Fig. 2. The measurement of the part

Fig. 3 .

 3 Fig. 3. The deformation model For example, the displacement in y direction due to bending in the (O, 𝒙, 𝒚) plane on point Pi(xi,yi,zi) is:

Fig. 4 .

 4 Fig. 4. Deviations Ei of the plane part (a) and deviations ei after digital straightening (b) Figure 4a shows that the part can be viewed as a twisted beam with maximum deviation Ei max = -0.616 mm. Fig. 4b was obtained after digital straightening with an angle a = -0.0192 rad ≈ -1.1 ° and a

Fig. 7 .

 7 Fig. 7. Tolerancing of the part

  Fig. 11. The blade model

Fig. 12 .

 12 Fig. 12. The fitting of the cloud of points

Fig. 14 .

 14 Fig. 14. Construction of the sample for the calculation of the straightening parameters

Fig. 15 .

 15 Fig. 15. Comparison before / after digital straightening depending on the mode selected

Fig. 16 .

 16 Fig. 16. Blade tolerancing by digital straightening with ISO 10579

The authors wish to thank Snecma for welcoming them in an internship at Gennevilliers and Villaroche in 2012 and for providing the measurements used in their study. This work was carried out under the Quick GPS project, labeled System@tic and Astech, funded by the "Region Ile de France".

Authors

Mathieu PETITCUENOT has a doctorate in mechanical engineering. He has worked in the field of geometric tolerancing and manufacturing specifications of complex surfaces, as turbine blade, at Snecma company.

Bernard ANSELMETTI is a professor in mechanical engineering internationally recognized for his works in ISO specification, metrology and manufacturing tolerancing. He notably developed software for automatic specification of mechanical parts based on contacting feature analysis.