Qingtai Xiao

Jianxin Pan
email: jianxin.pan@manchester.ac.uk

Junwei Huang

Jianxin Xu

Hua Wang

Discrepancy measure for uniformity qualification of bubble movement in a direct contact heat exchanger

Keywords:

is presentation aims to introduce three novel methods for assessing the distribution uniformity of bubbles and determining their mixing time to uniformity. ese uniform design based methods are illustrated through image analysis in a direct-contact heat exchanger (DCHE). e novelties include the li le constraint of final evaluation value and the efficient algorithm of bubble centroid, leading to fast and accurate characterization of mixing uniformity. To determine the effects of resolution aspect (considering two types of aspect ratios: 16:9 and 4:3) and local region on the measurement, various evolutions of bubbles movement and three mixing measures in the DCHE are investigated experimentally. Real experiments and simulations are conducted, showing that the proposed methods outperform the existing methods such as Be i number. It is also shown that the location effects of bubbles can be measured successfully using these discrepancy measures, which brings a new insight into the comparison of mixing state of different systems.

INTRODUCTION

e purpose of mixing is to obtain a homogeneous mixture or a certain degree of uniformity of mixtures in non-reactive systems [START_REF] Xu | Multiphase mixing quantification by computational homology and imaging analysis[END_REF][START_REF] Dong | antifying nanoparticle mixing state to account for both location and size effects[END_REF][START_REF] Fei | A novel approach for measuring bubbles uniformity and mixing efficiency in a direct contact heat exchanger[END_REF]. Mixing uniformity has a decisive impact on the overall performance of mixing processes and can sometimes serve as a surrogate for other properties, such as the quality grade of print defects [START_REF] Lee | antification of perceived macro-uniformity[END_REF], the heat transfer performance of fluids [START_REF] Xiao | Hypothesistesting combined with image analysis to quantify evolution of bubble swarms in a direct-contact boiling heat transfer process[END_REF][START_REF] Xiao | antifying the evolution of flow boiling bubbles by statistical testing and image analysis: toward a general model[END_REF], the mechanical behavior of materials [START_REF] Berthiaux | Principal component analysis for characterising homogeneity in powder mixing using image processing techniques[END_REF], the drug content of monolithic devices [START_REF] Karami | A novel image analysis approach for evaluation of mixing uniformity in drug-filled silicone rubber matrix[END_REF], and so on.

ere is an increased desire for measuring and comparing mixing uniformity which is required from a practical point of view and for validation of theoretical models as well in various fields [START_REF] Paulo Davim | A novel approach based on digital image analysis to evaluate the delamination factor a er drilling composite laminates[END_REF][START_REF] Xiao | Measure of bubble non-uniformity within circular region in a directcontact heat exchanger[END_REF][START_REF] Kim | Measurement of mixing of two miscible liquids in a stirred vessel with electrical resistance tomography[END_REF].

In the literature, the image processing technology has been widely used for feature extraction and multiphase mixing quantification in single phase, gas-liquid, solid-liquid and gas-liquid-solid systems [START_REF] Xie | Statistical image analysis reveals features affecting fates of Myxococcus xanthus developmental aggregates[END_REF][START_REF] Hamzah | antification of moisture sensitivity of warm mix asphalt using image analysis technique[END_REF][START_REF] Xiao | Nonuniformity quantification of temperature and concentration fields by statistical measure and image analysis[END_REF]. A large number of experimental studies have been also devoted to demonstrating its efficiency in mixing time estimation for configurable optimization [START_REF] Xu | Multiphase mixing quantification by computational homology and imaging analysis[END_REF][START_REF] Dong | antifying nanoparticle mixing state to account for both location and size effects[END_REF][START_REF] Xiao | antifying the evolution of flow boiling bubbles by statistical testing and image analysis: toward a general model[END_REF][START_REF] Karami | A novel image analysis approach for evaluation of mixing uniformity in drug-filled silicone rubber matrix[END_REF].

ere is however no universally accepted image-based method for determination of mixing time mainly because each one has its own limitations, such as conductivity method, pH probe method, dual indicator system method [START_REF] Rodriguez | On the measurement and scaling of mixing time in orbitally shaken bioreactors[END_REF], tracer concentration method, coloration decoloration method [START_REF] Cabaret | Mixing time analysis using colorimetric methods and image processing[END_REF] and box-counting method [START_REF] Coent | An original image-processing technique for obtaining the mixing time: e box-counting with erosions method[END_REF], etc. Some limitations have been reported in details, e.g., in [START_REF] Xiao | Hypothesistesting combined with image analysis to quantify evolution of bubble swarms in a direct-contact boiling heat transfer process[END_REF][START_REF] Xiao | antifying the evolution of flow boiling bubbles by statistical testing and image analysis: toward a general model[END_REF][START_REF] Meng | Statistical measures of two dimensional point set uniformity[END_REF].

In earlier publications, Be i number was reported in [START_REF] Xu | Multiphase mixing quantification by computational homology and imaging analysis[END_REF] as one of the most efficient methods in determining criti-cal mixing time of mixing process [START_REF] Xu | Accurate estimation of mixing time in a direct contact boiling heat transfer process using statistical methods[END_REF] and acquiring more spatial evolutions information of flow field associated with heat transfer performance of fluids [START_REF] Huang | antifying the synergy of bubble swarm pa erns and heat transfer performance using computational homology[END_REF]. However, we found that mixing time and uniformity determined by Be i number have an issue of spatio-temporal limitation. Particularly, the mixing uniformity of bubbles in a rectangle area with the same Be i number (β 0 when white pixels are the background or β 1 in other situations) is very different. e issue of measuring the space-time uniformity of random bubbles in a recorded image was addressed by the uniformity coefficient (UC) method [START_REF] Fei | A novel approach for measuring bubbles uniformity and mixing efficiency in a direct contact heat exchanger[END_REF] and modified UC method [START_REF] Xu | A modified L 2 -star discrepancy method for measuring mixing uniformity in a direct contact heat exchanger[END_REF] in certain extents. But in the definition of reported UC methods, it is reasonable that the final evaluation value should be between 0 and 1. Otherwise, the definition may be not sensible. Fang et al. pointed out that one should choose a set of experimental points with smallest discrepancy among all possible designs of a given number of factors and experimental runs [START_REF] Fang | Uniform design: eory and application[END_REF], which satisfies the limit condition. On the other hand, the finding procedure of bubble centroid can be made efficiently using the Image Processing Toolbox (IPT) in Matlab and its in-house function, which performs be er in computational complexity than traversal operation of elements by [START_REF] Fei | A novel approach for measuring bubbles uniformity and mixing efficiency in a direct contact heat exchanger[END_REF] and [START_REF] Xu | A modified L 2 -star discrepancy method for measuring mixing uniformity in a direct contact heat exchanger[END_REF].

Inspired and motivated by the previous work, we aim to investigate the space-time feature by the definition of discrepancy directly and present an analysis of numerical simulations and experiments of the mixing process. Discrepancy concept and thresholding algorithm are applied to quantify the space-time mixing uniformity of random bubbles in a direct contact heat exchanger (DCHE) in this paper.

METHODS

Suppose there are n bubbles of interest (i.e., black or white elements) over a binary image C 2 = [1, P]×[1, Q] (i.e., a P×Q matrix mathematically). e goal here is to assess whether or not the bubbles are uniformly sca ered on the processed image. Let F u (x) = x 1 x 2 be the uniform distribution function with a given point x = (x 1 , x 2 ) where 1

x 1 P and 1 x 2 Q, and F X (x) be the empirical distribution function of bubbles set

X = x T 1 , x T 2 , • • • , x T n T , i.e., F X (x) = 1 n n i=1 I [x i ,∞] (x) = ♯(X [1, x]) n ( 1
)
where

I [x i ,∞] (x) is the indicator function and ♯(X [1, x]) denotes the number of bubbles of X falling into the local region [1, x] = [1, x 1 ] × [1, x 2 ]
. e L ∞ -star discrepancy of X on C 2 (star discrepancy for short) is defined as

D * ∞ (X) = F X (x) -F u (x) ∞ = sup x ∈C 2 F X (x) -F u (x) (2)
where x decides the area and controls the number of the given points within the area (i.e., the number of local regions), and sup denotes the superior/maximum of the local discrepancy functions (LDFs).

Statistically, the centered discrepancy (CD) and the wraparound discrepancy (WD) exhibit some advantages such as permutation invariance, rotation invariance (reflection invariance) and the ability to measure projection uniformity, but star discrepancy does not. Assume that x i1 , x k1 are the x-axis values and x i2 , x k2 are the y-axis values of ith, kth bubbles,respectively, the analytical expressions for CD and WD are given respectively as follows:

CD(X) = 13 12 2 - 2 n n i=1 2 j=1 1 + 1 2 |x i j - 1 2 | - 1 2 |x i j - 1 2 | 2 + 1 n 2 n i=1 n k=1 2 j=1 1 + 1 2 |x i j - 1 2 | + 1 2 |x k j - 1 2 | - 1 2 |x i j -x k j | 1 2 (3)
and

WD(X) = - 4 3 2 + 1 n 3 2 2 + 2 n 2 n-1 i=1 n k=i+1 2 j=1 3 2 -|x i j -x k j | + |x i j -x k j | 2 1 2 (4)
Obviously, the greater the value of the discrepancy, the more non-uniform the bubbles in a given region. In order to measure the effects of robustness and assess the performance of three measures, experiments are carried out in a DCHE. As an effective way to use energy, DCHE has been applied and researched extensively in energy recovery from industrial waste. It can transfer the heat of continuous phase (i.e., heat conduction oil) to dispersed phase (i.e., organic fluid), which can be influenced by four factors with three levels. Let E 1 -E 9 denote nine different experimental levels based on orthogonal design table L 9 (3 4 ). e images of bubbles are captured using a high-speed video camera with brand PRAKTICA of Germany and resolution 4 million pixels with no LED light. Generally speaking, the objective is bubble (i.e., gas phase) and background represents the liquid phase. For the convenience of engineering calculation and programming, we set P=1280 and Q=720 while adopting equal interval sampling of 300 which corresponds to 5 min. A processed bubbles image is selected randomly for representing the binary mixture, as shown in the inserted illustration of Fig. e novelty of the presented contribution are two-fold: using the mathematical measure of discrepancy and applying the IPT function regionprops directly, leading to low computational cost and high degree of accuracy. us, the effects of resolutions aspects and numbers of local regions on the definitions of three discrepancies in this paper are all investigated. To save space, Fig. 1 only gives the D * ∞ (X) trends over time when the pixels size reduces from 16:9 to 4:3. e reason why the results are influenced by the aspect ratio of the image is that the LDFs ratio relies on the area ratio of images with different sizes. In terms of the mixing time determination, it shows that the influence of homogenization curve by our method does not change a lot with the resolution size. Likewise, CD(X) and WD(X) are both dependent on the coordinate axes established for the bubbles image. As magnification factor (px/mm) is still similar, the aspect ratio does not have an influence on the results if adequate equipment is used. As shown in Fig. 2, there is subtle influence of number of local regions on the trend of uniformity curve.

