Charlotte Derouet
email: charlotte.derouet@espe.unistra.fr

Gaëtan Planchon

Thomas Hausberger

Reinhard Hochmuth

Gaetan Planchon²

Thomas Hausberger²

Alexander Grothendieck

Bridging probability and calculus: the case of continuous distributions and integrals at the secondary-tertiary transition

Keywords: Transition to and across university mathematics, Teaching and learning of analysis and calculus, Teaching and learning of probability, Anthropological Theory of the Didactic

published or not. The documents may come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Continuous probability distributions (CPD) and integral calculus (IC) are two topics that are taught in France during the last year of high school (grade 12 of the scientific track). They constitute two sectors (in the sense of the Anthropological Theory of the Didactics, ATD) that belong to the two different but closely related mathematical domains of probability theory and calculus respectively. Indeed, the continuous probability of an event and the definite integral with respect to a non-negative function are both defined as areas of suitable two-dimensional domains in the syllabus (in France), and the formula 𝑃(𝑎 ≤ 𝑋 ≤ 𝑏) = ∫ 𝑓(𝑥)𝑑𝑥 𝑏 𝑎 is the key for solving several standard tasks in CPD, where 𝑋 represents a random variable and 𝑓 its associated density function.

IC is the focus of extended studies in mathematics education: for instance, research [START_REF] Schneider | Des objets mentaux "aire" et "volume" aux calculs de primitives[END_REF][START_REF] Luong | dx et ∫ ostensifs pour questionner l'enseignement des intégrales. Montpellier : IUFM de Montplellier[END_REF][START_REF] Haddad | Que retiennent les nouveaux bacheliers de la notion d'intégrale enseignée au lycée ?[END_REF] was conducted in the context of Belgian, French, Vietnamese and Tunisian secondary education, including the secondary-tertiary transition [START_REF] Haddad | Que retiennent les nouveaux bacheliers de la notion d'intégrale enseignée au lycée ?[END_REF]. By contrast, there is hardly any literature on CPD and, to make things worse, the available studies mostly put the emphasis on normal distributions [START_REF] Batanero | From data analysis to inference: a research project on the teaching of normal distributions[END_REF][START_REF] Wilensky | What is normal anyway? Therapy for epistemological anxiety[END_REF][START_REF] Batanero | Students' reasoning about the normal distribution[END_REF][START_REF] Pfannkuch | Reasoning about distribution: A complex process[END_REF]. Therefore, the teaching-learning phenomena generated by the interrelationship between CPD and IC are still to be investigated.

A first stone was laid by [START_REF] Derouet | How can histograms be useful for introducing continuous probability distributions?[END_REF] see also [START_REF] Derouet | La fonction de densité au carrefour entre probabilités et analyse en terminale S. Etude de la conception et de la mise en oeuvre de tâches d'introduction articulant lois à densité et calcul intégral[END_REF], who studied possible ways to introduce the density function at grade 12 so that students may construct this concept starting from considerations regarding histograms and therefore might relate continuous probability to the integral. By an analysis of textbooks, Derouet could show that the two sectors CPD and IC are very much compartmentalised in the French curriculum (Derouet, 2016, pp. 127-190). For instance, the above formula is seldom justified by a thorough discussion of the definitions involving areas, which certainly hinders the bridging of the two sectors by students.

In this paper, we regard this "compartmentalisation" of knowledge as an institutional phenomenon and therefore use ATD as the theoretical framework (see below). Our goal is to study the impact of this compartmentalisation on the learning of mathematics: are French students able to mobilize the sector of IC to solve tasks in the sector of CPD and vice versa at the transition from high school to higher education?

After the presentation of theoretical constructs from ATD used in this research, we will describe the reference epistemological model (REM) that we elaborated for the types of tasks in CPD and IC with regard to studying interrelations of the two sectors. We will then describe our methodology that builds on the elaboration of a questionnaire, based on the REM, that has been submitted to students at the entrance of university. We finally present results of a primary analysis of the data from the questionnaire and draw some conclusions and perspectives opened up through this study.

THEORETICAL CONSTRUCTS

ATD "postulates that any activity related to the production, diffusion or acquisition of knowledge should be interpreted as an ordinary human activity, and thus proposes a general model of human activity built on the key notion of praxeology" [START_REF] Bosch | Introduction to the Anthropological Theory of the Didactic[END_REF]. The praxeology 𝛱 is represented by a quadruple [𝑇/𝜏 /𝜃/𝛩]: its praxis part (or know-how) consists of a type of tasks 𝑇 together with a corresponding technique 𝜏 (useful to carry out the tasks 𝑡 ∈ 𝑇 in the scope of 𝜏). The logos part (or know-why) includes two levels of description and justification: the technology 𝜃, i.e. a discourse on the technique, and the theory 𝛩, which often unifies several technologies.

The elaboration of a reference epistemological model [START_REF] Florensa | The epistemological dimension in didactics: Two problematic issues[END_REF] as sequences of praxeologies, for a given body of knowledge, is an important step in any research carried out in the ATD framework. It is the tool that will be used by the researcher to describe, analyse, put in question or design the specific contents that are at the core of a teaching and learning process. In order to build such a model, "mathematical praxeologies are described using data from the different institutions participating in the didactic transposition process, thus including historical, semiotic and sociological research, assuming the institutionalized and socially articulated nature of praxeologies" (loc. cit. p. 2637).

Our study relies on an overview of standard textbooks used at grade 12 in France, as well as the official syllabus, in order to identify the standard praxeologies in CPD and IC that may be related and the nature of this relationship at the praxeological level. An epistemological investigation of the historical development and the interrelation of both domains have previously been carried out in (Derouet, 2016, pp. 67-85). In order to test the effect of the institutional compartmentalisation of knowledge on the learning of mathematics, we need to check the availability of the identified praxeologies in the praxeological equipment of students, and then submit tasks to students which need to bridge CPD and IT as mutually interdependent sectors that share techniques or technologies, borrow them from or lend them to each other. Special care must be taken in the phrasing of these bridging tasks, taking into account the effect of ostensives [START_REF] Bosch | La sensibilité de l'activité mathématique aux ostensifs[END_REF], that is to say the role of signs. Indeed, ostensives contribute to the activation of the specific sectors to which they belong and therefore direct students toward specific techniques.

REFERENCE EPISTEMOLOGICAL MODEL

Even though problems of quadratures arose in ancient Greece, IC finds its roots as a systematic method in the 17th century. The emergence of continuous probability may be situated in the 18th century with the theory of errors in physical measurements. Various functions were introduced to model the distribution of errors and the area under the curve permitted to evaluate the "theoretical frequency" (so the probability) of the deviation from the "true" value. CPD was thus naturally connected to IC in its historical roots. The gaussian distribution was proposed later by Gauss in 1809.

To identify the different praxeologies, we have analysed 12 textbooks of the grade 12 of the scientific track (edition 2012). We focused on the exercises with a given solution in the textbooks to have access to the usual techniques for the different tasks.

In our study, we will focus on two main types of mathematical tasks 𝑇 𝐼 and 𝑇 𝑃 , which are related to the mathematical domains of integral calculus and the continuous probability respectively:

• 𝑇 𝐼 : compute a value for an integral ∫ 𝑓(𝑥)𝑑𝑥 𝑏 𝑎 for a positive continuous function 𝑓;

• 𝑇 𝑃 : determine the probability 𝑃(𝑎 ≤ 𝑋 ≤ 𝑏) for a random variable 𝑋 endowed with a density function 𝑓. The type of tasks 𝑇 𝐼 may further be split into two subtypes of tasks, depending on the expected result: an exact value (𝑇 𝐼,𝑒𝑥𝑎𝑐𝑡 ) or an approximation (𝑇 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 ).

The most useful technique 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 to solve 𝑇 𝐼,𝑒𝑥𝑎𝑐𝑡 is to compute a primitive of the function and apply the fundamental theorem of calculus. The corresponding technology 𝜃 𝐼,𝑒𝑥𝑎𝑐𝑡 is given by the fundamental theorem of calculus that relates integrals and primitives: ∫ 𝑓(𝑥)𝑑𝑥 𝑏 𝑎 = 𝐹(𝑏) -𝐹(𝑎) with 𝐹′ = 𝑓. The theory 𝛩 𝐼,𝑒𝑥𝑎𝑐𝑡 includes the definition of the definite integral for a continuous positive function as an area and properties of areas that may be formalised into a local axiomatic theory 1 .

The technique 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 thus resorts to praxeologies dedicated to the computation of primitives. The standard technique at high school level is to use the "tabular of primitives" (deduced from the tabular of derivatives). The technology comprises the properties of the derivative and the theory is that of differential calculus. In the case of piecewise affine functions, an alternative technique to 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 is to interpret the integral as the area of an elementary surface (or a union of these).

The type of tasks 𝑇 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 may be solved using two main techniques: using a calculator (or software), more or less a blackbox, or applying the "rectangle method". The latter technique 𝜏 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 consists in considering the integral as an area, taking a subdivision of the interval of integration and computing the sum of rectangular areas. The technology 𝜃 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 comprises the definition of the integral and properties of areas. A further theoretical level 𝛩 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 is mainly non-existent at high school level (cf. endnote 1).

Regarding the type of tasks 𝑇 𝑃 , two cases need to be distinguished, depending on whether a primitive of the density function f is known (𝑇 𝑃,𝑝𝑟𝑖𝑚 ) or not. The latter case is reduced to that of the normal distribution (𝑇 𝑃,𝑛𝑜𝑟𝑚 ), which is dealt with using the implementation of such a distribution in a calculator or software. Computer scientific tools are mainly used as a blackbox by students, which hinders the possibility for students to make connections with IC.

The generic technique 𝜏 𝑃,𝑝𝑟𝑖𝑚 for 𝑇 𝑃,𝑝𝑟𝑖𝑚 is to compute ∫ 𝑓(𝑥)𝑑𝑥 𝑏 𝑎 , in other words to resort to the praxeology 𝛱 𝐼,𝑒𝑥𝑎𝑐𝑡 . The technology 𝜃 𝑃,𝑝𝑟𝑖𝑚 is given by the formula

𝑃(𝑎 ≤ 𝑋 ≤ 𝑏) = ∫ 𝑓(𝑥)𝑑𝑥 𝑏 𝑎
, and the theory 𝛩 𝑃,𝑝𝑟𝑖𝑚 comprises the definition of a continuous probability (as the area of the corresponding domain) and the definition of a probability density function. At high school, two particular cases are emphasised and lead to local techniques, as concrete formulas are available for 𝑃(𝑎 ≤ 𝑋 ≤ 𝑏) in the case of the exponential and uniform distributions. For instance, the technique 𝜏 𝑃,𝑒𝑥𝑝 in the case of the exponential may be reduced to computing 𝑒 𝜆𝑎 -𝑒 𝜆𝑏 with the technological argument 𝑃(𝑎 ≤ 𝑋 ≤ 𝑏) = 𝑒 𝜆𝑎 -𝑒 𝜆𝑏 .

Let us recall that our model is based on the study of standard textbooks used in grade 12 classes in France and is dedicated to the description of the teaching-learning of CPD and IC as it actually is (we are not planning task-design at this stage of the research).

In this model, we note the following links between CPD and IC: at the level of the theoretical blocks, praxeologies in both sectors are anchored on the empirical notion of area. At the level of the technique, 𝛱 𝑃,𝑝𝑟𝑖𝑚 uses 𝛱 𝐼,𝑒𝑥𝑎𝑐𝑡 , so that, from an ecological point of view (Bosch & Gascon, 2014, p. 72), CPD contributes to the thriving of such IC praxeologies. By contrast, 𝛱 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 does not seem to be reinvested in CPD (whereas the computation of a probability of the Gaussian distribution could be an opportunity to mobilize 𝛱 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 ). Conversely, we did not detect elements of the praxis of CPD in the IC sector. This isn't a surprise: IC is regarded as a prerequisite to CPD and precedes the teaching of CPD in all textbooks. Nevertheless, the normal distribution is a prototypical example of a function whose primitive cannot be expressed in terms of available elementary functions. This fact explains the choice of techniques in 𝛱 𝑃,𝑛𝑜𝑟𝑚 and contributes also to the logos of 𝛱 𝐼,𝑒𝑥𝑎𝑐𝑡 (by complementing the statement that every continuous function admits a primitive). What about the type of tasks 𝑇 𝐼,𝑛𝑜𝑟𝑚 : compute ∫ 𝑒 -(𝑥-𝑚) 2 /2𝜎 2 𝜎√2𝜋 𝑑𝑥 𝑏 𝑎 ? It could appear in the IC sector through 𝛱 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 only but we didn't find it in any textbook and it is never stated as such in the CPD sector. The most efficient technique requires to use the formula 𝑃(𝑎 ≤ 𝑋 ≤ 𝑏) = ∫ 𝑓(𝑥)𝑑𝑥 𝑏 𝑎 from right to left: although the equality is symmetrical as an equivalence relation, it isn't symmetrical as a sign which denotes a succession of operations in performing a computation. How would students react to such a task that asks to bridge CPD and IC in an unusual way? This question came to us as a starting point for the elaboration of our questionnaire dedicated to the investigation of the educational effects of the institutional phenomenon of compartmentalisation of knowledge, in the case of CPD and IC.

THE QUESTIONNAIRE

Our main goal is to test whether students are able or not to connect CPD and IC, and especially mobilize the CPD sector to solve an IC task when ostensives do not indicate explicitly the probability domain. To do so, we have elaborated a questionnaire both to check the availability of standard praxeologies of CPD and IC in the praxeological equipment of students and the capacity of students to complete such bridging tasks.

Bridging tasks appear at the very end of the questionnaire and are stated as follows: . It suffices to give an idea of the method if its implementation is too complicated. Question 6 contains an instance of the type of tasks 𝑇 𝐼,𝑛𝑜𝑟𝑚 discussed in the REM. The task is stated in an opened way, asking for every method that students may know to compute an exact approximate value for the Gaussian integral. Praxeologies 𝛱 𝐼,𝑒𝑥𝑎𝑐𝑡 and 𝛱 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 should therefore also show up. In question 7, we intend to check whether students can say that A is both a primitive for f and a probability associated with the density function f, or restrict to the IC sector with an interpretation in terms of areas.

Previous questions intend to "activate" both sectors CPD and IC equally. In this respect, question 1 offers a routine task of type 𝑇 𝐼,𝑒𝑥𝑎𝑐𝑡 in the case of a straightforward exponential function. Analogously, the first part of question 4 is a routine task of type 𝑇 𝑃,𝑒𝑥𝑝 (compute 𝑃(1 ≤ 𝑋 ≤ 5) when X has an exponential distribution of parameter 3). In its second part, students are asked for a graphical interpretation of the probability 𝑃(1 ≤ 𝑋 ≤ 5) that is in fact defined as an area in the high-school syllabus, as well as the integral: this interpretation is therefore essential in order to link the logos of 𝛱 𝑃 and that of 𝛱 𝐼 .

Question 2 activates the CPD sector by soliciting an element of the logos of 𝛱 𝑃 , namely the properties that define a density function. This logos is crucial in the bridging question 7, which is stated in the IC sector without any reference to CPD. Question 3 tests if students are able to retrieve the definitions of both the exponential and normal distributions by specifically asking for those in the case of simple parameters (the reduced centered gaussian law). The latter is an element of the logos of 𝛱 𝑃,𝑛𝑜𝑟𝑚 : we wish to check if students are able to recognize the normal distribution in the statement of the bridging question 6, while taking care not to direct them towards a specific technique (hence the order of questions).

In question 5, we rather activate the IC sector, more precisely elements of the logos of 𝛱 𝐼,𝑒𝑥𝑎𝑐𝑡 (primitives), but we have in mind the praxeology 𝛱 𝑃,𝑛𝑜𝑟𝑚 in relation to the bridging question 6: students are asked to provide an example of a continuous function, if it exists (or justify the impossibility), that a) doesn't possess primitives b) admits primitives but expressions for these are not "explicitly known" (expected answer: the density function of the normal law).

Summarizing, by the questions 1 to 4 we want to investigate whether the students master techniques 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 in a calculus context and 𝜏 𝑃,𝑝𝑟𝑖𝑚 in a probability context and whether they know technologies related to IC and CPD. Then, by questions 6 and 7, we focus on relationships between CPD and IC and the previous questions: we analyse links between the questions 3, 4 (second part), 5 and 6, on the one hand, and links between the questions 2 and 7, on the other hand.

DATA ANALYSIS AND RESULTS

The questionnaire was used at the beginning of September 2017 (the first week of classes) in two classes of first year CPGE (French engineers school preparatory classes) students, which is in fact at the transition between secondary and tertiary levels. The first class (called class N) is a class of MPSI (Mathematics, Physics and Engineering Science) and the second class (called class R) is a class of PCSI (Physics, Chemistry and Engineering Science) of a rather prestigious establishment. The students working on the questionnaire are in selective classes, so we can assume that they are "good" scientific students, and in particular, if they meet difficulties then these are shared by the other students. We only analysed answers from students who studied in French high school during the past year because we constructed the questionnaire taking into account the context of the French high school institution. We retrieved 82 questionnaires (40 of the class N and 42 of the class R). Except for a few students (less than 5), the students didn't use a calculator during the test.

From the 82 students, only one does not mobilize the technique 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 to resolve the routine task concerning IC (question 1). 85% find a correct expression of the primitive and 78% obtain the correct result, which means that this technique is mastered quite well by students. Regarding the computation of the probability for an exponential distribution (question 4), 63% of the students get a correct result. 82% of the students use the technique 𝜏 𝑃,𝑝𝑟𝑖𝑚 and the other students directly apply a formula. More than 70% of the students could identify the probability as an integral. Summarizing, except for some errors regarding the primitive or the computation, the majority of students are able to pass from a probability to an integral and, moreover, know the fundamental theorem of calculus. In praxeological terms, they are able to mobilize the technique 𝜏 𝑃,𝑝𝑟𝑖𝑚 that implies the use of the technique 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 in the case of an exponential distribution.

Regarding CPD and neglecting formulation and formalization issues, only 39% of the students know the definition of the density function (question 2) and only 27% recollect the density function of the normal distribution (question 3b). In view of question 5, only 32% mobilize the theorem claiming that all continuous functions on an interval admit a primitive. 52 % give an example of a function for which they don't explicitly know a primitive (although it might exist and be expressed in terms of standard functions 2 , for instance 𝑙𝑛(𝑥)). Among these, 15% mention the density function of the Gaussian distribution (or a function of the type 𝑒 -𝑥² ). Of the 22 students who know the density function of the normal distribution, half propose it as an example in question 5 (13% of all the students).

9% of the students do not answer question 6. The method most often proposed is technique 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 (46%). Less than 20% (16 students) mention the normal distribution and 16% propose the rectangle method (𝜏 𝐼,𝑎𝑝𝑝𝑟𝑜𝑥 ). A few students propose the technique "integration by parts", which is beyond the curriculum in grade 12. Moreover, only 23% of the students who propose the rectangle method are able to illustrate the method by drawing the graph of the Gaussian curve (the other students draw a wrong curve or do not consider any graph). Only 11 of the 16 students (69%) who mention the normal distribution in this question write that the integral I is equal to the probability 𝑃(-0,5 ≤ 𝑋 ≤ 1) with 𝑋 a random variable of a reduced normal distribution and 7 of them state more precisely that they have to use the calculator to evaluate this probability (𝜏 𝑃,𝑛𝑜𝑟𝑚 ). So, our results indicate that the ostensive

∫ 𝑒 -𝑥 2 /2 √2𝜋 𝑑𝑥 1 -0,5
without indication invite students to stay in the IC sector and even more particularly in the praxeology 𝛱 𝐼,𝑒𝑥𝑎𝑐𝑡 (even if it is not possible here), which reflects that 𝛱 𝐼,𝑒𝑥𝑎𝑐𝑡 is the praxeology most developed in grade 12. To pass from an integral to a probability and to change the sector in this direction does not seem to be natural for students. Moreover, of the 22 students who know the expression for the density function of the reduced normal distribution (question 3b), only 10 (45%) recognize the density function of the normal distribution is this context. Of the 15 students who say that the Gaussian function does not admit primitives expressed by using standard elementary functions (question 5b), around 27% proposes to use the technique 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 nevertheless and only 53% recognize the normal distribution in question 6. This means that most of the students were not able to mobilize 𝛱 𝑃,𝑛𝑜𝑟𝑚 in the IC-context of question 6, that is 𝛱 𝐼,𝑛𝑜𝑟𝑚 . In particular the application of 𝜏 𝐼,𝑒𝑥𝑎𝑐𝑡 , although question 5b is answered correctly, demonstrates the strong compartmentalisation between CPD and IC.