RESULTS AND DISCUSSION

antification of mixing state

As shown in Fig. 3, quantitative comparisons of the uniformity curves determined by the proposed techniques are conducted with our previous experimental data. Specially, here there is a clear and distinctive different at the early stage of mixing process of experimental cases. Obviously, any difference in numerical performance has a deep significant role in quantifying the mixing state. e results show that good agreements of the critical point obtained by existing and proposed methods are made, see Table 1. It is observed that the effectiveness of those techniques are valid. However, the Be i number method fails to take into account the location effects of bubbles with respect to mixing transient.

Account for the location effects

It may not be correct to conclude that the best transients are that present the highest Be i number for special cases, as the purely numerical approach does not take position details into account. In practical applications, one can more accurately calculate the value of discrepancy for choosing the best experimental condition for space-time feature of bubbles, which is the main advantage of our methods.

In order to verify the feasibility of our methods, two processed images which both present 194 bubbles are taken from an integrated experiment, as shown in Fig. 4 (a) and(b). e right-hand plot in Fig. 4 (c) presents the difference of mixing uniformity in the special transient with the same Be i number β 1 =194. Comparisons show that different experimental cases with the same Be i numbers can be detected by the discrepancy measure. e combined results in Fig. 4 indicate that the proposed evaluations may be enough to obtain a good determination and quantification of mixing state of gas-liquid two-phase flow. In addition, the movement rules of bubbles in DCHE are also discussed by the method of numerical simulation. According to Fig. 5, the green regions in the three inserted images refer to the bubbles in DCHE. In particular, three sets were randomly generated and fixed with β 1 = 234 by means of the Matlab so ware. Just like usual, each image has the dimension of 1280×720. Fig. 5 also depicts that the star discrepancy would be influenced by the initial position for calculating local discrepancy function, including top-le (TL), bo om-le (BL), bo om-right (BR) and top-right (TR), but CD and WD do not. It is also interesting to point out that when one wants to assess the mixing uniformity more accurately, the CD may outperform the WD and is more sensitive to practical engineering application in some sense.

Verification of three properties

Assume centroid coordinates of bubbles are wri en as, X =

x 11 x 12 . . . . . .

x n1 x n2 (5)
where n is the number of bubbles in an image processed by classical thresholding algorithm (for example, Otsu method) mentioned earlier. As mentioned early, the first column elements x 11 , x 21 , • • • , x n1 of X correspond to the x-axis values and the second column elements x 12 , x 22 , • • • , x n2 of X correspond to the y-axis values. Let CD(X) and WD(X) denote the two different measures of uniformity individually. eoretically, CD(X) and WD(X) have their own advantages. To begin with, they are both invariant to disrupted order of the experimental points. Second, they are invariant if x i1 and x i2 are replaced individually by 1-x i1 and 1-x i2 , 1≤i≤n.

ird, the projection uniformity over all sub dimensions is taken into account, and CD(X) and WD(X) are both invariant on the subspaces. e resulting data verification are depicted in Table 2, 3, and 4, respectively. Hence, we conclude that the different experimental pa erns with the same Be i numbers can be identified by the proposed discrepancy measures. Interestingly, those provides a tool to investigate the difference of mixing state quantification between 2-dimensional area and 3-dimensional spaces.

CONCLUSIONS

In this paper, the mixing state of gas-liquid two-phase flow in a DCHE for waste-heat utilization is quantified through three novel statistical measures with reference to discrepancy. e proposed methods do not rely on the perceived uniformity (resulting from relaxing the range constraint of final evaluation values, and the finding algorithm of bubble centroid based on [START_REF] Fei | A novel approach for measuring bubbles uniformity and mixing efficiency in a direct contact heat exchanger[END_REF] and [START_REF] Xu | A modified L 2 -star discrepancy method for measuring mixing uniformity in a direct contact heat exchanger[END_REF]. Summarily, three assertions for this narrating work are list as follows.

(1) With respect to the star discrepancy, the local discrepancy function of a set of bubbles seems to be a useful concept to measure the irregularity of the distribution of bubbles within a rectangle region. e scheme is based on partitioning of image data obtained from a high-speed video camera, thus leading to a new consideration in the interpretation of flow visualization.

(2) Using the three image analysis techniques (star discrepancy, CD and WD) processed in the ImageJ or Matlab so ware, the evolution of bubbles movement is tracked experimentally. e influences of expressions, iteration steps and pixels on the mixing efficiency are also discussed. Compared to [START_REF] Fei | A novel approach for measuring bubbles uniformity and mixing efficiency in a direct contact heat exchanger[END_REF] and [START_REF] Xu | A modified L 2 -star discrepancy method for measuring mixing uniformity in a direct contact heat exchanger[END_REF], the proposed approaches are easy to implement and are computationally low cost.

(3) e CD and WD are more reasonable than the star discrepancy in terms of mixing state quantification in this current study. at is to say, the former two could satisfy the three properties including invariance to permutation, invariance under reflection, and projection uniformity, but the later one does not.

Based on the aforementioned observations, it is believed that the proposed approaches can be applied to measure mixing uniformity with a good accuracy in pharmaceutical, chemical, metallurgical, printing and medical industries, etc.

Figure 1 .

 1 Figure 1. Effect of different resolution on star discrepancy

1

 1

4 Figure 2 .

 42 Figure 2. Effect of different local regions on D * ∞ (X)

Figure 3 .

 3 Figure 3. Comparison of three discrepancies curve for E 2 .

Figure 4 .

 4 Figure 4. Mean absolute discrepancy evolutions of different experimental images with the same β 1 =194.

Figure 5 .

 5 Figure 5. Different discrepancy measures for quantifying the particle mixing state in three different synthetic images with the same β 1 =234.

Table 1 .

 1 Mixing time obtained by different methods.

	methods β 1	D * ∞ (X) CD(X) WD(X)
	E 1	156	151	151	151
	E 2	93	97	95	96
	E 3	168	172	170	170
	E 4	225	224	225	224
	E 5	122	120	120	120
	E 6	84	84	82	83
	E 7	262	259	195	198
	E 8	117	118	118	118
	E 9	128	127	127	127

Table 2 .

 2 Verification of invariance to permutation

	Discrepancy	original order	disordered bubble coordinates	both
	CD(X)	0.0249	0.0249	0.0249	0.0249
	WD(X)	0.0277	0.0277	0.0277	0.0277

Table 3 .

 3 Verification of invariance under reflection

	Discrepancy	no reflected	x = 1 2	reflected y = 1 2	both
	CD(X)	0.0249	0.0249 0.0249 0.0249
	WD(X)	0.0277	0.0277 0.0277 0.0277

Table 4 .

 4 Verification of projection uniformity improved the presentation. is work was financially supported by National Natural Science Foundation of China (Nos. 51666006, 51406071 and 51706195), Joint Funds of the National Natural Science Foundation of China (No. U1602272), Scientific and Technological Leading Talent Projects in Yunnan Province (No. 2015HA019) and Academician Workstation of ZHANG Wenhai in Yunnan Province (No. 2015IC005).

	Discrepancy	no projected	y=0	x=0	projected to origin
	CD(X)	0.0249	0.6025 0.6009	0.9446
	WD(X)	0.0277	0.4720 0.4717	0.6812

ACKNOWLEDGMENTS

e authors wish to extend special thanks to anonymous reviewers for numerous detailed questions and constructive


Alik Palatnik
email: umapalatnik@gmail.com

Tommy Dreyfus
email: tommyd@tauex.tau.ac.il

Students' reasons for introducing auxiliary lines in proving situations

Keywords: High school geometry, proof, auxiliary lines, reasoning, task-based interviews

We study proving situations in geometry where students introduce auxiliary lines. We examine students' reasons for introducing auxiliary lines. The overarching theme of the tasks proposed to the students is the comparison of areas of triangles and/or parallelograms. The students gave a wide spectrum of reasons when introducing auxiliary lines. Two main groups of reasons have been discerned for modifying a given situation by introducing auxiliary lines: recalling a learned procedure using some known results or definitions, and anticipating more information from the modified situation. We present a case study of a pair of students who combine reasons of recalling and anticipating nature when introducing auxiliary lines.

Introduction

The objective of this study is to contribute to research on proving in geometry. Specifically, we investigate the introduction of auxiliary elements by students while proving and characterize reasons students give for their actions. Our interest in this process stems from the following facts: First, proving in geometry is difficult for students in general. Second, the introduction of auxiliary lines contributes substantially to students' difficulty with proofs in geometry (e.g. [START_REF] Hsu | Cognitive complexity of mathematics instructional tasks in a Taiwanese classroom: An examination of task sources[END_REF][START_REF] Senk | How well do students write geometry proofs[END_REF]. Our focus on students' reasons for introducing auxiliary elements led us to build on the seminal work of [START_REF] Pólya | How to solve it[END_REF], who highlighted the role of auxiliary lines in problem solving.

Theoretical background

It is widely recognized that in secondary school, proof and proving are still mostly taught in the area of geometry [START_REF] Herbst | Interactions with diagrams and the making of reasoned conjectures in geometry[END_REF][START_REF] Fujita | Reasoning-and-proving in geometry in school mathematics textbooks in Japan[END_REF]. According to [START_REF] Battista | The development of geometric and spatial thinking[END_REF] school geometry presents students with a complex network of concepts, representations and ways of reasoning "…that is used to conceptualize and analyze physical and imagined spatial environments" (p. 843). Combining deductive and spatial reasoning is pertinent to proving activities in geometry. Indeed, in his study of American geometry classes, [START_REF] Herbst | Interactions with diagrams and the making of reasoned conjectures in geometry[END_REF] demonstrated, among others, a generative type of students' interaction with diagrams in proving activities. Such interactions are characterized by deductive reasoning interrelated with modification of a given diagram. Such interactions are virtually absent from geometry classes, as studies of Japanese secondary schools confirm [START_REF] Miyazaki | Flow-chart proofs with open problems as scaffolds for learning about geometrical proofs[END_REF] but are central for proofs that require the introduction of an auxiliary element.

According to [START_REF] Pólya | How to solve it[END_REF], "an element that we introduce in the hope that it will further the solution is called an auxiliary element" (p. 46). In our study, a slightly expanded version of this definition is used. Namely, we consider also the modification of an already existing element as an auxiliary element (see Task 1 in the Methodology section). Moreover, an auxiliary element may be only imagined rather than physically drawn.

According to [START_REF] Senk | How well do students write geometry proofs[END_REF] the introduction of auxiliary lines is a critical part in the solution of geometric proof problems. She concluded that the difficulty with auxiliary lines "exemplifies the need to teach students how, why, and when they can transform a diagram in a proof" (p. 455). [START_REF] Ding | Teaching geometry in lower secondary school in Shanghai, China[END_REF] also identified the introduction of an auxiliary line as a source of great difficulty. [START_REF] Hsu | Geometric calculations are more than calculations[END_REF] claims that this difficulty is caused by students' need to perceive the diagrams dynamically and apply transformational observation to visualize a solution achievable with auxiliary lines.

Students' difficulty with introducing auxiliary lines when proving can be partially explained by a form of instruction, which does not encourage students to add new lines to the given diagram. For instance, students in the study of [START_REF] Herbst | Proving and doing proofs in high school geometry classes: What is it that is going on for students[END_REF] self-reported that they expect to play a passive role in interaction with diagrams when proving: They believe that they are usually not expected to introduce auxiliary lines, and if so, they expect to receive a hint to do so from the teacher. This difficulty may be overcome by the incorporation of various types of exploration in proving which makes the introduction of auxiliary lines more transparent (e.g. [START_REF] Fan | Does a transformation approach improve students' ability in constructing auxiliary lines for solving geometric problems? An interventionbased study with two Chinese classrooms[END_REF].

Considering these studies and Pólya's fundamental work, it makes sense to consider the introduction of auxiliary lines through the prism of reasons that led to their introduction. [START_REF] Pólya | How to solve it[END_REF] lists three main reasons for introducing auxiliary lines (p. 46): 1) Trying to use known results;

2) Going back to definitions;