Regarding the answers to question 7, we notice that 38% of the students (31) state that the function 𝑓 is a density function. 42% of them (13 students) could justify that it is a density function including 4 who forget to mention the positivity of the function (because they don't write this condition in their definition in question 2). Moreover, 12 of the 31 students as well as one additional one identify 𝐴(𝜆) as 𝑃(0 ≤ 𝑋 ≤ 𝜆). 8 students think that the distribution in question 7 is an exponential distribution of parameter 𝑥 and 7 students state that 𝐴 is a density function. Finally, more than 46% of the students mention "probability", which means that they manage to identify at least some link between the IC embedding of question 7 and CPD. Probably, the questionnaire itself influenced students and the percentage would be lower otherwise, i.e. if question 7 was asked independently. Overall correct results with justification are rare and of the 32 students who master the definition of a probability density function (one element of the technology 𝜃 𝑃,𝑝𝑟𝑖𝑚 tested in question 2), 31% (10 students) do not identify 𝑓 as a density function. 16% of the students do not at all answer question 7.

Again, we observe also with respect to question 7 that CPD praxeologies, although they are in principle available, often cannot be mobilized in the IC contextualization.

Summarizing, the data analysis shows that techniques related to one sector are available for the majority of students only when ostensives related to this sector are provided. Additionally, related technologies are much less mastered by students. Perhaps, this could be an explanation why it is not natural for students to mobilize praxeologies of the CPD sector for a task in the IC sector, in addition to the fact that these tasks are not taught in the classroom. The data analysis by all means shows a strong compartmentalisation between CPD and IC.

CONCLUSION AND PERSPECTIVES

The results of our primary data analysis clearly demonstrate a strong compartmentalisation between CPD and IC. In particular, techniques from CPD, although available in a CPD task, could not be mobilized in an IC-contextualized task. A next step in our research will be more detailed data analyses looking for correlations and interdependencies between techniques and technologies of CPD and IC. We further observed that available techniques were not accompanied by related technologies. One could claim that more elaborated technologies might support the transfer of techniques.

More generally, we think about studies investigating the impact of changes in the institutional setting, i.e. establishing innovative teaching sequences with less compartmentalisation. A teaching sequence articulating CPD and IC is proposed in [START_REF] Derouet | La fonction de densité au carrefour entre probabilités et analyse en terminale S. Etude de la conception et de la mise en oeuvre de tâches d'introduction articulant lois à densité et calcul intégral[END_REF]. The effect of this teaching on the answers of students to the questionnaire could be analysed by comparing the latter with the present results.

NOTES

1. This axiomatic remains implicit at the secondary level; it may be related to Measure Theory at university level.

2. We realised a posteriori that our question was not phrased properly: "a function that admits a primitive whose expression is not explicitly known to you" may be interpreted as a lack of techniques to actually compute the primitive and not the impossibility to provide an expression (in terms of elementary functions).

Question 6 :Question 6 :

 66 Expliquer toutes les méthodes que vous pouvez utiliser pour déterminer une valeur exacte et/ou approchée de l'intégrale suivante : On pourra se limiter à donner une idée de la méthode si sa mise en oeuvre est trop compliquée. Question 7: Soit A la fonction définie par 𝐴(𝜆) = ∫ 𝑓(𝑥)𝑑𝑥 𝜆 0 avec 𝑓(𝑥) = 𝑥𝑒 -𝑥 , pour tout  ∈ [0; +∞[. On peut démontrer que lim 𝜆→+∞ ∫ 𝑥𝑒 -𝑥 𝑑𝑥 𝜆 0 = 1. D'après ce résultat, expliciter tout ce que vous pouvez dire sur la fonction A et la fonction f.Translation: Explain all the methods that can be used to determine an exact and/or approximate value of the following integral:

Question 7 :

 7 Let A be the function defined by 𝐴(𝜆) = ∫ 𝑓(𝑥)𝑑𝑥 𝜆 0 with 𝑓(𝑥) = 𝑥𝑒 -𝑥 , for all  ∈ [0; +∞[. According to this result, note everything that you can say about the function A and the function 𝑓.

Figure 1 :

 1 Figure 1: bridging tasks submitted to students


Martin Bué
email: martin.bue@inria.fr

Diego Cattaruzza
email: diego.cattaruzza@centralelille.fr

Frédéric Maxime Ogier
email: maxime.ogier@centralelille.fr

Frédéric Semet
email: frederic.semet@centralelille.fr

A Two-Phase

HAL is

Introduction

In this paper, we focus on the warehouse management problem faced by our industrial partner Hap-pyChic1 . HappyChic is a French firm specialized in men's ready-to-wear. HappyChic gathers three leading brands: Jules, Brice, and Bizzbee. The logistic department of HappyChic is in charge of the three warehouses (one for each brand), and all the related logistics activities. In this paper, we consider operations in the warehouse dedicated to clothing products of the brand Jules, located in the North of France, close to the city of Lille. The warehouse supplies the shops of Jules in all France and some shops abroad. Thus, the usual customers to satisfy are the shops to which new clothing products are shipped.

The warehouse under consideration is depicted in Figure 1 and consists in four main zones: (1) the reception zone where products delivered by suppliers are received; (2) the storage zone where cartons are stored in racks; (3) the picking zone where pickers pick items from cartons and prepare boxes to be sent to customers; (4) the delivery zone where trucks are loaded with boxes shipped to customers.

All along the article, we make the distinction between carton and box. A carton contains items of a single product. Cartons are received by the warehouse on pallets and stored in the storage zone. They are then moved to the picking zone. There, items are picked from the cartons and put into boxes by pickers. Thus, a box contains customer demands and possibly items of different products.

Trucks carrying products from the suppliers arrive at the reception zone of the warehouse. Here dedicated employees unload the trucks, and pallets of cartons are prepared to be stored. Then, employees driving forklifts move cartons from the reception zone into the racks of the storage zone. Other employees are dedicated to supplying the picking zone. They pick up cartons from pallets in the storage zone and bring them, using forklifts, into the picking zone. They are allowed only to go through cross-aisles 1, 3 and 5 and to leave cartons at the top or the bottom of the different aisles (see solid blue lines in Figure 1).

The picking zone is composed of 4 blocks (see Figure 1). Let us call these blocks, block 1 to 4. Each block is equipped with racks where cartons with products are located. Contrarily to standard cases presented in the literature (see for example [START_REF] Scholz | Order batching and picker routing in manual order picking systems: the benefits of integrated routing[END_REF]), we distinguish the two sides of the racks: the front and the back side. Cartons are inserted into racks from the back side and products are picked from the front side. Racks are situated in such a way that a front (resp. back) side always faces another front (resp. back) side of another rack (except, of course, for the first and the last rack). This placement defines a front-aisle and a back-aisle.

Back-aisles are used by employees to bring cartons from the top or the bottom of an aisle to their exact location (see dotted blue lines in Figure 1). Front-aisles are used by employees called pickers to collect products from cartons.

The cross-aisles in the picking zone are not all walkable by pickers. For security reasons, crossaisles 1, 3 and 5 are dedicated to forklifts that supply the picking zone from the storage zone. Pickers are allowed to walk only across cross-aisles 2 and 4. Moreover, these cross-aisles are equipped with a conveyor belt that automatically brings boxes prepared by pickers from the picking zone to the delivery zone. The conveyor belt divides cross-aisles 2 and 4 into two sub-cross-aisles. Let us use the words up and down to distinguish between these two sub-cross-aisles. Pickers walk around the warehouse pushing a trolley that contains up to f boxes. Due to the particular warehouse layout, the pickers cannot cross the cross-aisle 3. As a consequence, picking tours are limited to blocks 1 and 2 or to blocks 3 and 4. Moreover, all the pickers follow the same circular path orientation and always walk following a clockwise walkway. This is due to space limitations in the aisles, and also facilitates picking operations since all routes are in the same direction. Thus, in a given route, the locations names follow a topological order. We call this particular policy the walking policy. It is formally defined in Section 3. In this paper, we focus on the optimization of the operations in the picking zone, namely the work of the pickers.

Let us illustrate the picking tours with an example. We consider blocks 1 and 2 and number the aisles as follows. Front-aisle and back-aisle numbered 1 are the bottommost of block 2. Other aisles are numbered clockwise. Enumeration of blocks 3 and 4 follows the same rule. Only enumeration of the front-aisles is reported in Figure 1 due to space limitation. A possible picking tour (solid green lines in Figure 1) starts at front aisle 5 and continues to front aisles 6 and 7. The picker then puts the boxes prepared on the conveyor belt and starts another picking tour. The tour may start from the front aisle 7, or the picker may walk back through cross-aisle up 2 and cross-aisle down 2 and starts a new picking tour.

Each working day a set of demands has to be prepared. The demands are associated with customers, i.e. shops located throughout France and abroad. Based on the sales, the shop demands are computed during the night. Demands usually consist of several items which may not be put into a single box. The demand of each customer has then to be split in different boxes. Clearly one box must contain only products associated with a single customer.

The problem under study consists in optimizing the productivity in the picking zone by simultaneously determining: (1) the products that constitute each box to be delivered to customers; and (2) the batch of boxes to constitute a picking tour.

Literature Review

Warehouse management involves different optimization problems such as conceiving picking tours, batching orders, storage assignment, layout design, and zoning (see, for example, de [START_REF] De Koster | Design and control of warehouse order picking: A literature review[END_REF]). Decisions to optimize productivity in a warehouse can have different horizons and be classified as strategic, tactical and operational. Strategic decisions determine the layout of the warehouse, the disposition of each zone (receiving, storage, picking and delivery) with respect to the others, the decisions that involve the storage policy as well as picking policy (automated versus human). Tactical decisions can involve the location of products based on their forecast demand in the storage and picking zone. At the operational level, storage and picking tours need to be efficiently determined and coordinated. Storage tours move products from the storage to the picking zone, while picking tours aim at collecting items to satisfy the demand of customers [START_REF] De Koster | Design and control of warehouse order picking: A literature review[END_REF], [START_REF] Marchet | Investigating order picking system adoption: a casestudy-based approach[END_REF]).

Picking operations can be classified based on the automatism introduced in the process. As in [START_REF] Marchet | Investigating order picking system adoption: a casestudy-based approach[END_REF] five classes can be identified: picker-to-parts, parts-to-picker, pick-to-box, pick-and-sort, and completely automated picking. In the first category, pickers move around the warehouse to pick items while in the parts-to-picker class, automated devices bring loads to pickers that are in charge of picking the right quantity required by the order under consideration. In pick-to-box systems, pickers are assigned to different zones and boxes containing customer orders move using a conveyor through different zones to be filled. In pick-to-sort systems, orders are batched and pickers collect items of a certain product to satisfy orders in the whole batch. A conveyor brings all the items to the sorting area where items are assigned to orders. Finally, humans do not intervene in full automated picking systems.

Picking systems where the human presence is necessary for operations still involve the majority of the warehouses (de Koster et al. ( 2007)). Recent studies (Ecommerce Europe (2017), [START_REF] Marchet | Investigating order picking system adoption: a casestudy-based approach[END_REF]) show how the cost of the automation would be too high to be profitable in a short or mid-term horizons. For these reasons and since in the problem under study picking operations can be classified as picker-to-parts, we concentrate the review on papers related to this area.

In the HappyChic case, our interest is focused on the operational decisions for picking operations. It concerns order batching and picker routing problems. Usually, order batching is performed when the size of orders is smaller than the capacity of the trolley used to collect the items. Order batching consists in grouping a set of orders into picking tours to minimize the total distance traveled by pickers. This problem is NP-hard [START_REF] Gademann | Order batching to minimize total travel time in a parallelaisle warehouse[END_REF]). Given a set of locations to visit, the picker routing problem aims to decide the sequence of locations to visit to minimize the traveled distance. These two problems have been studied extensively as separate problems [START_REF] De Koster | Design and control of warehouse order picking: A literature review[END_REF]). Some recent works have addressed the joint order batching and picker routing problem. [START_REF] Valle | Optimally solving the joint order batching and picker routing problem[END_REF] propose a branch-and-cut algorithm based on a mixed integer programming formulation, and [START_REF] Scholz | Order batching and picker routing in manual order picking systems: the benefits of integrated routing[END_REF] proposed an iterated local search approach.

In the HappyChic case, we face a specific batching problem. Indeed, the size of an order is usually larger than the capacity of a trolley. Hence, the batching problem has two levels. First, items of a given order have to be batched into boxes. Second, boxes have to be assigned to trolleys. It is noteworthy that these two batching levels are mutually dependent.

Problem definition

In this section, we formally introduce the problem and the notation used in the paper. All along the following sections, we suppose that each working day is divided into T periods. We indicate by T = {1, . . . , T } the period set. each product p ∈ P is associated the volume v p and the weight w p of a single item.

Picking zone

The picking zone of the warehouse consists of four blocks each containing different aisles having racks on one or both sides. Racks are divided in several locations that contain cartons. Locations are filled from the back and emptied from the front. This defines what we call front-aisles and back-aisles (see Figure 2). Only front-aisles are involved in the picking operations.

An order is imposed on the front-aisles of the picking zone. First, aisles of blocks 1 and 2 and aisles of blocks 3 and 4 are ordered starting from the bottommost aisle of blocks 2 and 4 to the bottommost aisle of blocks 1 and 3 following the clockwise sense. Let A be the set of front-aisles. Given two front-aisles a 1 , a 2 ∈ A, aisle a 1 precedes aisle a 2 and we write a 1 ≺ a 2 , if a 1 comes before a 2 in this given order.

The set of all possible locations in the aisles of the picking zone is indicated by L. A particular order is imposed on locations. Let us define the function a : L → A that associates each location with its aisle. Given two locations l 1 and l 2 such that a(l 1 ) = a(l 2 ), we say that l 1 precedes l 2 , and we indicate it by l 1 ≺ l 2 , if a(l 1 ) ≺ a(l 2 ). On the other hand, if a(l 1 ) = a(l 2 ), we define the following rule. Entering a front-aisle from the cross-aisle, locations on the left are ordered from the beginning to the end of the front-aisle, while locations on the right are ordered on the opposite way. Thus, when a(l 1 ) = a(l 2 ), we write l 1 ≺ l 2 if:

• l 1 is on the left-side of the aisle, while l 2 on the right;

• l 1 and l 2 are both on the left-side and l 1 is nearer to the beginning of the aisle than l 2 ;

• l 1 and l 2 are both on the right-side and l 1 is nearer to the end of the aisle than l 2 . This order is called walking order or walking policy. An example is given in Figure 2, where the numbers on each location follow the walking order. Given two locations l 1 and l 2 , distances d l 1 ,l 2 and d l 2 ,l 1 are determined. They represent the distance of the path that goes from one location to the other following the walking order. Due to logistic constraints in the warehouse, it can be impossible to reach location l 2 from location l 1 (and/or vice versa). In this case, the distance is infinite. Location 0 represents a location where a tour starts and ends. Picking tours start and end at any place along the conveyor belt. This can be approximated by the entry point or exit point of an aisle. Hence location 0 is a fictive point that represents entrance or exit of an aisle. We then consider distance d 0,l as the distance from the entrance point of aisle a(l) to location l, and distance d l,0 as the distance from location l to the exit point of aisle a(l).

At a given period t ∈ T , each location contains items of the same product. However, different locations in the picking zone may contain the same product. For each period t, we assume that the storage level in the picking zone is sufficient to satisfy the total demand.

Pickers

A team of employees is in charge of collecting items from the cartons in the picking zone and to fill the boxes that will be delivered to shops. These employees are called pickers. It is supposed that the number of pickers is large enough to ensure the preparation of the total demand to be prepared within a given period.

A picker walks around the warehouse respecting the walking policy. The picker pushes a trolley with up to f boxes available. He/she picks products from cartons and put them into the boxes. Once the all required items are picked, he/she walks to the conveyor belt on which he/she places the boxes. He/she then starts another picking tour with a trolley containing empty boxes. Then, the picker is allowed to walk against the walking policy to locate him/herself at the beginning of the next tour. Note that, due to the particular walking policy considered in the warehouse, the set of locations to visit to pick the items directly implies the path that a picker needs to follow to collect the items.

Demands

A set of demands D is to be prepared during a working day. The set D is partitioned in T subsets/periods, in which the working day is divided. Subset D t , t ∈ T , contains all the demands that need to be prepared during period t.

For each demand d ∈ D, a set of products P d = {p d 1 , . . . , p d n d } ⊆ P needs to be collected. A quantity q d j and a location l d j are associated with each product p d j in demand d. This corresponds to items of the same product, situated at the same location, and needed to prepare d. Note that different items of the same product could be picked from cartons located in different locations in the picking zone. We assume that given product p d j in demand d, its location l d j is set. This assumption is satisfied for the HappyChic case since they use a FIFO rule to empty the different locations associated with the same product.

Whenever useful for the presentation of this work, a set notation is used as follows. Given a demand d, p ∈ d, p ∈ P means that product p is required by demand d. Moreover, if an index is not necessary, it is removed to ease reading. For instance, q p d indicates the quantity of product p required in demand d.

Boxes

Boxes to be delivered to shops must respect a maximal volume v max box and a maximal weight w box . Moreover, each box must satisfy a minimal volume v min box of products. A box, that is not filled at least as much as the minimal volume, can be crushed by other boxes during transportation. Filling up a box increases its strength. Note that it may happen that a demand is too small to ensure the minimal volume capacity constraint. It is then the decision of the managers of the warehouse to decide if this demand is prepared (with violation of the capacity constraint) or delayed to another day. Here, we minimize the number of boxes that violated the minimal volume constraint and consider that they are integrated into the proposed picking tours.

Objective and constraints

The problem consists in simultaneously determining for each period t ∈ T : (1) the partition of each demand d ∈ D t into boxes; (2) the batching of boxes into groups of up to f boxes, i.e., the determination of a picking tour in order to minimize the number of boxes that violate the minimal volume constraint and the distance covered by the pickers while satisfying the maximal volume and weight constraints, the walking policy imposed into the picking zone, and the trolley capacity.

Note that the objectives are prioritized. First, we minimize the number of boxes that violate the minimal volume constraint, and second, the walking distance. Moreover, even if the objective is to minimize the total distance walked by the pickers, the industrial partner wants to ensure the number of tours remains reasonable. Indeed, since there is not a unique starting location (a tour can start at any place along the conveyor belt), minimizing the distance could lead to generate a large set of small tours. This is not acceptable since starting a new tour requires some time for the picker, that can be viewed as a fixed cost. Hence, we impose a fixed distance penalty on all tours.

Algorithmic approach

In this section, we describe the algorithm that has been developed to tackle the problem. Our algorithm includes two phases. The first phase splits each demand into different boxes (Section 4.1) while the second phase generates picking tours once the boxes have been determined (Section 4.2).

Splitting demands

Given a particular demand d, it may happen that the items cannot be contained into a single box or it may be advantageous to consider several boxes to increase the pickers' productivity. In these cases, the demand needs to be divided into several boxes able to contain it.

To perform this task, we apply the following algorithm inspired from the procedure split proposed by [START_REF] Beasley | Route first-cluster second for vehicle routing methods[END_REF] and successfully used by [START_REF] Prins | A simple and effective evolutionary algorithm for the vehicle routing problem[END_REF] in the context of the vehicle routing problem. This procedure is based on the computation of shortest paths on acyclic graphs. The algorithm we propose, seeks for a distribution of a given demand in several boxes such that volume and weight constraints are respected while increasing productivity, i.e. yield to operational plans with the shortest possible walking distance.

Let us suppose that products {p 1 , p 2 , . . . , p n } are required to satisfy demand d. Let us associate a quantity q i with each product p i included in demand d. It follows that demand d requires a total of I = n i=1 q i items. Let us assume that each item i ∈ I, I = {1, . . . , I = n i=1 q i } composing demand d has to be collected from the location l i of the picking zone and let us suppose, without loss of generality, that items are ordered following the walking policy, namely such that i, j ∈ I, i < j implies l i ≺ l j (see definition in Section 3.1).