3) Expecting to make the original problem "fuller, more suggestive, more familiar".

Pólya makes a distinction between the first two reasons and the third one. He states that in the third case "we scarcely know yet explicitly how we shall be able to use the elements added. We may just feel that it is a 'bright idea' to conceive the problem that way with such and such elements added" (ibid.). In our study we sought both, to obtain descriptive reports about students introducing auxiliary elements in proving situations, and to identify the types of reasons they provide for introducing them. We adapt Pólya's classification to our aims by introducing the terms recalling (Pólya 1, 2) and anticipating (Pólya 3) reasons.

Given this background, our research aims to answer the following questions:

A. What reasons do the students give when introducing an auxiliary line in a proving situation? B. To which type do these reasons belong, recalling or anticipating?

Methodology

Students and Tasks

The participants in our study were high school students studying geometry in Grade 10-the first grade of Israeli high school. Eighteen students were nominated by their mathematics teachers; ten of them volunteered for the study. The students were grouped into five pairs, following [START_REF] Schoenfeld | Purposes and methods of research in mathematics education[END_REF] suggestion that students working in pairs will show reliable and natural behavior.

We sought to create tasks where students would be inclined to draw auxiliary lines in order to make progress. We anticipated that tasks with multiple solutions, each based on different auxiliary lines, will elicit a variety of reasons for introducing auxiliary lines. The topic of areas was suitable as it was already known to the students, but proofs on areas had not been taught in class at the time of the intervention. We used two tasks with different levels of complexity.

Task 1: Quadrilateral ABCD is a parallelogram. Where on side AB do you have to place point E, so that the area of triangle DEC is a half of the area of the given parallelogram ABCD? Explain your reasoning.