The total volume of demand It means that the items of such demands may be split into several boxes. We construct a graph G = ( Ī, E) where Ī is the node set while E is the arc set. Set Ī contains 1 + I nodes, numbered 0, 1, . . . , I = n i=1 q i . Node 0 is a dummy node, while node i ∈ Ī, i > 0 represents the i-th item in I forming demand d.

d is v d = n i=1 q i v p i ,
With each node i ∈ Ī, we associate the quantities v i (resp. w i ) the volume (resp. the weight) of the item represented by node i.

Then, for each couple of nodes i, j ∈ Ī we calculate three quantities:

• d ij = d 0,l i+1 + j h=i+2 d l h-1 ,l h + d l j ,0 , i.e.
, the walking distance needed to pick items i + 1, . . . , j from their locations l i+1 , . . . , l j ;

• w ij = j h=i+1 w h , i.e., the weight of items i + 1, . . . , j;

• v ij = j h=i+1 v h , i.e., the volume of items i + 1, . . . , j.

Based on these quantities, we can determine the arcs of graph G in such a way that each arc represents a feasible box. We first determine the arc set E 1 defined as:

E 1 = {(i, j)|i < j, i, j ∈ Ī, v min box ≤ v ij ≤ v max box , w ij ≤ w box }. (1)
Second, in order to guarantee the graph connectivity, we consider the arc set E 2 made of all arcs (i -1, i), i ∈ Ī \ {0} even if the associated box does not respect the constraint on the minimal volume. Thus,

E 2 = {(i -1, i)|i ∈ Ī \ {0}, v i-1,i < v min box }. ( 2
)
We finally define the arc set

E = E 1 ∪ E 2 .
The cost c ij associated with each arc (i, j) ∈ E 1 is given by the quantity α BOX + d ij , where α BOX is a parameter that penalizes the creation of new boxes. This parameter α BOX permits to limit the number of boxes, hence to propose a reasonable number of tours. The cost c ij associated with each arc (i, j) ∈ E 2 is given by δ BOX + d ij , where δ BOX is a large-enough value.

To determine the distribution of the demand d into boxes, we compute the shortest path from the dummy source node 0 to node I with respect to costs c ij . Each arc belonging to the shortest path represents a box. The shortest path provides the best split of products into boxes given the order of products and the cost structure.

The value of δ BOX should be such that an arc in E 2 is selected only if needed to guarantee the existence of a path between nodes 0 and I. Whenever the shortest path contains an arc in E 2 the corresponding box does not respect all the constraints.

The algorithm is repeated for each demand d ∈ D to obtain its distribution into boxes. All the boxes that are created are inserted in set B.

An example of the split procedure is given in Figures 3a and3b. Here the demand is composed of six different products noted a, b, c, d, e, f . They are located in the warehouse as illustrated in Figure 3a. Values on the ordinate axis represent the distance from the beginning of the aisle to the product location, while values on the abscissa represent intra-aisle and inter-aisle distances. Following this disposition, we have that a ≺ b ≺ c ≺ d ≺ e ≺ f . We suppose that only one item per product is required. For simplicity, let us suppose that all the products have a volume and a weight equal to 10 units, the minimal volume for a box equals 10, while the maximal weight and volume equal 30 units.

The graph that is built with respect to this demand is reported in Figure 3b. A box that contains only product a is associated with a distance of 50 + α BOX units. The picker needs to get into the aisle and walk for 20 units to reach product a. He/she then needs to cross and to get out from the aisle for an additional distance of 30 units. This is represented by arc (0, a). A box that contains products b and c is associated with a distance of 90 + α BOX , that is the distance to get into the first aisle, walk until the position in front of b, cross the aisle to reach b, get out from the aisle, and walk to the next aisle for a total of 50 units. Then, the picker needs to get into the second aisle and walk until c. He/she then crosses the aisle and gets out of it for an additional 40 units of distance. Thus the total distance to compose the box is 90 + α BOX . This is represented by arc (a, c).

Note that if α BOX ≤ 50 the shortest path that goes from node 0 to node f is made by arcs (0, b), (b, e) and (e, f ) for a total cost of 220+3α BOX . The demand is then split in three boxes. One containing products a and b, another with products c, d and e, the last with only product f . On the other case, when α BOX > 50, the best solution is to create only 2 boxes. The shortest path in this case is made by arcs (0, c) and (c, f ) for a total cost of 270 + 2α BOX . The solution in this case considers two boxes, one with products a, b and c, the other with products d, e and f .

Batching boxes for picker trips

A picker pushes a trolley containing up to f boxes. We thus need to batch the boxes in order to obtain the trips that pickers have to perform. Let us suppose that a set B = {1, . . . B} of boxes have to be prepared. Let a S b be the starting aisle associated with box b ∈ B, namely the aisle where is located the first item to pick to prepare box b (with respect to the walking policy). Let a T b be the ending aisle associated with box b ∈ B: the aisle where is located the last item to pick to prepare box b. Moreover, let d b be the walking distance to pick all items of box b. We suppose, without loss of generality, that boxes are numbered such that b

1 < b 2 , b 1 , b 2 ∈ B if and only if a S b 1 = a S b 2 and a S b 1 ≺ a S b 2 or a S b 1 = a S b 2 and a T b 1 ≺ a T b 2 or a S b 1 = a S b 2 and a T b 1 = a T b 2 and d b 1 > d b 2 .
To create batches of boxes of cardinality up to f , we adapt the split procedure described in Section 4.1. Similarly, a graph G = ( B, E) is constructed. The node set B contains 1 + B nodes, numbered from 0 to B. Node 0 is a dummy node, while node i = 1, . . . , B represents the i-th box in B. The arc set E contains all the possible batches of boxes, namely, all batches with up to f boxes. Thus

E = {(i, j)|i, j ∈ B, j -i ≤ f }
A cost c ij is associated with each arc (i, j) ∈ E and represents the total walking distance needed to collect all items in boxes i + 1, i + 2, . . . , j, plus a fixed penalty α T OU R to start a tour. Note that c ij can be computed keeping in memory all the aisles that contain items in the boxes under consideration and the location of the furthest item in each aisle.

The shortest path that goes from the dummy node 0 to node B provide a partition of the box set B into batches of up to f boxes that minimizes the walking distance covered by the pickers.

Homogeneous boxes

The boxes contain clothes that are exhibited in shops. When a box arrives at the shop, a salesperson carries it around the shop to prepare the displays. His/her job is made easier when boxes are homogeneous. Homogeneous means that a box contains clothes of a single class. Classes can be defined as needed, but an example of classes can be t-shirts, jeans, shirts, pulls, etc. For instance, a box containing 10 t-shirts is homogeneous, while a box containing one jeans and nine t-shirts is not.

Whenever homogeneity of the boxes is prior in the composition criteria, we propose to modify the split algorithm proposed in Section 4.1 to take this feature into account. The arc set defined by Equation (1) does not change. However, we modify the cost associated with each arc by increasing it by a factor β whenever it considers a non-homogeneous box. Thus, the cost of an arc (i, j) such that items i + 1, . . . , j do not belong to the same class is α BOX + β + d ij . When the value of β is increased, more homogeneous but also more numerous boxes are generated (especially when β and α BOX are similar).

Moreover, if the minimal number of boxes required to contain the demand of a given customer, namely,

n min box = max i∈I v i v max box , i∈I w i w max box
is greater than a certain parameter γ, the order on which items are considered to build the acyclic graph is modified and does not depend anymore on the walking policy. In particular, items are first ranked such that items of the same class follow each other, and inside each class, items are ranked according to the walking policy. This increases the chance of having homogeneous boxes with the disadvantage of deteriorating the total walked distance.

Computational results

Our industrial partner provided us with a set of 487 instances. period, one instance is created for the demands to be prepared in blocks 1 and 2, and another instance is created for the demands to be prepared in blocks 3 and 4. For each instance, the set of demands corresponds to a number of items to be picked that varies between 45 and 18 209, with an average of 2 469. We run the algorithms presented in Sections 4.1-4.2 on the set of 487 instances to evaluate its efficiency. The parameters are set to keep a reasonable number of boxes and tours in the solution. Their value are: α BOX = 30 000 and α T OU R = 15 000. Each instance is solved on a Intel(R) Core(TM) i7-4600U CPU@2.10GHz, in less than one second. The overall results are reported in Table 1. From Table 1 we can notice that the number of boxes and the number of picking tours are similar in the solutions computed and in the solutions performed by HappyChic. However, the number of boxes with a volume lower than the minimal volume is drastically reduced in the proposed solution. Moreover, our approach drastically reduces the distance walked by pickers and, consequently increases their productivity. When we consider all the 487 instances individually, the proposed solution reduces the total distance in 462 instances. The deviation from the current industrial solution ranges from -90% to +8%, with an average of -21.85%. Hence the proposed solutions lead to significant improvements for the industrial partner.

As explained in Section 4.3, the tasks of salespersons in shops are easier when boxes are homogeneous. In Figures 4a and4b we present results obtained by the modified split algorithm proposed in Section 4.3 when homogeneity is explicitly taken into account. The blue dot in the left part of the figures corresponds to the industrial solution which is the worst among the others in term of homogeneity. The solution obtained with the algorithm presented in Section 4.1 without taken homogeneity into account is represented by the red square in the left part of the figures. This solution already improves the homogeneity of the boxes with respect to the industrial solution.

We have tested several values for parameters β and γ introduced in Section 4.3. β ranges from 5 000 to 40 000 with a step of 5 000, and γ ranges from 1 to 10 with a step of 1. Each combination of β and γ corresponds to a point in Figures 4a and4b.

The homogeneity of the boxes can be improved by applying the changes described in Section 4.3 at the cost of degrading the total walking distance and the number of boxes. It is possible to reach 35% to 45% of homogeneous boxes with a small increase in the distance and the number of boxes. This is achieved with β ≤ 20 000 and γ ≥ 6. Homogeneity can reach up to 82%, at an expensive cost in the number of boxes. However, with 82% of homogeneous boxes, the distance is still lower than the distance covered in the current industrial solution.

Conclusions

In this paper, we presented with an integrated order batching and picker routing problem based on an industrial case encountered in a warehouse of the company HappyChic. The specific features of the studied problem are: (1) orders have to be split into several boxes, (2) boxes are batched into trolleys, and (3) the picking zone has a specific configuration, so the optimal routing policy is easy to determine when the locations to visit are known. To deal with large size instances, we proposed a heuristic algorithm based on a double split procedure: the first one to create boxes, the second one to batch boxes into picking tours. The proposed algorithm can easily be adapted to take into account several classes of products with a homogeneity criterion when creating the boxes. The proposed results are very positive for the company: over a set of 487 instances, the total picking distance can be decreased by 24.21%.

ForkliftFigure 1 :

 1 Figure 1: The warehouse.

Figure 2 :

 2 Figure 2: Front-aisles and back-aisles. Legend as in Figure 1. Locations are numbered following the walking order.

  and the total weight of demand d is w d = n i=1 q i w p i . If v d < v min box and w d ≤ w max box , all items of demand d are put into a single box that violates the minimal volume constraint. In the following, we only consider demand d such that v d ≥ v min box or w d > w max box .

Figure 3 :

 3 Figure 3: Example of the split procedure.

  boxes versus total distance.

Figure 4 :

 4 Figure 4: Solutions taking into account homogeneity: each point is a solution.

  The warehouse considered contains a set of different products P. With

	12	front-aisle	13	36	front-aisle	37	60	front-aisle	61
	11 10	a 1	14 15	35 34	a 2	38 39	59 58	a 3	62 63
	9		16	33		40	57		64
	8		17	32		41	56		65
	7 6		18 19	31 30		42 43	55 54		66 67
	5 4 3 2 1		20 21 22 23 24	29 28 27 26 25		44 45 46 47 48	53 52 51 50 49		68 69 70 71 72
			1

Table 1 :

 1 These instances come from 17 working days between January and November 2017. Working days are divided into 11 to 16 periods. For each Global comparison with the industrial solutions.

		Industrial Proposed	Deviation
		solution	solution in percentage (%)
	Nb. of boxes	47 201	47 354	+0.32
	Nb. of picking tours	8 105	8 103	-0.02
	Nb. box violate volume min	575	33	-94.26
	Total walking distance (m)	1 298 958 984 479	-24.21
	Nb. homogeneous boxes	12 463	13 978	+12.16
	Pct. of homogeneous boxes	26.40%	29.52%	-

http://www.happychicgroup.com


M Agnaou

D Lasseux
email: didier.lasseux@ensam.eu

A Ahmadi

On the Understanding of the Macroscopic Inertial Effects in Porous Media Investigation of the Microscopic Flow Structure

Inertial flow in porous media occurs in many situations such as flow in column reactors, in filters or near wells for hydrocarbon recovery. It is characterized by a deviation from Darcy's law. Darcy's law <v β > ∝ ∇<p β > β Darcy's law + inertial correction non-linear relation between <v β > & ∇<p β > β correction ∝ <v β > 2 correction ∝ <v β > 3 Darcy's regime weak inertia regime strong inertia (Forechheimer) regime regime above strong inertia creeping flow laminar stationary flow laminar unsteady flow turbulent flow Re d time averaged macroscopic models correction ∝ ?

The existence and origin of different regimes, function of the microstructure and flow orientation, are still not well understood. We provide an in depth analysis of the flow structure to identify the origin of the deviation from Darcy's law using a theoretical justification. The role of the recirculation zones is highlighted.

INTRODUCTION Main Objectives

1. Identify different inertial regimes with respect to:

• the shape of the solid inclusions,

• the structural disorder,

• the orientation of ∇ p * β β .

2. Identify the origin of the deviation from Darcy's law with respect to:

• the role of the recirculation zones,

• the flow structure to microscopic convective acceleration relationship.

Macroscopic Model & Determination of Macroscopic Coefficients

The macroscopic model was obtained from the volume averaging of the Navier-Stokes equations [START_REF] Whitaker | The Forchheimer equation: a theoretical development[END_REF] v K * , H * and F are determined from the solution of the closure problem on M * and m * that are periodic.

* β = -K * .∇ p * β β -F. v * β = -H * .∇ p * β β (1a) ∇. v * β = 0 (
Re * ∇ p * β β .M * T .∇M * = ∇m * -∇ 2 M * -I (2a) ∇.M * = 0 (2b) M * = 0 at A βσ (2c) M * = H * = H/l 2 (2d) Re * = ρ β l 3 µ 2 β ∇ p β β : Reynolds number (define the forcing intensity). ∇ p * β β
: unit vector (define the macroscopic flow orientation θ). As in [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF], the deviation from Darcy's law was analyzed in terms of the inertial correction vector

f c = -F. v * β / v * β (3)
and with respect to the Reynolds number

Re k = v * β √ k * Re * = √ k * d * Re d (4)

RESULTS

Inertial Correction

T ≈ 1.04 T ≈ 1.2 T 1.16 T ≈ 1.13 T ≈ 1.12 Re k f cx OS θ = 0 • OC θ = 0 • OS θ = 45 • OC θ = 45 • DS θ = 0 • DC θ = 0 • • Intensity of inertial effects de- pend on the microstructure & orientation of ∇ p * β β
as observed before [START_REF] Lasseux | On the stationary macroscopic inertial effects for one phase flow in ordered and disordered porous media[END_REF].

• Inertial effects increase with tortuosity T . • Contribution of structural disorder to inertial effects is larger than that due to the solid inclusions shape.

(@f cx /@Re k ) n OS ✓ = 0 OC ✓ = 0 OS ✓ = 45 OC ✓ = 45 DS ✓ = 0 DC ✓ = 0
weak inertia regime strong inertia regime

• Weak inertia regime (f cx ∝ Re 2
), is identified in all configurations over almost the same Re k interval except on the disordered structures where it is shorter.

• Strong inertia (f cx ∝ Re k )
depends on the microstructure and orientation of ∇ p * β β but is more robust in the presence of structural disorder.

• A regime above strong inertia is observed over large Re k intervals.

Role of the Recirculation Zones

θ = 0 • θ = 0 • OS θ = 45 • OC θ = 45 • Re k V βV ort /V β Weak
θ = 0 • θ = 0 • OS θ = 45 • OC θ = 45 • θ = 0 • DC θ = 0 • Re k T ortuosity
Weak inertia regime Transition Strong inertia regime Regime above strong inertia 1

• Kinetic energy lost in vortices is negligible as it is always < 0.5% of the total kinetic energy.

• Volume of recirculation zones (V βV ort ) is important & its variation with Re k leads to streamlines deformation.

• Enlargement of recirculation zones beyond weak inertia regime decreases tortuosity (with different rates).

• Flow structure in OS. Recirculation zones delimited by red streamlines enlarge with Re k .

Relation between Convective Term & Flow Structure

Decomposition of the convective acceleration in the orthonormal basis (t, n) of the Frenet frame: • Whereas the 1st term of the RHS of Eq. 5 to streamlines curvature (κ) weighted by the local kinetic energy, the 2nd characterizes the variation of kinetic energy along the streamlines. The 1st term can be related to tortuosity, invoked heuristically in the literature. • Constant flow structure on weak inertia regime resembling that of creeping flow is observed.

v * β .∇ v * β = v * 2 β κn + 1 2 dv * 2 β ds t (5) 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 •10 -3 θ = 0 • θ = 0 • OS θ = 45 • OC θ = 45 • θ = 0 • DC θ = 0 • Re k v * 2 β κn β Weak inertia regime Transition Strong inertia regime Regime above strong inertia 1 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 0 1 2 3 4 5 6 7 •10 -4 θ = 0 • θ = 0 • OS θ = 45 • OC θ = 45 • θ = 0 • DC θ = 0 • Re k 1 2 dv *
• Beyond weak inertia, different variations occur explaining the various behaviors of the inertial correction.

CONCLUSIONS • As in [2],

weak inertia regime is always well identified, strong inertia regime is robust in the presence of structural disorder yielding a shorter weak inertia regime and stronger inertial effects. • The kinetic energy lost in vortices is insignificant.

• Streamlines shape induced by the recirculation zones implies variations in the flow structure which is itself correlated to the flow regimes. • The intensity of inertial effects and the occurrence of different inertial regimes is correlated to the flow structure based on a formal justification.

-The higher the tortuosity, the stronger inertial effects. -The weak inertia regime is characterized by a constant flow structure. -Beyond weak inertia, different and complex flow structure variations occur explaining the weakness of the Forchheimer model.

More details about this work can be found in [START_REF] Agnaou | The origin of the deviation from darcy's law: a further analysis of the microscopic flow structure[END_REF].

•

  Ordered (a and b) / disordered (c and d) model structures of parallel cylinders. Porosity = 75%. • Representative Elementary Volume (REV) made of 1 × 1 (a and b) / 30 × 30 (c and d) geometrical unit cells.

  1b) v * β : filtration velocity. ∇ p * β β : macroscopic pressure gradient. K * = k * I = K/l 2 : intrinsic permeability tensor (structure dependent only). F: inertial correction tensor. H * = (I + F) -1 K * = H/l 2 : apparent permeability tensor. H * and F depend on both the structure and the flow.


V Daru

M.-C Duluc

O Le Maître

P Le Quéré

NUMERICAL SIMULATION OF LOW MACH NUMBER LIQUID GAS FLOWS

This work is devoted to the numerical simulation of liquidgas flows. The liquid phase is considered as incompressible, while the gas phase is treated as compressible in the low Mach number approach. We present a model and a numerical method aimed at the computation of such two-phase flows. The numerical model uses a lagrangian front-tracking method to deal with the interface. To validate our model, we present preliminary results in the 1D case of a heated enclosure, inside which a liquid zone is inserted between two gaseous zones.

INTRODUCTION

This work is devoted to the numerical simulation of liquidgas flows with realistic physical properties and high density ratios between the phases. For confined flows (as is often the case in nuclear engineering or microfluidic applications), the variations of the thermodynamic pressure in the gaseous phase play an important role and must be taken into account. This means that the gaseous phase must be treated as a compressible fluid. However, in most of the applications the velocities in the gas are very small, and a low Mach number approach is well adapted for such flows. The liquid being considered as incompressible, an adequate model must be able to take into account both a truly incompressible flow and a low Mach number flow, separated by a dynamic interface across which the physical properties of the fluid are discontinuous.