There are multiple ways to show that E can be anywhere on side AB. Figure 1 hints at some of them, highlighting the role of appropriate auxiliary lines in the proofs. Task 2: Prove that the quadrilateral formed by the midpoints of the sides of a given quadrilateral is a parallelogram, and that its area is half of the area of the given quadrilateral (Varignon's theorem, Theorem 3.11 in [START_REF] Coxeter | Geometry revisited[END_REF]. Task 2 contains two parts, the second of which relates to area. This two-part structure and the large variety of possible auxiliary lines allowing to make use of properties of a triangle's midsegment and of a parallelogram to prove the claim make it especially suitable for our study (see students' sketches in Figure 2 in the Findings section).

Procedure

We used two interview techniques to obtain information on the introduction of auxiliary lines: stimulated recall (SR) and partial revealing (PR). SR is a reflective method used in qualitative research (e.g., [START_REF] Bikner-Ahsbahs | Towards the emergence of constructing mathematical meanings[END_REF]) in which students carry out a task undisturbed and are later prompted by the use of some visual or auditory stimulus to recall their thoughts while carrying out the task. A crucial advantage of SR is that it does not interfere with carrying out the task. Thus, we decided to use SR in the relatively simple first task. Also, the SR sessions which followed the first task familiarized students with the researchers' request to reflect upon their actions.

In the more complex Task 2, we wanted to make sure that students do not get stuck. For this purpose, we developed and used the idea of partial revealing (PR), a variation on task-based interviews [START_REF] Goldin | A scientific perspective on structured, task-based interviews in mathematics education research[END_REF]: We prepared a progressively more revealing sequence of 12 visual or verbal hints. The hints were descriptions/ sketches of an auxiliary line or some key idea/ definition, supporting different ways to prove Task 2 (cf. hints to the same problem in de Villiers, 2014). The hints helped students recalling relevant knowledge or focusing on specific information given in the problem. After obtaining a hint, the students were encouraged to reflect upon the new situation created by the hint. This approach helps the students (and the researcher) to understand and describe the possible reasons for introducing an auxiliary line, whether they tried to use a definition, a previously solved problem or just feel that the hint will advance the proof.

In a first meeting (about 20 minutes), each pair of students worked on Task 1 (Figure 1) without researcher intervention. The solution, explanation, or proof, agreed upon by the two students was collected by the researcher without comment. In the second meeting (about 45 min), two days after the first one, the student pairs returned to Task 1, using SR. In the third meeting (about 45 min), the students worked together on Task 2. The PR method was applied.

The qualitative analysis of the data was based on the methodology used by [START_REF] Jones | Deductive and intuitive approaches to solving geometrical problems[END_REF] for the exploration of students' solutions of geometric construction problems on the subject of interplay between students' intuitive and deductive reasoning. First, the video-recordings of all interviews and the students' written productions were searched for auxiliary lines drawn by the students. Then, the videos were transcribed; we read the transcripts of all the meetings line by line, and marked phrases that seemed to relate to reasons for adding auxiliary elements. Next, students' actions and utterances while introducing the auxiliary lines were coded according to our adaptation of [START_REF] Pólya | How to solve it[END_REF] list of reasons for introducing of auxiliary lines. For instance, the utterance "The subject of the problem is an area, we need to draw the height" which accompanied introduction of the line segment EF (Figure 1c) was coded as "recalling" (going back to definition/ concretizing a formula).

The utterance "It just seems right to draw this line" about the same auxiliary line was coded as "anticipating" (expecting to make the original problem more suggestive). These broad categories of reasons were further divided according to the type of auxiliary line and the type of reason.

When analyzing Task 1, we juxtaposed students' explanations during the SR sessions with their explanations to their partners during the problem-solving meetings and with written productions and sketches. We paid special attention to the general context of the proof and in particular to the change in students' actions and utterances after introducing the auxiliary line; in this, we followed [START_REF] Goldin | A scientific perspective on structured, task-based interviews in mathematics education research[END_REF] way of analyzing interactions between the subjects and mathematical structures of the tasks. This procedure helped us to confirm and elaborate our initial coding.

When analyzing Task 2 we paid special attention to the moments where specific PR hints (definitions, formulas or visual images) led to the introduction of an auxiliary line; again we followed [START_REF] Goldin | A scientific perspective on structured, task-based interviews in mathematics education research[END_REF], this time with respect to the role of hints in the analysis of students' problem-solving behavior.

Findings

We analyzed five pairs of students. Here1 , we present the case of Bob and Boris whose reasons were of both types: "recalling" and "anticipating".

Task 1

Boris soon drew the auxiliary line GF (Figure 1d) stating: "I want to calculate the area of CED and then the area of the two smaller triangles, to compare." In the SR session, he described GF as a height through the midpoint of AB, and gave the following reasons for the introducing GF: "Because they ask for area, and the height is a part of the area formula" (Bob&Boris SR10). When asked about the reason for drawing the height exactly at the midpoint of AB, Boris can't give a specific reason: "It was the first point to draw the height"; "it was an initial thought" (Bob&Boris SR 28, 36). In the SR session Bob recalled that immediately after Boris drew the auxiliary line, he thought that there probably more places where the area is the same.

The students' solution of Task 1 culminated in the following exchange:

Boris:

The height will always be the same! Bob: Let's assume that I place it, point E, together with point A [Points at the sketch, draws imaginary AC].

Boris:

[loud] Look, it is at every point! Bob:

[continues with confidence] So, it happens at every point.

Boris:

The base is the same base, and the height is the same height, so at every point it will be equal.

Bob: U-hu, if I place it there, at point A (draws imaginary line AC=EC), it is diagonal, it halves the parallelogram. If I place it at the point B, diagonal BD also halves the parallelogram.

Bob: So, it looks to me, that at every point. (Bob&Boris M1 35-41)

Analysis of Task 1. The auxiliary elements introduced by the students in the course of solving the task were: (1) GF-the height through the midpoint of AB;

(2) and (3) diagonals AC and BD (imagined). As follows from the students' explanation given above, the reason to use height was that it forms part of the formula for the area and thus of recalling nature. A possible explanation for the choosing the height through the midpoint of the side is the intuitive rule same A-same B [START_REF] Stavy | How students (mis)understand science and mathematics: Intuitive rules[END_REF], namely, half of the area (A) -half way on the side (B). As for the diagonals AC and BD, Bob imagined them and gestured them in the air. In the SR session, he explained that if E is placed at the vertices A and B it's easy to see that the diagonals AC and BD halve the parallelogram (anticipating).

Task 2

The students started by drawing their own sketch (Figure 2a). Boris soon focused on the second part of the task: the area of the parallelogram is half that of the quadrilateral. Bob drew EG, a height in ABCD, as first auxiliary line, but the students realized that they do not know a formula for the area of a general quadrilateral. They received as hints (a) an accurate sketch of the task; (b) a drawing of a triangle with a midsegment. The midsegment immediately triggered progress. Boris drew the diagonals AC and BD, and Bob formulated a full formal proof to the first part of the task, using the property that a midsegment of a triangle is parallel to the third side of the triangle.

To prove the second part of the task, the students wanted to calculate areas and looked for a place to construct a height. A switch from calculating to comparing areas was initiated when they attended to pairs of congruent triangles and in particular to the triangles AHN and HDI (Figure 2b,2c). Bob suggested comparing their areas to that of quadrilateral HIMN, and Boris introduced the auxiliary line NI (Figure 2c). Later the students developed this idea into a proof (Figure 2d). We may assume that the auxiliary line EG, perceived by the students as a height, was introduced as a part of a recalled procedure which connects between height as an auxiliary line and height as a part of the area formula as previously mentioned by Boris.

The diagonals AC and BD were introduced immediately after the visual hint (a drawing of a triangle with a midsegment). We assume that the hint led the students to recall the definition and properties of a midsegment, which became a reason for introducing AC and BD which in turn helped them elaborate the sketch into a formal proof.

We interpret that the reason for introducing the crucial auxiliary line NI was an anticipating one. Despite the fact that Boris was not sure how this auxiliary line will advance his proof, he introduced it wishing to show that the areas of quadrilateral HIMN and the combined areas of HDI and AHN are equal. It is remarkable that in the collaborative problem-solving efforts of Bob and Boris in Tasks 1 and 2 the reasons to introduce auxiliary lines were both of recalling and anticipating nature.

Discussion

Our aim was to reveal and characterize students' reasons for introducing auxiliary lines in geometric proof situations. The five pairs of students who participated in the study introduced a total of 26 auxiliary elements of various nature: heights of triangles and parallelograms, diagonals of parallelograms, triangle midsegments, and medians. Some of these elements were actually drawn, others imagined. As we expected, Task 2 was not easy for students to prove (cf. [START_REF] Ding | Teaching geometry in lower secondary school in Shanghai, China[END_REF][START_REF] Hsu | Cognitive complexity of mathematics instructional tasks in a Taiwanese classroom: An examination of task sources[END_REF] about the complexity of tasks with auxiliary lines).

The reasons students gave when introducing auxiliary elements are of two types (1) recalling some known formulas, definitions and procedures, and (2) anticipating to receive more information from a modified situation. There were 14 auxiliary elements of the first type and 12 of the second.

In the tasks we used, most of the reasons of a recalling nature were related to the introduction of heights. The students drew heights as part of a learned algorithm for finding areas: height is explicitly present in the area formula; hence students recalled the formula, and almost immediately drew heights as auxiliary lines; this seems to be a pertinent feature of students' proving process in the context of areas. It was very salient in the interview with Bob and Boris.

Our data allow us to elaborate on characteristics of Pólya's third type of reason for introducing auxiliary lines: anticipating to make problem fuller, more suggestive, without knowing exactly how the auxiliary line will be used. Most of the examples of such reasons we found can be organized into two groups. In the first group (4 occurrences), we put reasons related to intuitive rules [START_REF] Stavy | How students (mis)understand science and mathematics: Intuitive rules[END_REF]. Evidence for this is provided by the intuitive "half side-half area" reasoning when introducing the auxiliary lines through the midpoint of the side (e.g., Bob and Boris, Task 1, line GF). In the second group of anticipating reasons for introducing auxiliary lines (6 occurrences), we put reasons related to the students' feeling that they exhausted the available resources and they needed to produce some new source of information (e.g., Bob and Boris, Task 2, NI).

All students used area formulas when proving Task 1; a majority also introduced heights in Task 2, and when this failed, they introduced other auxiliary lines. One possible explanation for this is, that when the students felt that they had exhausted the resources suggested by recall, they tried auxiliary lines based on some anticipation, while keeping in mind that they have to substantiate it (as happened to Bob and Boris when introducing NI in the second part of Task 2).

Based on these findings, we suggest that the introduction of auxiliary lines in proving situations can be expected to provide a rich context for the development of geometrical thinking. Teachers should encourage proactive and reflective approaches to the introduction of auxiliary elements. It is impossible to develop mastery related to the introduction of auxiliary lines without creating experiences of drawing, imagining and dynamically modifying auxiliary lines. This study of students introducing auxiliary elements when proving included only a small number of participants, and only two tasks from a single geometry topic, and it was conducted with paper and pencil only. Further research is required to better understand the complex dynamics of introducing auxiliary lines.

  Figure 1. (a-c) Task 1 with auxiliary lines, (d) students' sketch

  a diagonal [of HIMN]. And [points to the sketch] this angle equals this, 'cause it is parallel. Bob: No, you still have to prove this. Boris: Wait a minute [rotates the sheet with the sketch]. We have to prove that these areas are equal [points to the triangles AHN, HDI, NMI and IHN]. It will make it a midsegment (of triangle AMD) (Bob&Boris M3 134-140).

Figure 2 .

 2 Figure 2. The students' sketches for Task 2 (2c is our copy of a relevant part of students' sketch 2b)

Analyses of two more cases are presented elsewhere(Palatnik & Dreyfus,

2018)


Mohamed Benkirane

François Clautiaux

Jean Damay

Boris Detienne

A Hypergraph Model for the Rolling Stock Rotation Planning and Train Selection

Keywords: rolling stock rotation planning, train selection, mixed integer programming, hypergraph model

This paper is about an integrated optimization approach for timetabling and rolling stock rotation planning in the context of passenger railway traffic. Given a set of possible passenger trips, service requirement constraints, and a fleet of multiple heterogeneous self-powered railcars, our method aims at producing a timetable and solving the rolling stock problem in such a way that the use of railcars and the operational costs are minimized. To solve this hard optimization problem, we design a mixedinteger linear programming model based on network-flow in an hypergraph. We use this models to handle effectively constraints related to coupling and decoupling railcars. To reduce the size of the model, we use an aggregation and disaggregation technique combined with reduced-cost filtering. We present computational experiments based on several French regional railway traffic case studies to show that our method scales successfully to real-life problems.

Introduction

Railway production is based on the use of several resources such as infrastructure, staff, and rolling stock. A service planning assembles all these resources while respecting technical, organizational and legal constraints. Many working rules also have to be followed. Interactions between these different resources and constraints complicate the process of planning. This process is generally decomposed into sequential subproblems in the following order: line planning, timetable planning, rolling stock rotation, and crew scheduling. The process of line planning corresponds to the definition of lines based on the rail network. The line planning includes the number and the routes of lines. It also includes the choice of the line frequencies (see [START_REF] Schöbel | Line planning in public transportation: models and methods[END_REF]). In the timetable planning, the decision maker chooses when trips will be operated, using the lines designed in the previous step. In the rolling stock rotation problem, one defines the train composition for each trip defined in the timetable, i.e., one associates to each train a set of railcars and specifies when coupling and decoupling occurs. Finally, crew scheduling corresponds to the process of assigning crews of working agents to trains. Rolling stock optimization is key in our problem. This problem was largely discussed in the past recent years. Several different versions of the problem have been studied, depending on whether maintenance constraints, order and orientation of the railcars are considered or not. There are two main types of formulations in the literature: pathbased and flow-based formulations. Path based models integrate easily unit specific constraints (e.g maintenance constraints), but produce exponentially large formulations that have to be generated dynamically. Flow-based models are generally not able to follow individual railcars, but can be solved directly by a general-purpose MIP In this example we can operate a train from Bordeaux to Agen each hour, the defined cycle time is equal to one hour. For instance 00h02, 01h02, ... , 22h02 with an arrival at 01h01, 02h01, ..., 23h01 respectively solver. For a recent survey of the literature, see [START_REF] Jørgen T Haahr | A comparison of two exact methods for passenger railway rolling stock (re) scheduling[END_REF], where the authors compared the two types of models using testbeds from different countries. Among the most recent approaches, the path based model presented by [START_REF] Jørgen T Haahr | A comparison of two exact methods for passenger railway rolling stock (re) scheduling[END_REF] considers the order of train units in a composition. It extends the model presented in [START_REF] Thorlund Haahr | A branch-and-price framework for railway rolling stock rescheduling during disruptions[END_REF]. A branch-and-price algorithm generates the path variables dynamically by solving a sub-problem that corresponds to a shortest path problem with resource constraints. Many network-flow models are extensions of the seminal work of [START_REF] Fioole | A rolling stock circulation model for combining and splitting of passenger trains[END_REF]: [START_REF] Kjaer Nielsen | A rolling horizon approach for disruption management of railway rolling stock[END_REF], [START_REF] Cacchiani | Railway rolling stock planning: Robustness against large disruptions[END_REF] and [START_REF] Kroon | Rescheduling of railway rolling stock with dynamic passenger flows[END_REF] for example. One of the most recent flow-based model was presented by [START_REF] Kjaer Nielsen | A rolling horizon approach for disruption management of railway rolling stock[END_REF], based on the work described in [START_REF] Fioole | A rolling stock circulation model for combining and splitting of passenger trains[END_REF]. The authors present a multi-commodity flow model where each commodity is a train unit type, which takes into account railcar inventory constraints at the end of the day as well as changes in the composition of trains. Network-flow models in hypergraphs natively embed train composition constraints within the hypergraph, without adding additional linking constraints. In addition one can associate to multiple train units a cost that is non proportional to the cost of a train unit. [START_REF] Borndörfer | Integrated optimization of rolling stock rotations for intercity railways[END_REF] solve the rolling sock problem based on such an hypergraph model. They integrate in their models regularity constraints, maintenance constraints and train composition. Each node of the hypergraph represents a departure (or arrival) train station, a departure (or arrival) time, a composition and a train unit appearing in this composition and its orientation. The hyperarcs represent either trips or coupling/decoupling operations. Hyperarcs representing trains link a set of departure nodes and a set of arrival nodes, those nodes represent the same trip and the same composition. Hyperarcs representing train re-compositions link a set of arrival nodes of a trip and a set of departure nodes of another trip such that the number of train units for each type is conserved. The latter hyperarcs integrate also the regularity of the rolling stock. The mixed integer program integrates maintenance constraints through another flow commodity which is absorbed by trip hyperarcs and generated by service nodes. At the national French railway company (SNCF) the timetable planning is decomposed in two phases: netgraph design and train selection. The concept of netgraph is used to describe all possible trips that can be produced. The netgraph is made available by the infrastructure manager for all train operators, which justifies how the train selection problem arises at SNCF Mobility, which represents the train operator, and at all its potential competitors. The netgraph defines a set of train patterns (see Figure 1), i.e., an origin and a destination train stations, and a set of stops. For each train pattern, a set of possible departure and arrival minutes are given. Those minutes are translated to possible starting times through a given cycle time of the pattern (typically every one or two hours).

To operate a train, an operator needs to purchase a train path on a specific departure time from the infrastructure manager. The netgraph design focuses on the "clock-face" scheduling, which is known in the literature as the cyclic train timetabling problem (see e.g. [START_REF] Serafini | A mathematical model for periodic scheduling problems[END_REF][START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF][START_REF] Cacchiani | Nominal and robust train timetabling problems[END_REF]). The train selection selects in the netgraph a set of trains that satisfies the passenger demand and public authorities requests. In this work, the term train stands for a trip defined by specific stopping points between origin and destination at precise times; the term railcar stands for a self-propelled passenger railway vehicle, also called train unit. Given a netgraph and train frequency requirements, our approach aims at solving the train selection and rolling stock problems in an integrated manner. This problem occurs at the strategic level, when the company decides which trains will be operated. The objective is to minimize the cost of the railcars used, and some penalties for not satisfying demand constraints. An originality of our approach is to integrate the choice of the selected trains within the process of optimizing the rolling stock. Since we deal with a strategic problem, precise operational constraints and maintenance are not considered, although an estimation of production costs is accounted for. Trip demand is also defined in an original way, which arises from the opening up to competition and from the netgraph formalism. We use a network-flow model on a hypergraph to model and solve our problem. This model allows to effectively deal with coupling/decoupling constraints and costs. Our hypergraph formulation differs from [START_REF] Borndörfer | Integrated optimization of rolling stock rotations for intercity railways[END_REF] for several reasons: one flow unit can be related to a single train unit, but also to several combined units. Likewise, nodes are associated with combined train units. We do not integrate maintenance constraints.

In terms of solving approach, we use a variant of the aggregation technique presented in [START_REF] Borndörfer | Integrated optimization of rolling stock rotations for intercity railways[END_REF]. In our method, we first solve with an MIP solver a simplified problem where all railcars are considered equivalent. We use this first phase to filter out many variables from the initial problem and produce a first feasible solution. Both techniques are used to speed up the solution of the initial problem. We describe our problem in section 2, and propose a hypergraph-based mixed integer model in section 3. We present in section 4 our solving approach. Finally, we present our computational experiments in section 5.

Problem description

In the present work, we focus on combining the train selection and the rolling stock problems, i.e. we take the line planning and the netgraph as inputs, and do not consider constraints related to crew scheduling. The goal of our optimization problem is to evaluate for a specific set of lines and in a rapid manner the cost of the involved rolling stock. Note that the cyclic timetabling problem and the netgraph design are not addressed by our approach.

Physical network and rolling stock

Let S be the set of train stations. The infrastructure resources enable to operate train units from one station of S to another one at given times of the horizon.

Definition 1 (train pattern). A train pattern p is a tuple (s - p , s + p , S p ), where s - p , s + p ∈ S respectively represent departure and arrival stations, and S p the set of stations where stops occur.

The set of train patterns is denoted as P . Train patterns can only be operated at certain times, which are allocated by the infrastructure manager. A train pattern assigned to a certain time is called train in the remainder of this paper. The set of possible trains is denoted as V.

Definition 2 (train). A train i is a tuple (p i , t - i , t + i ), where p i ∈ P is a train pattern, t - i and t + i the departure and arrival times for p i .

For the rolling stock rotation problem, we consider that we have at our disposal a set of heterogeneous railcars. Each railcar has a railcar type r, which is characterized by a length, a cost ψ(r), and a payload µ(r). We denote by R the set of railcar types. Railcars are distributed over a set of warehouses, which are located at some given stations. For s ∈ S, r ∈ R, we denote respectively by n min r,s and n max r,s the minimum and maximum number of railcars of type r that begin the service from s in a feasible solution. The number of railcars beginning the service from a train station must be equal to the number of railcars ending the service at the same train station. This is required to allow a cyclic planning. A railcar composition c is an ordered list of railcar types. We denote by C the set of feasible railcar compositions. The number of elements in the list is denoted len(c). When two compositions c 1 and c 2 are concatenated, we use the classical notation c 1 + c 2 . With a slight abuse of the notation, we will denote µ(c) = r∈c µ(r) the payload of composition c. Individual railcars are considered as simple compositions. Compositions are created from other compositions that are coupled or decoupled. A composition can be used at a given station only if it is already present at this station, or the stock of compositions at this station is sufficient to create the composition. Operational constraints limit the creation of a new configuration to only one coupling or decoupling operation at train station platforms. Stationdependent constraints may also forbid such operations. However, any composition can be created from single units at the depots, and vice-versa. Each coupling (resp. decoupling) operation incurs a cost denoted by δ + and δ -. A composition is not allowed to stay idle at a station platform more than a given duration τ s , s ∈ S: after that time period, it is moved to the depot, as a set of individual railcars.

Constraints on train selection

Deciding whether each train is operated by a vehicle or not is the main decision in our problem. To be operated, a train has to be assigned a composition. A railcar composition can only be used to operate a train if feasibility conditions are met, i.e the length of the composition does not exceed the length of the train station platforms; tracks and powering types must be compatible to railcar types composing the train composition. As part of the problem definition, some pairs of trains are mutually exclusive. For each i ∈ V, we denote by E i the set of trains that cannot be operated if train i is (independently of the composition used). If two trains share a same portion of the trip on departure or on arrival, they can be operated together, which reduces the operational costs. In this case, they are coupled (or decoupled) at an intermediate station. When a coupling operation is performed, it means that two trains are first operated each with a given composition. After the intermediate station, the concatenation of the two compositions is used to operate both trains over their common portion. Pairs of trains that can be coupled or decoupled are specified in the data. We denote by V + the set of ordered pairs of trains sharing the arrival portion; V -represents the set of ordered pairs of trains sharing the departure portion. A railcar, appearing in two different railcar compositions, can only be used for two successive trains if the trains and their compositions are compatible, i.e. it is practically feasible. Some pairs of trains are not compatible in any composition. Some are only compatible in some compositions. Compatibility depends on the station where the successive trains are operated. Several practical constraints are taken into account (arrival times, time of coupling or decoupling railcars, stub-end stations etc.). Moreover, operating consecutive (set of) trains with some given compositions has a cost. The cost for operating one or two trains together with one or two compositions (in the case of coupling/decoupling within the trip) depends on several factors. To simplify the exposition, we introduce β((I k , c k ) k ) as the cost of operating trains in sets I k with their respective compositions c k . If this association is not possible, then β((I k , c k ) k ) = +∞. This notation is later used to express the cost of hyperarcs in our model.

Mission demands

The selection of trains is subject to passenger demand constraints coming from issues and public authorities requests. The frequency distribution of trains and the planned capacity have to meet demand requirements, which are expressed as missions.

Definition 3 (mission).

A mission m is a tuple (p m , ν min m , θ m , [t min m , t max m ]) where p ∈ P is a train pattern, [t min m , t max m ] a time interval, ν min m ∈ N the minimum payload capacity for m and θ m is a function determining the cost of underor overcovering the mission demand compared to the desired frequency.

We denote by M the set of missions. Since the possible trains are known a priori, for each m ∈ M, the set V m of trains which can cover the demand of mission m is given as an input. Train i covers mission m if p i = p m and

t + i ∈ [t min m , t max m ]
. Mission demands are expressed as follows. During interval [t min m , t max m ], the number of chosen trains that cover m must be in [f min m , f max m ]. Additionaly, the sum of the payload capacities of the train compositions used to operate them must be at least ν min m . For each mission m ∈ M, the following piecewise linear cost function θ m : N → R + expresses the penalty associated to not covering m with the desired frequency f c m .

θ m (f ) =        π +∞ m (f -f max m ) + π + m (f max m -f c m ) if f ≥ f max m π + m (f -f c m ) if f c m ≤ f ≤ f max m π - m (f c m -f ) if f min m ≤ f < f c m π -∞ m (f min m -f ) + π - m (f c m -f min m ) if f ≤ f min m where f min m ∈ N, f c m ∈ N, f max m
∈ N are respectively the minimum, ideal, and maximum frequency for the mission, π + m ∈ R + and π - m ∈ R + respectively the unit costs of over-and undercovering the demand, and π -∞ m ∈ R + , π +∞ m ∈ R + respectively the unit costs of critical over-and undercovering the demand. The first and last cases respectively happen when the frequency is higher than the maximum or lower than the minimum allowed. In both In some cases and for some political issues it is mandatory to operate some specific trains. The set of mandatory trains is denoted by V f .

A network flow formulation 3.1 Hypergraph description

In [START_REF] Borndörfer | Integrated optimization of rolling stock rotations for intercity railways[END_REF], the authors show that using hyperarcs allows to represent train compositions, coupling and decoupling operations in an efficient manner. We also use hypergraphs in our model: vertices correspond with logical/spacetime states where decisions can be made, while hyperarcs correspond with railcar operations (technical moves or operating trains). One unit of flow corresponds with one or several train units from one place to another, possibly after rearrangement of their configuration. Our directed hypergraph is defined by a couple of sets (V, A) where V is the set of nodes and A the set of hyperarcs. Each hyperarc a ∈ A is defined by two multisets: T (a) is the tail of a and H(a) its head. In this paper, we denote by card(u, N ) the multiplicity of element u in multiset N , and by support(N ) the set of unique elements in N . The sum over a set Ξ of multisets is the multiset ξ = N ∈Ξ N such that support(ξ) = ∪ N ∈Ξ support(N ) and card(u, ξ) = N ∈Ξ card(u, N ). Each hyperarc a has a cost q(a). The reader is referred to [START_REF] Cambini | Flows on hypergraphs[END_REF] for more definitions on hypergraphs. A node is characterized by a tuple (s, t, c, I), where s is a station, t a time, c a railcar composition, and I a set of trains. Basically, such a node corresponds with being in station s at time t with composition c, with the additional information that it corresponds with operating trains in I. For each train i ∈ V , two train nodes are created for each possible composition c that can operate i:

(s - i , t - i , c, {i}) and (s + i , t + i , c , {i})
, where c is the resulting composition if the train was operated with composition c. For each pair of trains (i, j) in V -, and each compatible composition c, a node (s + i , t + i , c, {i, j}) is created. Similarly, for each pair of trains (i, j) in V + , and each compatible composition c, a node (s - i , t - i , c, {i, j}) is created. There are also depot nodes, which are related to technical stays. For these special nodes, c = (r) is composed of a single railcar type r, and I = ∅. There is one depot node for each element of S × T × R. In addition, two artificial vertices α and Ω are added to V to represent respectively the source and the sink of the network. Hyperarcs are used to represent the decisions that can be made. We distinguish three types of hyperarcs: train hyperarcs, which represent the possibility of selecting trains with a specific composition, connection hyperarcs, to connect consecutive trains and/or rolling stock in the technical stations, and technical hyperarcs, to ensure flow conservation for each railcar type in the technical stations. Depending on the number of railcars used in the corresponding composition, each hyperarc can have one, two, or three elements in its tail/head. All hyperarcs have unit capacity except technical arcs, which are bounded by the quantity of available rolling stock. Train hyperarcs are related to the main decision variables (i.e. train selection). A hyperarc can represent operating one train, or two trains when a coupling/decoupling operation is performed during the travel. When one

v 1 v 2 A single train from V v 3 v 4 v 5
A pair of trains from

V + v 6 v 7 v 8
A pair of trains from V - Connecting an arrival with a departure

v 3 v 4 v 5
Connecting two arrivals with a unique departure

v 6 v 7 v 8
Connecting a unique arrival with two departures

v 1 v 2 v 3
Connecting one arrival with two identical railcars from the depot to ensure a departure

v 4 v 5 v 6
Connecting one arrival to a departure and sending one railcar to the depot

v 7 v 8
Combining three identical railcars from the depot to ensure a departure

v 9 v 10 v 11
Sending two railcars of different types, comming from an arrival, to the depot Figure 4: Square nodes are depots nodes, double circled nodes are arrivals nodes and unique circled nodes are departure ones train is involved, train hyperarcs are created as follows. For each i ∈ V, and each compatible composition c, there is a hyperarc a such that T (a) = { {v i } } and H(a) = { {w i } }, where v i = (s - pi , t - i , c, {i}), w i = (s + pi , t + i , c , {i}), and c is the arrival composition obtained if c is used to cover i. Note that c can be either c or its symmetric in case of stub-end intermediate stops. The cost q(a) is β(({i}, c)). Train hyperarcs also represent the selection of two trains i and j that are either coupled or decoupled. For all (i, j) ∈ V + and each possible combination c of compositions c i and c j , a hyperarc a is created, such that

T (a) = { {v i , v j } }, H(a) = { {w ij } }, where v i = (s - pi , t - i , c i , {i}), v j = (s - pj , t - j , c j , {j}), and
w ij = (s + pi , t + i , c , {i, j}). Similarly, for each pair of trains (i, j) ∈ V -, a hyperarc a is created, with T (a) = { {v ij } }, H(a) = { {w i , w j } }, w i = (s + pi , t + i , c i , {i}), w j = (s + pj , t + j , c j , {j}), v ij = (s - pi , t - i , c , {i, j}).
In the two latter cases, q(a) = β(({i}, c i ), ({j}, c j ), ({i, j}, c )). Figure 3 illustrates the different types of train hyperarcs. Connection hyperarcs are used to connect consecutive trains and/or rolling stock. We aggregate all technical requirements for connection between operations in a single function, which is defined as follows. Possible connections are pairs of multisets (e 1 , e 2 ) where e 1 and e 2 are multisets of V . We define a function γ, which takes a pair (e 1 , e 2 ) in parameter and produces a positive real value if the connection is possible, +∞ otherwise. Various technical constraints are embedded in the definition of γ, such as shunting (move to depot), maximum stay time at platforms τ s , coupling/decoupling and minimum required time for connection. For each pair of multisets (e 1 , e 2 ) (each of cardinality one, two, or three), there exists a ∈ A with H(a) = e 1 and T (a) = e 2 if and only if γ(e 1 , e 2 ) < +∞. Figure 4 illustrates different types of connection hyperarcs. Finally, technical hyperarcs are added to take into account the rolling stock constraints. For each depot node v, there is a hyperarc connecting α to v. This arc is used to position vehicles at each technical station s at the beginning of the time horizon. Its cost is equal to ψ(r). Additionally, for each t = 0, . . . , T -1, there is a hyperarc of cost zero connecting (s, t, (r), ∅) to (s, t + 1, (r), ∅) This arc represents the fact that railcars of type r stay at technical station s between time t and t + 1. Finally, there is a hyperarc of cost zero connecting each depot node v to Ω. This arc is used to measure the quantity of vehicles of type r at each station at the end of the time horizon. To reduce the number of nodes in the model, we apply the following preprocessing method. Two nodes related to trains arriving at the same time, in the same station entrance, with the same composition, are merged. Moreover, if a vertex has only one predecessor and one successor, then it is removed, and the two corresponding hyperarcs are merged.

Mixed integer programming formulation

Based on hypergraph (V, A), we define the following MILP model of the problem, which will be referred to as TSRP in the sequel. For i ∈ V, we denote by A(i) ⊂ A the subset of hyperarcs that are related to operating train i: a hyperarc a ∈ A(i) represents either i or a pair of trains belonging to V -∪ V + and including i. For i ∈ V, a ∈ A(i), we denote by N (a, i) = µ(c), where c is the composition related to train i in hyperarc a, the number of seats represented by the train composition modeled by hyperarc a to cover train i. Also, we denote by a α r,s the hyperarc linking the source node with the first 1 depot node representing the train station s ∈ S and the railcar type r ∈ R; likewise a Ω r,s is the hyperarc linking the last 1 depot node with the sink node. We associate to each hyperarc a ∈ A an integer variable denoted by z a representing the flow handled by a. For all v ∈ V , we respectively denote by Γ + (v) and Γ -(v) multisets of outgoing and ingoing hyperarcs. For each mission m ∈ M, we define four continuous variables y - m , y + m , y -∞ m and y +∞ m where (y - m + y -∞ m ) represents under frequency coverage and (y + m + y +∞ m ) represents over frequency coverage. We report below our mixed-integer programming formulation for the TSRP problem. ( 5)

a∈A(i) z a = 1 ∀i ∈ V f (6) i∈Vm a∈A(i) N (a, i)z a ≥ ν min m ∀m ∈ M (7)
a∈A(i)

z a ≤ 1 ∀i ∈ V (8)
a∈A(i)∪A(j)

z a ≤ 1 ∀i ∈ V, ∀j ∈ E i (9) z a α r,s = z a Ω r,s ∀r ∈ R, ∀s ∈ S (10)
n min r,s ≤ z a α r,s ≤ n max r,s ∀r ∈ R, ∀s ∈ S (11)
y + m , y - m , y +∞ m , y -∞ m ∈ R + ∀m ∈ M (12)
z a ∈ {0, 1} ∀a ∈ ∪ i∈V A(i) (13)
z a ∈ N ∀a ∈ A \ ∪ i∈V A(i) (14)
Objective function (1) minimizes the production cost and the penalty related to deviating from the desired demand covering (see figure 2). Constraints (2) represent the flow conservation. Constraints (3) express the demand covering, (y - m + y -∞ m ) and (y + m + y +∞ m ) are respectively equal to the amount of frequency under-and over-covering. Constraints [START_REF] Caprara | Modeling and solving the train timetabling problem[END_REF] impose that one hyperarc must be activated from a set of hyperarcs referring to the same train within V f . Constraints [START_REF] Fioole | A rolling stock circulation model for combining and splitting of passenger trains[END_REF] impose a minimum number of seats that must be available for each mission. Constraints (4) and ( 5) bound the value of (y - m ) and (y + m ). Constraints [START_REF] Jørgen T Haahr | A comparison of two exact methods for passenger railway rolling stock (re) scheduling[END_REF] impose that at most one hyperarc can be activated from a set of hyperarcs referring to the same train. Constraints (9) avoid activating two hyperarcs representing two conflicting trains. Constraints [START_REF] Kroon | Rescheduling of railway rolling stock with dynamic passenger flows[END_REF] impose for each train station s and for each railcar type r the number of railcars starting the service at s is equal to the number of railcars ending the service at s. Constraints [START_REF] Kjaer Nielsen | A rolling horizon approach for disruption management of railway rolling stock[END_REF] impose that the number of railcars, for each type, beginning the service in a specific station s must not exceed n max r,s and must not be lesser than n min r,s .

Solving approach

Our preliminary computational experiments have shown that model ( 1)-( 14) can be solved efficiently for real-life instances when only one type of railcars is considered. However real instances include several types (typically three), and including them in the model increases by a wide range the computing time needed to solve the model. Our method relies on an initial relaxation of the problem, which is refined in a second step to include the totality of the problem. In the relaxation we use, we consider a smaller number of railcar types. The goal of this first relaxation is twofold: it allows us to eliminate a large number of hyperarcs using reduced-cost fixing techniques, and it is used in a heuristic method to compute a first feasible solution.

Our relaxation is defined by an aggregation of the network: all nodes related with the same station, time, number and class of railcars are aggregated.

F ∈P A max a∈F card T (a) (W ) X F = F ∈P A max a∈F card H(a) (W ) X F ∀W ∈ P V (16) -y + m -y +∞ m + y - m + y -∞ m + i∈Vm F ∈P A :A(i)∩F =∅ X F = f cm m ∀m ∈ M (17) y + m ≤ f max m -f c m ∀m ∈ M (18) y - m ≤ f c m -f min m ∀m ∈ M (19) F ∈P A :A(i)∩F =∅ X F = 1 i ∈ V f (20) F ∈P A :A(i)∩F =∅ max b∈F N (b, i) X F ≥ ν min m ∀m ∈ M (21) F ∈P A :A(i)∩F =∅ X F ≤ 1 i ∈ V (22) F ∈P A :(A(i)∪A(j))∩F =∅ X F ≤ 1 ∀i ∈ V, ∀j ∈ E i (23) X F α R,s = X F Ω R,s ∀R ∈ P, s ∈ S ( 24
) r∈R n min r,s ≤ X F α R,s ≤ r∈R n max r,s ∀R ∈ P (25) y + m , y - m , y +∞ m , y -∞ m ∈ R + ∀m ∈ M (26) X F ∈ {0, 1} ∀F ∈ P A : F ∩ A V = ∅, X F ∈ N ∀F ∈ P A : F ∩ A V = ∅ (27) Proposition 1. Model (ATSRP) is a relaxation of Model (TSRP).
Proof. We first need an intermediate result to perform a variable substitution in constraints ( 3), ( 6), ( 8) and ( 9). Since P A is a partition of A, a∈A(i) z a = F ∈P A a∈A(i)∩F z a . Besides, for all F ∈ P A and a, b ∈ F , we have by definition of P A that agg V P (H(a)) = agg V P (H(b)) and agg V P (T (a)) = agg V P (T (b)). Hence, a ∈ A(i) if and only if b ∈ A(i), so that either A(i) ∩ F = F or A(i) ∩ F = ∅ and (28) follows.

For all i ∈ V,

a∈A(i) z a = F ∈P A :A(i)∩F =∅ X F ( 28
)
We also need to express the flow on hyperarcs going out of or into a vertex aggregate in terms of those new variables.

v∈W a∈Γ + (v) z a = v∈W a∈A card(v, T (a))z a (29) = F ∈P A a∈F v∈W card(v, T (a))z a (30) = F ∈P A max a∈F card T (a) (W ) X F (31)
Equality (29) follows from the definition of Γ + . The second line (30) is deduced from the fact that P A is a partition of A, while (31) comes from the definitions of X F and card T (a) . Remark that for all a ∈ F , the values of card T (a) (W ) are the same and are thus equal to the maximum of them. Note that the choice of the maximum over other possibilities such as the minimum is arbitrary. Similarly, we can prove the following relation: (2)

For all W ∈ P V , v∈W a∈Γ -(v) z a = F ∈P A card H(a),a∈F (W )X F (
⇒ v∈W a∈Γ + (v) z a = v∈W a∈Γ -(v) z a ∀W ∈ P V \ {{α}, {Ω}} ⇔ F ∈P A card T (a),a∈F (W )X F = F ∈P A card H(a),a∈F (W )X F ∀W ∈ P V \ {{α}, {Ω}} (from (31)-(32))
Based on the definition of P A , we can also make the following observation.

For all R ∈ P, for all s ∈ S, {a α r,s , r ∈ R } ∈ P A and {a

Ω r,s , r ∈ R } ∈ P A . ( 33
)
Following the same idea, constraints [START_REF] Kroon | Rescheduling of railway rolling stock with dynamic passenger flows[END_REF] and ( 11) are summed up over each railcar type aggregate:

(10) ⇒ r∈R z a α r,s = r∈R z a Ω r,s ∀R ∈ P, s ∈ S ⇔X F α R,s = X F Ω R,s ∀R ∈ P, s ∈ S (11) ⇒ r∈R n min r,s ≤ r∈R z a α r,s ≤ r∈R n max r,s ∀R ∈ P ⇔ r∈R n min r,s ≤ X F α R,s ≤ r∈R n max r,s

∀R ∈ P

with F α R,s = agg a P (a α r,s ) and F Ω R,s = agg a P (a Ω r,s ) for some arbitrary r ∈ R, for all R ∈ P (from (33))

Constraints [START_REF] Fioole | A rolling stock circulation model for combining and splitting of passenger trains[END_REF] are relaxed in the following way:

(7) ⇒ i∈Vm e∈A(i) max N (b, i) : b ∈ agg a P (e) z e ≥ ν min m ∀m ∈ M ⇔ F ∈P A :A(i)∩F =∅ max N (b, i) : b ∈ F X F ≥ ν min m ∀m ∈ M (from (28))
Note that disjunctive constraints (22) induce binary bounds on variables X F with F ∈ P A such that F ∩ A V = ∅ in constraints (27). We now consider the objective function.

a∈A q(a)z a ≥

F ∈P A a∈F min{q(b) : b ∈ F }z a = F ∈P A min{q(b) : b ∈ F }X F
We have shown that any solution for the initial model is a solution of the new model, and that for any solution, the new objective function is less than the initial one. Therefore we obtain an relaxation.

The ATSRP and TSRP models have the same types of constraints. Therefore, the ATSRP problem is seen as a TSRP problem based on different data inputs. We show that it is the case by defining from the ATSRP the new artificial railcar types (denoted by R a ) and train configurations (denoted by C a ). Each R ∈ P R defines a new railcar type r ∈ R a equipped with the most permissive characteristics among the types of R ; the set C a is defined likewise based on P C . Solving the ATSRP is then equivalent to solving a TSRP defined by R a and C a .

Our iterative algorithm

solving approach, sketched in Algorithm 1 has three main steps: solving a relaxation, computing a feasible solution, and solving exactly the problem using information obtained in the two first phases. Lines 1, 2 and 3 correspond with defining and solving the linear relaxation of the TSRP model with inputs (R a , C a , V, V + , V -, M); λ * is the value of the linear relaxation. We then compute in line 4 and 5 the value σ of a feasible solution by first solving the relaxed model to optimality. We extract from this solution the subset of trains used (J , J + , J -) (such that J ⊆ V, J + ⊆ V + and J -⊆ V -). Second, we solve in line 5 the TSRP restricted to that subset of trains, i.e. with inputs (R, C, J , J + , J -, M). We then remove all hyperarcs of the relaxed model by applying the reduced cost-based variable fixing technique described in Proposition 2. We finally solve the initial TSRP model restricted to the remaining hyperarcs.

Algorithm 1: The iterative algorithm based on the reduced cost variable fixing procedure Data: (R, C, V, V + , V -, M) 1 P ← Determine an initial aggregation of the railcar types ; 2 (R a , C a ) ← Compute the artificial railcar types and train configurations based on P; 3 λ * ← Compute the value of the linear relaxation of the TSRP(R a , C a , V, V + , V -, M); 4 (J , J + , J -) ← Solve the TSRP(R a , C a , V, V + , V -, M) and compute the selected trains; 5 σ ← Compute an upper bound by solving the TSRP(R, C, J , J + , J -, M); 6 Remove all hyperarcs of the TSRP(R a , C a , V, V + , V -, M) whose reduced cost is larger than (σ -λ * ); 7 Generate the filtered original TSRP(R, C, V, V + , V -, M); 8 Solve the filtered original TSRP(R, C, V, V + , V -, M);

We now show that hyperarcs of the original model whose variable has been removed by the reduced-cost fixing (line 6) in the aggregated model can be removed in the original model. By construction, model (ATSRP) is a relaxation of model (TSRP) and the model variables are linked by the following property, which is ensured by the constraint-wise process used in the previous section to build (ATSRP). Proposition 2. Let (RATSRP) be the linear relaxation of (ATSRP), and let σ be an upper bound on the optimum of (TSRP), (X * , y * ) be an optimal solution of (RATSRP), λ * its cost, and cF be the reduced cost of a variable X F such that X * F = 0. If λ * + cF > σ then for all solutions of (TSRP) with a cost smaller than or equal σ, we have z a = 0 for all a ∈ F . Proof. Let F ∈ P A such that λ * + cF > σ and consider a feasible solution (z, y) of TSRP such that z a ≥ 1 for some a ∈ F . From the proof of Proposition 1, (X, y) with X G = a∈G z a for all G ∈ P A is a feasible solution to ATSRP and γ(z, y) ≥ γ(X, y). Moreover, X F = a∈F z a ≥ 1, so that γ(X, y) > σ from standard reduced cost-based variable fixing result. It follows that γ(z, y) > σ. Finally, ẑa = 0 for any solution (ẑ, ŷ) of TSRP of cost less than σ.

Computational experiments

We present in this section our experimental results carried out on real life instances. The objective of our experiments is twofold: (i) to evaluate the quality of our formulation and (ii) to compare the filtering based algorithm and the initial one. All experiments were run using a 2.5 Ghz Haswell Intel Xeon E5-2680 with 128Go of RAM. The MIPs are solved using the commercial MILP solver IBM Ilog Cplex 12.6.

Instances are distinguished by the size of the set of trains, the set of mission demands, the set of railcar types and the set of compositions. We also study parameter τ s , s ∈ S, which reflects the maximum stay duration alongside a train station platform. This last one impacts the size of the hypergraph. We recall that flow units are merged in depot nodes depicting a move of railcars to depot areas. τ s is set to around an hour for instances with major train stations having several train departures and arrivals. It is set to around two hours for small train stations with restrained trips. The time period of the study is equal to one day for all instances and the size of R is bounded by three, in accordance with current practice of the French operator over a specific part of the network. Some trains have a huge demand of payload capacity, we therefore allow train compositions with up to three combined railcars. These compositions are not desired as they are practically difficult to manage, so we penalize them in the objective function. We measure the percentage of these compositions in C and denote it by %MU3. Table 1 summarizes the instances in our testbed. For each instance, we generated some variants, named in column Data, by changing parameters R, C and τ . Column TP represents the number of train patterns and τ = max s∈S (τ s ) is expressed in minutes.

In the following, the gap between a lower bound l and an upper bound p is denoted by Gap(l, p) and is equal to 100 * |l-p| |p| . The speedup between an initial CPU time i and a new one n is equal to 100 * i-n i and is denoted by S CPU , the speedup for the hypergraph generation phase is denoted by S build . To evaluate the impact of operating train compositions with three railcars, we vary the fraction ξ of the cost of those train compositions over the penalty π - m related to the demand covering.

Solving TSRP directly

Table 2 shows the results of model TSRP (subsection 3.2) using using a general purpose commercial MILP solver. In columns #Ctrs and #Vars, we respectively report the number of constraints and variables in the model. Column CPU1 represents the overall amount of time in seconds needed to solve the TSRP and %CPU Build represents the percentage amount of time needed to generate the hypergraph model. Columns Opt and LP respectively contain the optimum and the linear relaxation value for the problem. Here, the value of ξ is set to 0.7. Although the models have up to 2 000 000 variables, the quality of the linear relaxation to allow a fast resolution. All instances are solved to optimality (within the solver's default tolerance equal to 10 -4 ). The hardest instance is solved in 32 hours, while 15 out of 20 of them are solved within one hour. Compared to the overall computing time, the percentage of hypergraph generation time is small except for small instance I7 and instance I5. 3 represents the results where we significantly penalize the use of train compositions with three railcars. In both tables, CPU2 represents the overall computing time in seconds needed to solve TSRP using algorithm 1. %Gen represents the percentage of the generated hyperarcs compared with #Vars. LPA is the value of the linear relaxation of the aggregated problem, Primal is the value obtained using the primal heuristic presented in the algorithm 1, S CPU and S build are the speedup of the overall computing time and the build time. d build is the time difference between the initial build and the build of the filtered hypergraph.

Overall, we managed to speed up consequently the solution time, with quite better results when compositions with three railcars are more penalized (Table 3). This is explained by the large part of hyperarcs that are filtered in this case. For instance I3 (D5, D6 and D7) we speed up the solution time even if we generated the whole initial hypergraph. This is explained by the tight value of the primal solution. However, almost no variable is filtered because of the poor linear relaxation obtained for in this case. Also, data sets D16 and D20 draw no benefit from our algorithm. For D16, the quality of both the primal solution and the linear relaxation are poor. For D20, the number of possible configurations is very small, leading to a small hypergraph, which is not reduced further by the filtering algorithm. Note that the amount of time needed to solve them is already small, 8 seconds for D16 and 406 seconds for D20. Considering the primal heuristic, good solutions are obtained for most of the instances except for D16 where the fixed set of trains obtained in the aggregated model leads to a surplus of two railcars.

Conclusion

In this work, we have developed an optimization approach to deal with a new problem integrating the train selection problem and the rolling stock rotation problem (TSRP). The proposed algorithm is based on a network flow in a hypergraph model which made it possible to easily apprehend coupling and decoupling multiple units. An originality of our model is that flow units represent both railcar units in the technical hyperarcs and train configurations in the other hyperarcs. Also, merging flow units in depots nodes reduces the size of the generated hyperarc set. To reduce the size of the model, we implemented an algorithm based on an aggregation and disaggregation technique combined with reduced-cost filtering. The computational experiments show that this algorithm is numerically efficient for real size instances. Future works may include a broader integration in the railway production problem, including netgraph design decisions, or maintenance.

Figure 1 :

 1 Figure 1: Train pattern in the netgraph.In this example we can operate a train from Bordeaux to Agen each hour, the defined cycle time is equal to one hour. For instance 00h02, 01h02, ... , 22h02 with an arrival at 01h01, 02h01, ..., 23h01 respectively

Figure 2 :

 2 Figure 2: Representation of function θ m , m ∈ M

Figure 3 :

 3 Figure 3: Double circled nodes are arrivals nodes and unique circled nodes are departure ones

v 1 v 2 Figure 5 :

 25 Figure 5: A technical hyperarc

Figure 6 :

 6 Figure 6: Steps leading to the hypergraph aggregation

Table 1 :

 1 Characteristics of instances

	Instance Data	TP	|V|	|V + | + |V -|	τ	|M|	|R|	|C|	%MU3
	I1	D1	30	426	0	50	63	2	14	57
		D2	30	426	0	50	63	3	47	70
	I2	D3	31	620	0	50	175	2	14	57
		D4	31	620	0	50

Table 2 :

 2 Solving directly the initial TSRP

	Instance	Data	#Ctrs	#Vars	CPU1	%CPU Build Gap(LP, Opt)
	I1	D1	13211	139255	924	0,54	4,57
		D2	38814	530663	12161	0,62	4,38
	I2	D3	19175	463959	29361	0,17	5,31
		D4	44655	1190217	79316	0,41	5,32
	I3	D5	17557	196889	1024	0,98	42,58
		D6	17557	310065	3787	0,61	42,71
		D7	52088	733240	15939	0,80	42,58
	I4	D8	19759	114598	821	0,49	1,51
		D9	19759	387908	3968	0,91	1,49
		D10	58795	1518731	112967	0,47	1,50
	I5	D11	20904	846961	425	39,06	2,03
		D12	43204	2102872	2984	34,01	6,73
	I6	D13	10178	145421	80	7,50	13,23
		D14	10178	247017	200	8,00	13,06
		D15	22214	690238	2284	4,90	9,64
	I7	D16	4648	80661	9	22,22	4,63
		D17	19888	729407	621	24,32	0,98
	I8	D18	70309	2245103	3474	36,87	0,45
		D19	22501	426573	1165	4,12	0,65
		D20	16633	168505	540	1,48	0,65

Table 3 :

 3 Comparing the direct solving and the filtering procedure with ξ from {6, 8, 10, 15} We now compare the reduced cost based filtering procedure with the direct solution. This time we tested instances of table 2 with different values of ξ between 0.7 and 15. Columns of table 3 and table 4 represent average values obtained by fixing ξ from the set {6, 8, 10, 15} for table 3 and fixing it from the set {0.7, 1, 2, 4} for table 4; the over-lines stands for averages. Basically, table

	Data Gap(Opt, Primal) Gap(LPA, Primal) Gap(LPA, Opt) %Gen	CPU2	S build (d build )	S CPU
	D1	0,19	5,22	5,04	49,35	111,25	61,16 (3,22)	84,35
	D2	0,16	5,00	4,85	41,60	1024,00	71,17 (49,09)	90,07
	D3	0,26	6,40	6,16	58,46	5185,25	51,07 (26,15)	72,18
	D4	0,25	6,40	6,16	55,75	15152,75 52,02 (171,85)	83,21
	D5	0,01	43,01	43,00	99,90	241,75	-5,25 (-0,51)	77,05
	D6	0,02	43,15	43,14	100,00	525,75	-1,95 (-0,45)	76,31
	D7	0,01	43,01	43,00	99,90	1660,75	-1,13 (-1,48)	82,69
	D8	0,02	1,82	1,80	27,79	123,50	79,78 (2,93)	79,47
	D9	0,02	1,79	1,77	24,81	129,25	92,01 (32,71)	97,39
	D10	0,02	1,79	1,76	17,59	858,50	96,37 (530,70)	99,24
	D11	16,95	54,98	45,79	18,26	154,00	96,17 (160,10)	71,57
	D12	24,44	45,21	27,48	14,67	1133,50	97,68 (990,99)	82,25
	D13	13,65	44,35	35,55	37,72	56,00	75,78 (4,49)	48,88
	D14	8,48	41,35	35,91	20,77	58,75	94,17 (15,02)	86,43
	D15	7,57	25,97	19,90	16,49	317,50	96,84 (126,89)	87,37
	D16	33,57	57,77	36,42	26,88	8,50	77,04 (1,41)	-3,47
	D17	0,86	11,99	11,22	15,19	136,50	97,35 (124,40)	81,86
	D18	2,12	17,33	15,54	18,93	2692,75 96,25 (1244,08) 24,66
	D19	1,48	5,93	4,52	39,33	703,00	82,67 (41,00)	27,76
	D20	1,45	5,90	4,52	99,64	974,75	-3,84 (-0,31)	-2,53
	5.2 Reduced cost filtering procedure

Table 4 :

 4 Comparing the direct solving and filtering procedure with ξ from {0.7, 1, 2, 4}

	Data Gap(Opt, Primal) Gap(LPA, Primal) Gap(LPA, Opt)	%Gen	CPU2	S build (d build )	S CPU
	D1	0,19	5,21	5,04	98,39	221,75	-3,82 (-0,20)	74,89
	D2	0,14	4,98	4,85	98,41	980,50	0,05 (0,04)	90,33
	D3	0,25	6,65	6,41	97,33	2812,50	4,01 (2,06)	85,79
	D4	0,25	6,65	6,41	98,24	11157,00	3,06 (10,07)	76,56
	D5	0,01	43,01	43,00	99,90	248,00	-4,01 (-0,39)	76,41
	D6	0,02	43,15	43,14	100,00	601,50	-2,06 (-0,47)	76,30
	D7	0,01	43,01	43,00	99,90	1764,25	-1,18 (-1,51)	87,11
	D8	0,02	1,82	1,80	98,56	143,50	-17,59 (-0,65)	78,28
	D9	0,02	1,79	1,77	98,70	294,00	1,01 (0,36)	92,41
	D10	0,02	1,79	1,76	99,48	4057,00	1,88 (10,30)	95,82
	D11	16,92	54,97	45,79	78,59	455,00	25,76 (42,90)	-1,20
	D12	23,69	44,66	27,48	60,56	2594,75 47,50 (482,46) 16,82
	D13	15,03	45,24	35,55	100,00	114,75	-3,64 (-0,22)	-15,29
	D14	8,48	41,35	35,91	99,79	164,50	-0,08 (-0,01)	7,37
	D15	7,68	26,05	19,90	78,85	603,00	24,08 (31,15)	70,94
	D16	33,61	57,79	36,42	100,00	12,75	-8,41 (-0,15)	-45,49
	D17	0,86	11,99	11,22	19,26	146,75	96,69 (136,00) 78,71
	D18	2,06	17,28	15,54	78,82	3654,50 25,87 (332,17) 22,48
	D19	1,39	5,85	4,52	54,23	800,00	63,28 (31,97)	50,34
	D20	1,48	5,93	4,52	99,41	1032,25	-4,72 (-0,39)	-54,63

The first and the last depot nodes are seen with respect to an increasing hourly order of depot nodes for each train station and for each railcar type.

Acknowledgments

This work was partly funded by an ANRT (http://www.anrt.asso.fr/) CIFRE scholarship between University of Bordeaux and SNCF. Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from Bordeaux INP, LABRI and IMB and other entities: Conseil Régional d'Aquitaine, Université de Bordeaux, CNRS and ANR in accordance to the programme d'investissements d'Avenir (see https://www.plafrim.fr/).

This type of aggregation is tighly related to the method used in [START_REF] Borndörfer | A Coarse-To-Fine Approach to the Railway Rolling Stock Rotation Problem[END_REF]. The authors use simultaneously three models with different levels of precision, called layers. In the most precise model, which solves the complete problem, arcs and hyperarcs are considered, and the orientation and position of each railcar in the configuration is recorded. The intermediate model disregard the position and the orientation of the railcars. The less precise model does not consider hyperarcs anymore, which leads to a minimum-cost flow problem. The three layers are used in a column-generation method, where constraints of the fine level are kept in the master problem, and coarse levels are used to price out new columns. This decomposition into layers is well adapted to the version of the problem of [START_REF] Borndörfer | Integrated optimization of rolling stock rotations for intercity railways[END_REF]: the coarse problem is easy to solve and produces a good lower bound on the cost of operating railcars. The intermediate layer allows to produce feasible parts of solutions to be added to the master problem. The same idea is used differently in our work. We also identified the fact that solving the precise model directly is computationally challenging, and removing fine grained constraints does not hurt too much in terms of objective function. However, two main differences between our problem and the one solved by [START_REF] Borndörfer | Integrated optimization of rolling stock rotations for intercity railways[END_REF] explain why we used another variant of aggregation. First, since we have to solve the train selection and the rolling stock problems simultaneously, the analog of the coarser relaxation of [START_REF] Borndörfer | A Coarse-To-Fine Approach to the Railway Rolling Stock Rotation Problem[END_REF] would remain an NP-hard minimum cost flow problem with additional demand covering constraints. Second, we do not consider maintenance constraints and thus our flow problem is easier to solve even when hyperarcs are considered.

The relaxed problem

Our relaxed problem is defined by an aggregation of the set of railcar types, which cascades down to an aggregation of the vertices and the hyperarcs. In the model, each set of aggregated vertices/hyperarcs is represented by a representative, single element. Figure 6 illustrates the successive aggregations. Let P be a partition of set R, with P = {R 1 , . . . , R |P| }. For all railcar types r 1 ∈ R, we denote its part -or aggregate -by agg r P (r 1 ) = R i such that r 1 ∈ R i . • Train configurations. The aggregate of a train configuration is the set of configurations with the same aggregated railcar types. The aggregate of c = (

The aggregate of a vertex u = (s, t, c, I) ∈ V is defined as the set of vertices sharing the same characteristics with u except for their configurations, which share the same aggregate: agg v P (u) = {(s, t, c , I) ∈ V : c ∈ agg c P (c)} ⊆ V . The aggregate of α and Ω are defined respectively by agg v P (α) = {α} and agg v P (Ω) = {Ω}. • Hyperarcs. In order to define the aggregation of the hyperarcs, we need to introduce the minimal set of vertex aggregates that covers a set U of vertices: agg V P (U ) = {agg v P (u) : u ∈ U }. We have to define an aggregate cardinality function as well: given a multiset N and a set of vertices U , card N (U ) = u∈U card(u, N ) is the cumulated cardinality over elements of U in multiset N . The aggregate of a multiset is then:

We produce a relaxation of our problem by replacing the vertices of our hypergraph model with vertex aggregates, and the hyperarcs with hyperarc aggregates. Let P V = {agg v P (u) : u ∈ V } and P A = {agg a P (a) : a ∈ A} be the aggregates of V and A. This results in a smaller size hypergraph model, but requires to relax some constraints and the objective function. Formally, let us introduce new variables X F = a∈F z a , F ∈ P A , which correspond with the total flow on hyperarcs in each aggregate F . For each railcar aggregate R ∈ P and each train station s ∈ S, we denote by F α R,s = agg a P (a α r,s ) and F Ω R,s = agg a P (a Ω r,s ) with r an arbitrary element of R. A mixed integer linear programming model for the relaxed problem is obtained from (1)-( 14) by replacing z-variables with X-variables, partially relaxing some constraints in the process. This yields the following MILP model of our relaxed problem, which we call (ATSRP).

(ATSRP) : min γ(X, y) =