We present such a model together with a numerical method aimed at the effective and efficient computation of such two-phase flows. The model is based on the solution of the Navier-Stokes equations, and it is developed by using a one-field approach, where the liquid or the gas is indicated by a Heaviside function. The low Mach number model in the gas involves two pressures, the usual dynamic pressure and a thermodynamic pressure which is defined only in the gaseous regions and is obtained using integral relations. The numerical model uses a lagrangian front-tracking method to deal with the interface. This method has shown to be effective for complex multidimensional situations including interface merging/breakup and phase change (Tryggvason et al. (2001) or [START_REF] Shin | Modeling three-dimensional multiphase ow using a level contour reconstruction method for front-tracking without connectivity[END_REF]).

To validate our model, we present an application in a one-dimensional case where the effects of the thermodynamic pressure are significant. Comparison with a reference solution shows that the front tracking model gives accurate solutions, and is adequate for treating multidimensional application involving compressibility effects.

PHYSICAL MODEL

We consider a two-phase liquid-gas flow, the liquid phase being incompressible and the gaseous phase being compressible. The Mach number in the compressible phase is supposed to be always small. The flow can then be represented using a low Mach model as proposed by [START_REF] Chenoweth | Natural convection in an enclosed vertical air layer with large horizontal temperature differences[END_REF]. In this model, the pressure is decomposed in two components, a thermodynamic pressure which is uniform in space, and the usual hydrodynamic pressure which is involved in incompressible flows. Let us introduce a Heaviside function H(x,t), which plays the role of the characteristic function of the gaseous phase (H is equal to 1 in the gas, and 0 in the liquid phase). The function H is simply advected with the fluid velocity, and follows a transport equation :

0 = ∇ ⋅ + ∂ ∂ = H v t H dt dH r
The two-phase single-field model can be written, expressing the conservation of mass, momentum and energy:

( ) 0 1 1 1 1 ( ) p p v t v v v t T d v T k T H t c c ρ ρ π τ ρ ρ ρ ρ ⎧ ∂ ⎪ + ∇ ⋅ = ∂ ⎪ ⎪ ∂ + ⋅∇ + ∇ = ∇⋅ ⎨ ∂ ⎪ ⎪ ∂ + ⋅∇ = ∇⋅ ∇ + ⎪ ∂ ⎩ r r r r r P d t (1)
where P(t) is the thermodynamic pressure, π is the hydrodynamic pressure and τ is the viscous tensor. The effect of gravity and surface tension is not included here. The gas is assumed to follow a perfect gas law:

g P r T ρ = (2)
with r = c p -c v . The liquid being considered as incompressible has a constant density ρ l . For our twophase model, the generalized equation of state then writes:

(1 ) (1 )

g l l P H H H H rT ρ ρ ρ = + - = + -ρ (3)
The thermodynamic pressure must be calculated by using integral relations that will be derived in the following. In equations (1), the coefficients c p and k are discontinuous across the interface, and are obtained in the same way as in (3) from their values in the gas and in the liquid.

Using equation (3), the temperature equation in (1) can also be written in the following form, using the grouping (1/P)(dP/dt) rather than dP/dt :

1 1 ( ) p dT dP k T T H dt c P dt 1 γ ρ γ - = ∇⋅ ∇ + (4)
Let us now derive the integral relations for P(t). To this aim, we insert (3) in the mass equation in (1) obtaining:

1 1 ( ) dT dP v H T dt P dt ⎛ ⎞ ∇ ⋅ = - ⎜ ⎝ ⎠ r ⎟ (5)
Let us suppose that the gaseous phase in enclosed by walls or the liquid phase (case of bubbles for example), and let us denote by Ω g (t) its spatial extension, Σ g (t) the bounding surface. We then integrate (5) over Ω g (t) to obtain:

( ) ( ) ( ) 1 1 1 g g g t t t dP dT H dx v nd P dt T dt Hdx Ω Σ Ω ⎛ ⎞ = - ⎜ ⎝ ∫ ∫ ∫ r r S ⋅ ⎟ ⎠ (6)
In this way we calculate (1/P)(dP/dt) to obtain the source term in (4). Once (1/P)(dP/dt) is known, the pressure P(t) can be obtained directly from a time integration, i.e.:

δ + ⎛ ⎞ + = ⋅ ⎜ ⎝ ⎠ ∫ ' ⎟ (7)

EFFECT OF THE THERMODYNAMIC PRESSURE ON THE DYNAMICS OF THE FLOW

The model described above does not take into account the effect of thermodynamic pressure that can be encountered for example if several bubbles are enclosed in a liquid. In this case, each bubble has its own pressure, and the net effect of these pressure forces result in an acceleration of the liquid. Let us illustrate this on a simple 1D model case: a gaseous zone enclosed inside two liquid moving zones. The incompressible motion of the liquid in 1D is a uniform velocity, named V 1 (t) and V 2 (t) respectively. If we assume that the gas density ρ is uniform in space (isothermal case), the formulation of the problem in the gas reads (neglecting viscous effects):

1 u x u u u t x ρ ρ x π ρ ∂ ⎧ = - ⎪ ∂ ⎪ ⎨ ∂ ∂ ∂ ⎪ + =- ⎪ ∂ ∂ ∂ ⎩ & (8)
where the notation ρ & denotes the time derivative of ρ.

The density can be calculated directly considering the conservation of total mass M 0 in the gaseous enclosure of instantaneous length l(t), that is ρ(t)=M 0 / l(t). An origin being fixed, let us denote by l 1 (t) and l 2 (t) the positions of the left and right liquid-gas interfaces respectively. One can easily show that the local velocity in the gas is given by:

1 2 ( , ) 2 1 x l x l u x t V V l l - - = - (9)
Replacing this expression in the momentum equation, one obtains the value of the pressure gradient as:

1 2 1 2 1 u u l ll l u x t x x l l π ρ ∂ ∂ ∂ - + =- = + ∂ ∂ ∂ && && && l (10)
On the interfaces, this gives:

1 1 1 2 2 1 l V x l x l V x l π ρ ⎧ 2 = = ∂ ⎪ - = ⎨ ∂ = = ⎪ ⎩ && & && & (11)
The velocities being uniform inside the liquid zones, this means that the hydrodynamic pressure gradient must be uniform in the liquid in order to incorporate the thermodynamic pressures exerted by the gaseous zones enclosing the liquid. Thus we must have:

1 1 l P x L π ρ ρ ∂ Δ = ∂ (12)
in each liquid zone of length L, ΔP being the thermodynamic pressure difference exerted on the liquid. This will be taken into account in the numerical model through a source term added in the momentum equations, as will be developed below.

NUMERICAL SOLUTION

The model described above must be coupled to a specific method for tracking the liquid-gas interface. Among the existing methods, we have chosen here to use the lagrangian front tracking method developed in Tryggvason ( 2001) and subsequent papers. This method uses a lagrangian discretization and movement of the interface. A Heaviside function, defined on the Eulerian grid, is built from the lagrangian definition of the interface. This H function is smoothed over 3 to 4 cells for numerical stability reasons. This leads to the existence of a "mixing" zone between the two fluids, where the model has no physical meaning.

The spatial discretization is based on centered finite differences for both the convection and diffusion terms. A staggered mesh is used, where hydrodynamic pressure and temperature are located at the center of the cells, while the components of the velocity are located on the faces. We have used here a first order explicit temporal discretization, but the extension to semi-implicit second order should not be a problem. A classical prediction-projection algorithm is used to compute the velocities. In order to take into account the effect of the thermodynamic pressure, a source term is added in the prediction step of the velocity that allows the hydrodynamic pressure to reproduce the behavior established in the previous section. The predicted velocity is calculated from:

* v r Π ∇ - + ⋅ ∇ + ∇ - = - l n n n n H v v t v v ρ τ ρ δ 1 ) 1 ( 1 . * r r r r (13)
The added source term is the last term in (13), and is the multidimensional generalization of the right hand side of equation ( 12). It is null in the gas, where H=1, and it is constant in liquid zones. The function P is constructed such as to be defined everywhere in the spatial domain. It extends the thermodynamic pressure defined in the gaseous zones in the liquid zones where it is not defined.

To construct it we impose the following conditions: -it should be equal to the thermodynamic pressures in the gaseous zones, -it should be harmonic in the liquid zones, in order to leave the velocity field in the liquid divergence free, -it should verify a Neumann boundary condition on the walls.

P can thus be interpreted as the potential in the liquid of the forces due to the thermodynamic pressures in the gaseous phases. Let us suppose that there exists N closed gaseous zones in the domain, each marked by a characteristic function H j . The global characteristic function H is equal to the sum of the N functions H j . We can express the above conditions through the following equation that must be verified by P:

( * 1 1 n v v t π ρ δ + ⎛ ⎞ ∇ ⋅ ∇ = ∇ ⋅ -∇ ⋅ ⎜ ⎟ ⎝ ⎠ r r ) 1 ( 14
)
in which stands for the desired divergence of the projected velocity, which is .computed using the mass conservation equation ( 5), after the temperature and the integral in (6) have been calculated. Once the hydrodynamic pressure is obtained, the velocity is calculated from:

1 n v + ∇ ⋅ r Π ∇ - - ∇ - = + l n H t t v v ρ δ π ρ δ 1 ) 1 ( 1 * 1 r r (15)
The last term is subtracted from because the thermodynamic pressure effect is already taken into account in the hydrodynamic pressure field, and otherwise would be taken into account twice.

* v r

To summarize, the complete algorithm reads as follows: -Initialization -Begin the time loop -Solve Eq. ( 4) for T n+1 .

-Calculate H n+1 using the front-tracking method -Calculate ( ) , H n+1 and T n+1 in ( 6) and ( 7).

1
-Calculate ρ n+1 using (3) and update the values k, c p .

-Calculate P and the source term in the momentum equations.

-Calculate

1 n v + ∇ ⋅ r
for the Poisson equation using the mass equation ( 5).

-Calculate the predicted velocities from ( 13).

* v r

-Solve ( 14) for the hydrodynamic pressure π (multigrid).

-Project the velocity by ( 15).

-End of time loop.

For closed domains, an important point to emphasize is the compatibility relation that must be verified when processing the integration of (6) to calculate + . In this case, if we integrate (6) over the whole domain, the surface term must be zero. However, in the front-tracking method, we have pointed out the existence of a mixing zone, typically a few cells thick, inside which the velocity is only defined as an interpolation between the liquid velocity, which is divergence free, and the gas velocity which is not divergence free. Thus there is no chance that the global balance of the velocity on the interfaces be zero if the integration for each gaseous volume ends inside this mixing zone. In fact the integration must cross it and end in the liquid, or equivalently take into account all values of H which are not zero. This is a crucial point to get a correct balance of the thermodynamic pressure effects on the liquid zones.

REFERENCE ALE 1D NUMERICAL METHOD

In order to get a reference solution, we also have developed, in the 1D case, a code based on a Arbitrary Lagrangian Eulerian (ALE) method (moving mesh). In this way, the interfaces between liquid and gas are real discontinuities, and there are no errors that could be ascribed to the front tracking method and the existence of a mixing zone. Using an accurate discretization, we can consider the results given by this code as reference results. Let us emphasize that the ALE method is limited to small distortions situations in the multidimensional case, and is not adequate for treating complex multiphase flows, as opposed to the front-tracking method.

We here briefly recall the ALE method. The domains moving with time, linear transformations are used to map the physical domains to a reference mathematical domain. The spatial discretization uses a high order spectral element for each sub-domain, with order p for the temperature, density, velocity and order p-2 for the hydrodynamic pressure. The mapping being time dependent, it yields additional correction terms in the governing equations to account for the mesh velocity.

The time integration of the flow involves diverse time discretizations and numerical methods for the inversion of the resulting discrete operators. The complete numerical methodology is detailed elsewhere (see [START_REF] Duluc | Numerical study of liquid inclusion study of liquid inclusion oscillations inside a closed 1-d microchannel filled with gas[END_REF]), and we simply mention here that an iterative procedure is used at each time step in order to get converged values of the thermodynamic pressure and the temperature, velocity and density fields.

NUMERICAL RESULTS FOR A 1D PROBLEM

We consider a one-dimensional fluid system consisting of two layers of gas (air) enclosing a layer of liquid (water). The system is closed by two walls. The total length of the system is 100 µm, the length of the liquid layer is 10 µm, the latter being initially situated at the center of the system. The initial thermodynamic conditions are P 0 =101325 Pa, T 0 =293.15K. At initial time, the left wall is heated to T w =373.15K. The right wall is insulated. After a transient evolution, a steady state establishes where the initial positions of the liquid-gas interfaces are recovered, due to mass conservation. In the gas, the density at steady state is unchanged, T f = T w everywhere, and P f =P 0 .T f / T 0 =128976.37 Pa following the perfect gas equation of state. The domain is discretized using 100 grid points. The time step is equal to 10 -8 s, due to stability limits of the explicit time discretization. We represent in Fig. 1 trajectories of the two liquid-gas interfaces. The first remark that can be done is that those trajectories are parallel (the liquid has a solid body movement), and that the initial positions are well recovered. This is a very important and demanding case to test for mass conservation for the method. A method we developed previously based on a single pressure approach [START_REF] Daru | Modélisation et simulation numérique du changement de phase liquide-vapeur en cavité[END_REF]) revealed very serious problems with this criterion.

A logarithmic representation of the interface trajectory (Fig. 1b) reveals that in the first moments after the heating of the left wall, an oscillatory movement of the liquid zone is observed. This effect which is entirely due to compressibility of the gas was analyzed in details by [START_REF] Duluc | Numerical study of liquid inclusion study of liquid inclusion oscillations inside a closed 1-d microchannel filled with gas[END_REF]. The corresponding thermodynamic pressure history is shown in natural (Fig. 2a) and logarithmic (Fig. 2b) scales. As the pressures in the two gaseous domains are not equal, the liquid is accelerated. In Fig. 3 is shown the fluid velocity in the whole domain, for two successive times. One can remark that the velocity is constant in the liquid, showing that the liquid behaves as a solid in this 1D case, and that compressibility effects produce a non-uniform velocity in the gas. In the first microseconds, large velocities can be attained. The hydrodynamic pressure at t=1 µs is shown in Fig. 4. The temperature fields for several successive times is shown in Fig. 5 and reveals that after an initial transient where a piston effect can be observed in the left gaseous zone, a purely conductive regime takes place in the left gaseous zone while the temperature is nearly uniform in the right gaseous zone due to the adiabatic wall condition. The relative total mass variation in the cavity never exceeds 0.5% throughout the unsteady evolution, demonstrating good mass conservation properties of the algorithm.

Finally, in Figs. 6 and7 we compare the results given by the front-tracking algorithm and the ALE approach, the latter being considered as the reference. It can be seen on the trajectories and the pressures history, that the results are very close (Figs. 6 and7), although the front-tracking algorithm is only first order in time. This validates the low Mach compressible/incompressible approach in the front-tracking framework.

CONCLUSION

We have constructed a model aimed at the numerical simulation of liquid-gas flows, where compressibility effects in the gas are taken into account using a low Mach approach. The interface is treated using a fronttracking approach, which has shown to be effective in complex multidimensional situations. The model has been validated in the one-dimensional case by comparison with reference results obtained using an ALE approach. Work is currently in progress to develop multidimensional applications and to extend the model to liquid-vapor phase change problems.

  boundary conditions on the walls if the domain is closed. This implies that we have to solve a Helmholtz equation to calculate P. The hydrodynamic pressure is then calculated by solving the following Poisson equation:

γ

  NOMENCLATUREc p specific heat at constant pressure J/kgK c v specific heat at constant volumetric J/kgK k thermal conductivity W/mK H specific heat ratio π hydrodynamic pressure Pa ρ density kg/m 3 τ viscous tensor Pa Subscripts g related to the gaseous phase l related to the liquid phase

Fig. 1a

 1a Fig. 1a Trajectories of interfaces in natural time scale


Juste Raimbault
email: juste.raimbault@polytechnique.edu

Multi-modeling the morphogenesis of transportation networks

This research introduces a multi-modeling approach to the growth of transportation networks. More precisely, we implement and compare several models, based on biological network growth, cost-benefit rules, and gravity potential breakdown. The resulting multi-modeling framework is calibrated on observed topological data for the European road network. We show that different heuristics are complementary to cover the feasible topological space and that all are necessary to approach existing configurations, what suggests the superposition of corresponding processes in territorial systems.

The growth of transportation networks in territorial systems bears stunning similarities with biological networks, and the understanding of processes driving their spatial extension can both have practical planning applications but also bring theoretical insights into complex morphogenetic systems. Network growth models have been proposed in several disciplines, of which [START_REF] Xie | Modeling the growth of transportation networks: A comprehensive review[END_REF] give a large overview, taking into account diverse processes such as economical, geometrical, geographical processes for example. There exists to our knowledge no systematic comparison of different network generation heuristics. We develop here such a comparison in a multi-modeling paradigm.

Modeling network growth

We introduce a general model of road network growth, conditioned to a fixed population density on a grid, in a sequential way with the following steps. A common core ensures the positioning of new centers preferentially to population density and a direct connection to the existing network. Nodes and links are then added following different potential processes corresponding to the selected heuristic: (i) no supplementary growth (baseline); (ii) random links; (iii) biological network growth following the model of [START_REF] Tero | Rules for biologically inspired adaptive network design[END_REF], where new links are taken as the links with the highest capacity in the stationary flow network obtained by evolving capacities of the full network of potential new links; (iv) cost-benefit network growth as proposed by [START_REF] Louf | Emergence of hierarchy in cost-driven growth of spatial networks[END_REF]; (v) stochastic gravity potential breakdown as used by [START_REF] Schmitt | Modélisation de la dynamique des systèmes de peuplement: de SimpopLocal à SimpopNet[END_REF]; and (vi) a novel model based on thresholded deterministic potential breakdown. The network is grown until reaching a fixed number of new nodes to ensure comparability between the different heuristics.

A network is then quantified by topological indicators, namely average closeness and betweenness centralities, diameter, efficiency, and average path length. To compare the generated network to real data, we collected the European street network from the open data provided by Open-StreetMap, and extracted the value of topological indicators on 50km side moving windows. Typical configurations were then selected according to morphological classes for population density. The corresponding real network measures are compared to network generated on the same population density configurations, for 50 different configurations.

Results

The heterogeneity of models suggested the choice of an implementation in NetLogo which handles easily heterogenous multi-agent systems. Source code and results are available on the open repository of the project1 . The parameter space of the model is explored using the workflow engine Open-MOLE [START_REF] Reuillon | Openmole, a workflow engine specifically tailored for the distributed exploration of simulation models[END_REF] which allows an easy distribution of tasks on a computation grid.

We sample for each density configuration 1000 parameter points (including the generation heuristic and specific parameters for each heuristic in the random LHS sampling) and repeat the stochastic model 5 times for each point. Results are shown in Fig. 1 below. We illustrate for a synthetic population density configuration the different types of networks obtained. The feasible space covered by the sampling in the reduced indicator space (two first components of a principal component analysis) is such that heuristics are complementary to cover the full area, and overlap between point clouds quantified with a concentration index are low (Herfindhal index with a first quartile at 0.54 on a grid of size 20 separating the reduced indicator space). This complementarity is confirmed when comparing to the indicator values for real networks: some points fall relatively far from the point cloud, but each heuristic captures some real configurations. The distribution of distances of simulated points to real points confirm this low minimum for each heuristic. When conditioned by morphological class for population density, we obtain varying behavior regarding the most performant heuristic, conforming that different urban regimes are associated to specific network morphogenesis processes.

Therefore, our work demonstrates that, among the heuristic we compared, there is no best model to reproduce real configurations and that all processes included are complementary. Possible developments include a more refined calibration to find for typical configurations of real networks the heuristic and parameters producing the closest networks, or a comparison between heuristics controlling for the number of parameters to take overfitting into account (which could be here at the origin of the low performance of the biological heuristic e.g.). We also suggest that similar benchmarks should be more systematically performed for the modeling of complex territorial systems.

Figure 1 :

 1 Figure 1: (Top left) Examples of network configuration generated by the different heuristics; (Bottom Left) Feasible space of topological indicators for networks generated by the model, shown in a reduced dimension after principal component analysis, with a different color for each heuristic; (Right) Feasible space with red points corresponding to real networks, and histograms of indicator distance to data, for all densities configuration and conditioned by the type of urban form.

at https://github.com/JusteRaimbault/ CityNetwork/tree/master/Models/MesoCoevol