François Cheviré

Franck Tessier

Roger Marchand

Optical properties of rare-earth oxides and oxynitrides, …towards "second generation" UV absorbers

Keywords: UV-absorber, absorption edge, optical properties, oxynitride

Second-generation" UV-blocking materials are needed. Using the fluorite-type rare-earth tungstates R 6 WO 12 as starting reference compounds, anionic or cationic substitutions were achieved in (R,R')-(W,Mo)-(O,N) systems, thus delimiting large solid solution domains. The continuous shift of the absorption edge in diffuse reflectance spectra towards lower energy values makes it possible to tune it at the required value of 3.1 eV (400 nm), with in several cases a steep enough slope. Attempts were also carried out in perovskite (A,R)-Ti-(O,N) systems (A = alkaline-earth, R = rare-earth) to isolate suitable oxynitride compositions.

Introduction

It is well known that UV radiations can cause irreversible damages to organic materials, human tissues included. The so-called "UV B" (wavelengths within the range 280-320 nm) are the most energetic, that often results in short term effects (sunburns). However, the "UV A" (320-400 nm), about which one has not cared for a long time, are not so harmless because of cumulative effects. So, an effective protection in the whole UV spectrum is more than necessary for outer applications. The filtering system can be incorporated into the material to be protected (plastics) or can cover it (varnish, skin protection).

To date, inorganic anti-UV are not totally satisfactory. Mainly TiO 2 -based, rutile or anatase phase, they are often not effective at high wavelengths and engender whitening phenomena due to the high value of their refractive index (rutile: n = 2.72, anatase: n = 2.5) [START_REF] Yabe | Cerium oxide for sunscreen cosmetics[END_REF]. Potential challengers, with lower indices, are ZnO (n = 2.2) and CeO 2 (n = 2.05) whose absorption edge is also close to 400 nm (or 3.1 eV), the frontier between the UV part and the visible part of the electromagnetic spectrum.

The powder morphology is a crucial parameter which has to be adapted to each case. In particular, to avoid whitish aspect, nanosized particles are needed for TiO 2 , while sizes of about 100 nm are sufficient for ZnO or CeO 2 .

Furthermore, a major drawback of these "first-generation" UV-absorbers is their photocatalytic activity: under irradiation they generate oxidizing species which induce secondary reactions at the surface of the host matrix, with noticeable degradation.

The purpose of this study is to finalize the preparation of "second-generation" UV-blocking materials which are effective for the whole UV spectrum. A good candidate for such an application should satisfy the following specifications:

1/ absorption of all radiations with < 400 nm 2/ transmission of all the lower energy radiations 3/ as a consequence of 1) and 2), a steep absorption edge, as much as possible, defining the spectral selectivity 4/ a low refractive index in the Visible (however, a high refractive index in the UV gives a maximum opacity) 5/ no photocatalytic activity. Modification in the chemical composition of a solid is often accompanied by changes in its physical or chemical properties. In case of a progressive modification inside a solid solution domain, it is possible to conceive new materials with properties particularly well fitted for given applications by adjustment of the composition. Considering the optical absorption properties of semi-conducting compounds, a slight modification in composition can affect the width of the bandgap (and the position of the absorption edge) through a modification of the electronic structure.

In this paper, two types of solid solutions have been studied. Large composition domains have been evidenced in fluorite-type systems R-(W,Mo)-(O,N) (R = rare earth) [START_REF] Larcher | Recherche de nouveaux matériaux anti-UV[END_REF] and in perovskite-type systems (A,R)-Ti-(O,N) (A = alkaline-earth) by involving cationic and/or anionic substitutions. In the case of fluorites, starting from a rare-earth tungstate absorbing in the UV, the width of the bandgap can be decreased to 3.1 eV (400 nm) by modifying the position either of the valence band (N/O substitution) or of the conduction band (Mo/W substitution). In the case of perovskites, N 3-substituting progressively for O 2-in a titanate ATiO 3 induces an increase in the anionic formal charge which can be compensated, according to a cross-substitution principle, by a trivalent rare-earth cation replacing the divalent alkaline-earth. Such a substitution results in a shift of the absorption edge towards higher wavelengths and in the possibility to adjust it to 400 nm.

Experimental

Preparation of rare-earth-based oxide precursors by the ceramic method has shown limits regarding homogeneity and reactivity of the products. Single phases are often difficult to prepare in a pure state and the difficulty increases with multinary oxides, as numerous other compositions are potentially accessible.

Several chimie douce-type processes have been developed to prepare oxide powders in order to improve their quality (purity, chemical homogeneity, etc..) and also their reactivity. Herein, the "amorphous citrate" route was preferentially used to synthesize oxide precursors. It involves citric acid as a complexing agent and is not, strictly speaking, a classic sol-gel process in the sense that no gel is formed by a metal-oxygen-metal network, but rather from metal-organic complexes which yield after calcination ultrafine reactive powders with a good chemical homogeneity. The principle consists in preparing a solution containing the complexed cations in the desired ratios, afterwards the solution is evaporated until formation of a viscous liquid. A progressive heating to 250°C follows, then a calcination up to 500-600°C to get rid of all organics. The citrate route can be generalized to a large number of compositions which cannot be obtained via traditional methods. The use of citric acid presents several advantages such as formation of very stable solutions of more or less complex stoichiometries. The cation stoichiometry being the same in the solution and in the powdered residue, a priori any designed composition may be easily and rapidly obtained by this aqueous process.

Nitridation reactions were carried out in an alumina boat containing the oxide precursor powder placed inside an alumina tube through which ammonia gas flowed at a rate of 30-40 l.h -1 . The temperature was raised in the range 500-1000°C, depending on the precursor and system studied.

XRD powder patterns were recorded using a Philips (PW 3710 X'PERT) diffractometer operating with Cu K  radiation ( = 1.5418 Å). Nitrogen and oxygen contents were determined with a LECO TC-436 analyser using the inert gas fusion method, nitrogen as N 2 by thermal conductivity measurement, oxygen as CO 2 by infrared detection. Diffuse reflectance spectra were collected using a UV-Visible Varian Cary 100 spectrometer in which the absolute reflectance of the sample is compared with a "spectralon" standard (Labsphere Cie).

Results and discussion

Our first approach of solid solutions has concerned the systems R-(W,Mo)-(O,N), more particularly the rare-earth tungstates R 6 WO 12 . Y 6 WO 12 is a representative term of this anion-defect fluorite-type series (Y 6 WO 12  2 or Y 3.43 W 0.57 O 6.85  1.15 ), giving a white powder absorbing at 3.65 eV.

 Thermal ammonolysis of Y 6 WO 12 between 500 and 800°C for different times forms a series of oxynitrides with a defect fluorite-type structure [START_REF] Diot | Macaudière: Rare-earth and tungsten oxynitrides with a defect fluorite-type structure as new pigments[END_REF][START_REF] Tessier | Ternary and higher order rare-earth nitride materials: synthesis and characterization of ionic-covalent oxynitride powders[END_REF]. The progressive anionic substitution of nitrogen for oxygen according to: 2N 3-+ 1 = 3O 2- reaches the nitrogen-rich limit composition Y 3.43 W 0.57 O 4.3 N 1.7  2 . It is accompanied by a continuous shift of the absorption edge from 3.65 to 2.36 eV, as a consequence of an increased covalent character that leads to a decrease in the optical bandgap value (Table 1, Fig. 1). Hence, it is possible, by adjusting the nitrogen content, to obtain an oxynitride composition cutting off exactly at 3.1 eV. The nitrogen incorporation causes a decrease in the maximum of diffuse reflectance. Besides, less steep than for the starting oxide, the slope of the absorption edge does not change with nitrogen content. However, the spectral selectivity is not good enough for anti-UV application.  As molybdenum is more electronegative than tungsten, the energy level of the 4d(Mo) orbitals is lower than that of the 5d(W) orbitals. Consequently, a progressive cationic substitution of Mo for W induces a shift of the absorption edge, from 3.65 eV for Y 6 WO 12 to 2.80 eV for Y 6 MoO 12 (Fig. 2). The non linear variation law of the optical bandgap along the solid solution suggests that the 4d(Mo) and 5d(W) levels form two bands relatively distinct from each other [START_REF] Larcher | Recherche de nouveaux matériaux anti-UV[END_REF]. We have succeeded in preparing an intermediate composition cutting off precisely at 3.1 eV. Moreover, after powder annealing, a steep absorption edge could be observed, giving a spectral selectivity of 0.49 nm, comparable to 0.29 nm for TiO 2 nanorutile). Consequently, Y 6 W 1-x Mo x O 12 compositions with x ~ 0.8-0.85 seem to be attractive materials for anti-UV purposes, as they have in addition a low refractive index (n ~

2) [START_REF] Marchand | Nouveaux tungstates de terres rares, leur procédé de préparation et leur utilisation comme agent de protection contre les rayons ultraviolets[END_REF]. Electron energy loss spectrometry (EELS) measurements, which are capable of distinguishing absorption from diffusion, are now in progress. This study is funded, in particular, by the French Ministry of Trade and Industry under the acronym AUVIB [START_REF] Auvib | Absorbeurs UV Inorganiques pour le Bois[END_REF]. The project aims at emerging new varnish and wood stain formulations containing stable UV-absorbers, in order to prevent or limit wood photodegradation.

Our second approach of solid solutions has concerned domains ranging between two perovskite-type limits: a white oxide ATiO 3 (A = Ca, Sr, Ba), absorbing in the UV, and a strongly coloured oxynitride RTiO 2 N (R = La, for ex.) [START_REF] Marchand | Handbook on the Physics and Chemistry of Rare Earths[END_REF]. Applying the cross-substitution principle:

R 3+ + N 3-= A 2+ + O 2-
we have attempted to synthesize intermediate compositions showing an absorption edge located between those of the two limits of the solid solution domain.

As illustrated in Fig. 3, the powder colour acts as an excellent first indicator to estimate the width of the bandgap and, consequently, the corresponding absorption edge position. 2). The absorption properties change progressively when decreasing x, the powder colour going from red-brown (x = 1: LaTiO 2 N) to yellow (x = 0.10). The yellow colour of the Sr 0.90 La 0.10 TiO 2.90 N 0.10 composition corresponds to an absorption edge located at around 500 nm, a lower x value is therefore necessary. However, the powder corresponding to x = 0.05 shows after ammonia treatment an undesired yellow-green colour, indicative of a black component and of a reduction phenomenon of Ti IV , which is no more totally stabilized under the used reducing ammonia synthesis conditions. This would indicate that the inductive effect of strontium ions is weaker than that of lanthanum ions. This study, still in progress, focuses now on the different possibilities to modify the alkalineearth/rare-earth couple, and also the (alkaline earth + rare-earth)/titanium ratio, in order to enhance the inductive effect and thus stabilize Ti IV during thermal ammonolysis.

Conclusion

The amorphous citrate route proves its efficiency to prepare multinary oxide compositions. By this way, we have succeeded in modifying optical absorption properties by adjusting both oxide and oxynitride chemical compositions. Potential oxide UV-absorbers have been isolated. Preliminary results of a study of their photocatalytic activity (phenol photodegradation test) [START_REF] Colbeau-Justin | [END_REF] are encouraging. Another aspect concerns optimization of the performances of oxynitride compositions, especially the steepness of the absorption edge. In all cases, the powder morphology is an important parameter to be controlled before envisaging an industrial opening for these oxide and oxynitride materials.

Fig. 2 .

 2 Fig. 2. Diffuse reflectance spectra of Y 6 W 1-x Mo x O 12 compositions, compared with TiO 2 nanorutile.

Fig. 3 .

 3 Fig. 3. Colour of the oxynitride powder as a function of the bandgap energy value.

Table 1 .

 1 Y 6 WO 12-3x N 2y oxynitrides synthesized from an oxide precursor issuing from the citrate route. Diffuse reflectance spectra of several Y 3.43 W 0.57 O 6.86-3x N 2x (0 ≤ x ≤ 0.52) compositions.

	NH 3 treatment	wt% N	Fluorite-type formulation	Colour
	500°C/5h	0.4	Y 3.43 W 0.57 O 6.63 N 0.17	white
	550°C/5h	1.2	Y 3.43 W 0.57 O 6.23 N 0.40	yellow
	600°C/60h	2	Y 3.43 W 0.57 O 5.77 N 0.74	yellow
	800°C/15h	2.9	Y 3.43 W 0.57 O 5.31 N 1.03	mustard yellow

Table 2 :

 2 Colour of nitrided powders in the (Sr,La)-Ti-(O,N) system.

	Cationic composition	Observed color
	LaTi	red brown
	Sr 0.25 La 0.75 Ti	red
	Sr 0.50 La 0.50 Ti	orange
	Sr 0.75 La 0.25 Ti	yellow orange
	Sr 0.90 La 0.10 Ti	yellow
	Sr 0.95 La 0.05 Ti	yellow green
	SrTi	white


Patrick Pujo
email: patrick.pujo@univ-amu.fr

Ilham El Khabous
email: ilham.elkhabous@gmail.com

Fouzia Ounnar
email: fouzia.ounnar@univ-amu.fr

Experimental assessment of the productivity improvement when using U-shaped production cells with variable takt time

Keywords: U-shaped production cell, straight production line, production variable throughput, linear programming for optimization, assembly line balancing

Experimental assessment of the productivity improvement when using U-shaped production cells with variable takt time

Purpose -The aims of this research is to discuss the benefits of U layout for production cell operating in variable takt time. Different experiments were conducted using benchmarks, in order to highlight the performance gap between a straight cell and a U-cell.

Design/methodology/approach -The implementation of the production cell, in U-shaped or in straight line, is optimized through linear programming based on the number of operators. The two corresponding programs, in Mosel language, use the same approach, in order to not introduce defects in the comparison of results. The study used our own datasets and well known academic benchmarks.

Findings -A comparison between the obtained takt times, with equivalent operating conditions, in U-cell and straight cell was conducted. A significant increase of the production rate was observed. This increase has often exceeded 10 per cent to reach 32 per cent. All the experiments show that, with the same number of operators, a cell in a U layout is always at least as efficient, in terms of reachable production rates, than an equivalent cell in a linear layout. 96 per cent give an improvement of production rate. Moreover, the dispersion of the U-cell results is weaker, which suggests that the U-layout gives, in more robust manner, better performances.

Research limitations/implications -The results were obtained through a study of various academic benchmarks. The results must be validated on industrial situations.

Practical implications -This paper will be very useful for researchers and practitioners in order to understand lean implementations and its derived benefits. This paper will allow them evaluating and analysing the expected benefits of the implementation of the production cell in Ushaped (operating in variable takt time).

Originality/value -U-Cells constitute an appropriate solution for a layout of any kind of production cells with variable structure. When facing a significant variation in the demand, the response consists in modulating the number of operators assigned to the cell. This study addresses jointly the problem of U-cells layout and the operation in variable takt time.

Introduction

Recent developments towards lean management [START_REF] Holweg | The genealogy of lean production[END_REF] reflect an increasing trend by companies to use this approach. This is often found in the literature in the form of testimonials: the effectiveness of lean management methods and tools are typically expressed through case studies [START_REF] Detty | Quantifying benefits of conversion to lean manufacturing with discrete event simulation: a case study[END_REF], [START_REF] Motwani | A business process change framework for examining lean manufacturing: a case study[END_REF], [START_REF] Modarress | Kaizen costing for lean manufacturing: a case study[END_REF], [START_REF] Kumar | Implementing the Lean Sigma framework in an Indian SME: a case study[END_REF]. For example, tools such as value stream mapping (VSM) have been the subject of evaluation 'before & after' comparison, used alone [START_REF] Singh | Lean implementation and its benefits to production industry[END_REF] or in combination with other techniques in the case of production lines [START_REF] Alvarez | Redesigning an assembly line through lean manufacturing tools[END_REF]. Based on value stream mapping and milkrun, using lean metrics, such as dock-to-dock time and lean rate, these works analyse the internal materials flow and its improvement.

Other tools, such as U-shaped cells producing or takt time variable, have not been subjected to comparison 'before & after'. U-shaped production cells (figure1(b)) emerged in the Western industrialized world with the introduction of the Just in Time philosophy [START_REF] Monden | Toyota Production System[END_REF] and are currently associated with the Lean Production approach [START_REF] Shah | Lean manufacturing: context, practice bundles, and performance[END_REF].

The majority of the production lines not identified as falling within one of these two approaches are organized linearly (figure 1(a)). The optimization of their performance is obtained using load balancing techniques which are dealt with in numerous studies [START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF].

Please insert here Figure 1. Straight production cell (a), U-shaped production cell (b) U-Cells offer many advantages over traditional linear arrangements, in terms of flexibility, standardization, decrease in work in progress and the maintaining of high quality products. They also provide improvements in productivity.

In the case of high variations in production volumes, it is easy to adjust the number of operators to adapt the production rate to the instantaneous demand (takt time). Through variable takt time, the production lines enable new strategies, so that lean product realisation can become a reality for variable productions [START_REF] Weston | Next Generation Manufacturing Systems[END_REF]. This is common in U-Cells, but can also be achieved in straight production lines. It becomes then interesting to compare both layout approaches and evaluate their respective performances through experiments using benchmarks.

After defining the basic principles of the operation of a variable rate production line, a few fundamentals regarding U-Cells are recalled. Then, the issue of the layout optimization for such systems is addressed and linear programming models are developed aiming to obtain the best possible balancing of the operator tasks.

These models are then implemented into a solver and tested on academic benchmarks which are set to be the most impartial and the least specific possible.

Finally, the results of these experiments are analysed and discussed.

Choice of a layout architecture for a production cell

Cellular production: some reminders

Production lines were deployed at the beginning of the 20th century, after Henry Ford installed the first assembly line in 1913 [START_REF] Sward | The legend of Henry Ford[END_REF]. After being used for a century in mass manufacturing, this concept is still evolving today [START_REF] Hirst | Flexible specialization versus post-Fordism: theory, evidence and policy implications[END_REF], being revisited [START_REF] Bartholdi | A production line that balances itself[END_REF]. This history, combined with a pragmatic quest for efficiency [START_REF] Magee | How Toyota became #1: leadership lessons from the world's greatest car company[END_REF], led the Japanese industry to develop the idea of production cells: a production cell is a small production unit which specializes in the manufacturing of a family of products or parts. This type of organization should offer greater usage flexibility, better performances in terms of quality, cost and time, as well as a fluidization of flows combined with reduction in the volume of work in progress and greater responsiveness.

In fact, a cell is a small production line whose conditions of use can vary according to needs: change of rate, of products, of the number of machines, etc. The line architecture is related either to the grouping together of all the processes related to the fabrication of a family of products, or to the grouping resulting from a frequent series of operations common to the manufacturing of different types of products.

When this approach is generalized to the whole production system, we can refer to it as cellular manufacturing [START_REF] Irani | Handbook of cellular manufacturing Systems[END_REF]: the variety of products to be manufactured, their estimated volumes and the fluctuations in demand will determine the set of cells which will need to be formed and designed. Numerous studies about cellular production focus on the characterization of this system using different approaches such as group technology [START_REF] Kulak | A complete cellular manufacturing system design methodology based on axiomatic design principles[END_REF], integer linear programming (Arifar and Ismail, 2009), metaheuristics of genetic algorithm type [START_REF] Yin | Similarity coefficient methods applied to the cell formation problem: A taxonomy and review[END_REF] or ant colony optimization algorithms [START_REF] Sabuncuoglu | Ant colony optimization for the single model U-type assembly line balancing problem[END_REF].

Another emerging problem is related to the dynamics of the control of this set of cells [START_REF] Wenming | Cell staffing and standardized work design in Chaku-Chaku production lines using a hybrid optimization algorithm[END_REF]. For instance, if we consider a workshop where many neighbouring cells respond to varying and desynchronized customer loads, it is possible to adjust the number of operators per cell and transfer operators within cells (communicating vessels) so that the workshop is fully controlled and maximum productivity is achieved, preserving responsiveness and dynamics whatever the variation of customer demands and providing a high potential increase in production rate. For example, a cell functioning with n operators (nominal operation) can be adjusted to n-1 or n+1, but it is possible to make it function even faster, with up to n+k operators, the additional operators being taken from neighbouring cells or even make it function very slowly (with only one operator) or stop it temporarily.

In any case, the search for the optimal cell productivity requires the balancing of the production loads: the key is to set up the cell so it can produce at the lowest cost possible while complying with the technical constraints.

Balancing of a production cell

This issue was first raised in automotive assembly lines, and is now referred to in the scientific literature as SALBP (Simple Assembly Line Balancing Problem) [START_REF] Baybars | A survey of exact algorithms for the simple assembly line balancing[END_REF]. A typology of the resolution methods for this type of problems was suggested by [START_REF] Boysen | A classification of assembly line balancing problems[END_REF]. Generally speaking, these methods consist in splitting the task as effectively as possible between the operators so that the workload is equally distributed and the production rate is optimal.

The general features of this problem are the following:  a single operation cannot be shared between different operators;

 the operations are partially ordered;  all operations must be executed;  durations of operations must not depend either on the operator or on previous operations;

 an operator can carry out any operation;

 any operation can be executed by any operator;

 the products pass one by one in front of operators in a given order;  the operators themselves are positioned in a given order;  one single product is produced per line.

Operations have to be split between operators so that waiting times are minimal.

The waiting time of each operator is determined by the difference between the cycle time of the line (the load of the busiest operator) and the actual load (the sum of all the times of operation assumed by the operator). The solution in which the operator waiting time is the shortest provides the lowest cost.

Two cases are possible. First, if the operators are assigned to more or less versatile workstations (for example, in the case of an assembly workshop), the balancing will consist in assigning them a maximum number of operations. Then, if the operators have to move between different workstations, the balancing will consist in assigning them various consecutive workstations, with a displacement cycle between these workstations.

Even though the mathematical formulation of these two cases is equivalent, our interest will be mainly focused on the second case. The displacement of operators between different workstations is called Chaku-Chaku [START_REF] Wenming | Cell staffing and standardized work design in Chaku-Chaku production lines using a hybrid optimization algorithm[END_REF].

Let us note that the so called "caravan system" [START_REF] Spengler | Zum Einsatz von Chaku-Chaku-Systemen in der Montage konsumentennaher Erzeugnisse -eine Fallstudie bei Rahmenauftragsfertigung[END_REF] It is also important to detail the nature of tasks Ti which refer to the operations carried out by operators Ok on the different workstations Pi. We should distinguish between durations Doi of tasks Ti and displacement times between workstations Ddxi. The latter are usually small compared to the cumulative durations of tasks performed by each workstation and they can be either neglected or integrated to the corresponding durations: Doi includes Ddxi. Besides the displacement time, the duration includes product loading/unloading tasks on the workstation, possible product control tasks, and either the execution of this task -in the case of a purely manual task -or the initiation of the execution of the task followed by a short period of monitoring to visually check whether the task has started correctly -in the case of an automated workstation. This refers to the apparent time component Di of the operation on workstation Pi.

For automated tasks, the continuation and completion of the task progress in masked time, on a period called Dmi. Note that Doi = Di + Dmi. To each workstation Pi corresponds a chronogram of all the tasks of the operation [START_REF] Monden | Toyota Production System[END_REF] -For a given operator k, to which nk operations i have been assigned:

∑ 𝐷 𝑖 𝑗 𝑛 𝑘 𝑗=1 < Tk (1)
-For each given operation i:

Di + Dmi < Tk (2)
-If T is the cycle time of the cell and np the number of operators, then:

 k  [1, np], max(Tk) = T (3)
These two overall variables strongly influence the resolution method, depending on whether they are fixed or to be minimized. SALBP can be expressed in 4 ways [START_REF] Baybars | A survey of exact algorithms for the simple assembly line balancing[END_REF], [START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF]:

Please insert here Figure 3. Classification of line balancing problems SALBP-F is the feasibility study of a configuration, and not a search for optimization. SALBP-1 studies the minimization of the number of operators while respecting a given production rate. SALBP-2 aims at minimizing the cycle time so as to achieve a maximum production rate through a fixed number of operators. SALBP-E combines those two options and maximizes the efficiency of the production line

These SALBP problems deal with the balancing of straight cells. The layout of U-shaped production cells is the subject of specific studies collectively known as SULBP (Simple U-Line Balancing Problem) [START_REF] Aase | U-shaped assembly line layouts and their impact on labour productivity: An experimental study[END_REF] or UALBP (U-Assembly Line Balancing Problem) (Scholl and Klein, 1999b).

Balancing of a U-Cell

For U-shaped cells, the classification of the resolution methods is similar [START_REF] Kriengkorakot | The U-line Assembly Line Balancing Problem[END_REF]. [START_REF] Boysen | A classification of assembly line balancing problems[END_REF] identifies a number of jobs, denoted by [ |u| ].

Regarding SULBP-1, [START_REF] Miltenburg | Balancing U-lines in a multiple U-line facility[END_REF] proposed first a dynamic programming approach, then an integer linear programming approach (Miltenburg, 2001b). More recently, [START_REF] Sabuncuoglu | Ant colony optimization for the single model U-type assembly line balancing problem[END_REF] suggested a resolution using an ant colony metaheuristic.

Other studies relate to SULBP-E: for example, [START_REF] Chiang | The stochastic U-line balancing problem: A heuristic procedure[END_REF]) used a heuristic procedure; [START_REF] Kara | Binary fuzzy goal programming approach to single model straight and U-shaped assembly line balancing[END_REF] applied a multi-objective fuzzy integer programming method.

Finally, a number of studies address the Mixed Model U-Line Balancing and Scheduling Problem (MMULBS) regarding the production of a variety of products in a U-Cell, simultaneously considering balancing and scheduling [START_REF] Miltenburg | Balancing and scheduling mixed-model U-shaped production lines[END_REF], [START_REF] Kim | An endosymbiotic evolutionary algorithm for the integration of balancing and sequencing in mixed-model U-lines[END_REF], [START_REF] Kazemi | A novel two-stage genetic algorithm for a mixed-model U-line balancing problem with duplicated tasks[END_REF].

Performance comparison between the balancing of a straight cell and a Ushaped cell

Numerous studies were undertaken simultaneously on both approaches [START_REF] Ağpak | A chance-constrained approach to stochastic line balancing problem[END_REF]Gökçen, 2007, Kara et al., 2009), but focusing on the way the production cells were organized and thus on the cell architecture generation algorithm, without comparing their performances with equivalent datasets.

As our study is performed from a perspective of variable takt time, we can compare both types of layout by changing the number of operators and comparing the respective production rates obtained. For that, we should move towards SALBP-2 and SULBP-2 approaches.

Modelling of the SxLBP-2 balancing problem for a production cell

In this section, the implemented linear programming models are described. The two models have to be homogeneous in terms of structure and design in order to avoid any methodological biases. As there is comparison, we chose to consider one of the layouts as a reference point, in terms of modelling, implementation and results analysis.

Based upon the literature review, we chose to start from an existing SALBP-2 approach and extend it to a U-Cell.

Modelling of the SALBP-2 problem

The SALBP-2 model described here results from [START_REF] Guéret | Applications of optimization with Xpress-MP Dash Optimization Ltd, Translated and revised[END_REF]. This model is relatively simple, gives excellent results and its implementation is already using Xpress solver. The following notations are considered:

 nt: number of operations,

 Di: duration of operation i (i  [1…nt]),
 np: number of operators,

 na: number of arcs (sequence of operations). An arc (u, v) links operation u to operation v if u is the immediate predecessor of v.

In terms of modelling, xik binary variables will be used to place the operations:

 xik = 1 if and only if operation i, is assigned to operator k (see (4) & (8)).
For all valid assignments, in accordance with (4) and (8), it is also necessary to respect precedence constraints: for any arc (u, v), the rank of operation u must be lower than or equal to that of operation v (see ( 5)).

Finally, a real variable T≥0 is used (see (9)) for the cycle time of the cell, which according to (6) is equal to the cycle time of the busiest operator. The objective function consists in minimizing T (see ( 7)).

The formulation of the SALBP-2 problem is as follow:

 i=1…𝑛 𝑡 : ∑ 𝑥 𝑖𝑘 𝑛 𝑝 𝑘=1 = 1 (4)  (u,v)[1…𝑛 𝑡 ]: ∑ 𝑘. 𝑥 𝑢𝑘 𝑛 𝑝 𝑘=1 ≤ ∑ 𝑘. 𝑥 𝑣𝑘 𝑛 𝑝 𝑘=1
(5)

 k=1…𝑛 𝑝 : ∑ 𝐷 𝑖 . 𝑥 𝑖𝑘 𝑛 𝑡 𝑖=1 ≤ T (6) Min z = T (7)  i=1…𝑛 𝑡 ,  k=1…𝑛 𝑝 : 𝑥 𝑖𝑘  0,1 (8)
T ≥ 0 (9)

Under these conditions, the first operator handles the cell inputs while operator np manages the outputs.

Extending the modelling of SALBP-2 to SULBP-2

The U shaped-cell layout involves treating simultaneously incoming and outgoing branches of the U-cell. Let j be the branch index, with j = 1 for the incoming branch and j = 2 for the outgoing branch.

Modelling the SULBP-2 problem thus requires the binary variables to adjust to this notion of treating simultaneously the two branches to position the operations.

 xijk = 1 if and only if operation i is assigned to operator k on branch j (see ( 10) & ( 15)).

For all valid assignment, in accordance with (10) & (15), it is also necessary to respect precedence constraints: for any arc (u, v), the rank of operation u must be lower than or equal to that of operation v. However, the calculation of this rank can differ according to the branch. In branch 1, we are back to the previous case (see (11)). In branch 2, the precedence must be reversed (see ( 12)).

The principle of minimization of the busiest operator cycle time remains unchanged (see ( 13)), and takes into consideration the two branches, ( 14) & ( 16)).

The formulation of the SULBP-2 problem is as follows:

 i=1…𝑛 𝑡 : ∑ ∑ 𝑥 𝑖𝑘𝑗 𝑛 𝑝 𝑘=1 2 𝑗=1 = 1 (10)  (u,v)[1…𝑛 𝑡 ]: ∑ 𝑘. 𝑥 𝑢𝑘1 𝑛 𝑝 𝑘=1 ≤ ∑ 𝑘. 𝑥 𝑣𝑘1 𝑛 𝑝 𝑘=1 (11)  (u,v)[1…𝑛 𝑡 ]: ∑ 𝑘. 𝑥 𝑣𝑘2 𝑛 𝑝 𝑘=1 ≤ ∑ 𝑘. 𝑥 𝑢𝑘2 𝑛 𝑝 𝑘=1 (12)  k =1…𝑛 𝑝 : ∑ ∑ 𝐷 𝑖 . 𝑥 𝑖𝑘𝑗 𝑛 𝑡 𝑖=1 2 𝑗=1 ≤ T (13)
Min z = T ( 14)

 i = 1…𝑛 𝑡 ,  k = 1…𝑛 𝑝 ,  j = 1, 2: 𝑥 𝑖𝑘𝑗  0,1 (15)
T ≥ 0 (16)

In this formulation of the SULBP-2 problem, operator np manages the inputs and outputs of the cell. Operator 1 is situated in the heart of the cell.

Implementation of the models

The implementation of these models is performed using Mosel language and Xpress solver, available under FICO optimization environment (www.fico.com).

Thereafter, only the programming of the SULBP-2 problem (figure 4) will be described, that of the SALBP-2 problem having already been described in [START_REF] Guéret | Applications of optimization with Xpress-MP Dash Optimization Ltd, Translated and revised[END_REF].

The binary assignment variable x is represented by a 3D table (nt, np, 2), the third dimension being used to separately treat the two branches of the U.

There are nt constraints 'OneOperationPerOperator(i)' which express the fact that an operation i cannot be shared between several operators.

There are np constraints 'Cycle (m)' which express the fact that the sum of operations assigned to an operator m is lower than or equal to the cycle time of the U-Cell.

Finally, the constraints 'Prec_branch_input(a)' and 'Prec_branch_output(a)' allow the precedence constraints between operations to be taken into consideration for each sequence of operations and whatever the assignment of an operation on one or the other of the U-Cell branches.

The number of operators is defined by the declared variation range of the variable 'OPERATORS'.

Please insert here Figure 4. SULBP-2 problem programmed in Mosel language

Experimentations with simple linear precedence graph

Establishment of the dataset

The dataset presented in Table 1 is an academic benchmark which was precisely elaborated for the specific purpose of controlling the variability of a number of parameters, in order to establish a performance comparison between U-sha ped and straight cells.

Therefore, all the 34 tested cases were built by altering only the duration of the tasks, according to a certain number of invariants: number of operations ( 17workstations), total task duration (303 Time Units), duration of the shortest task (5 TU), duration of the longest task (60 TU) and purely linear flow: the output flow of operation u goes directly into operation u+1, which induces a linear precedence graph. The above mentioned models would naturally support other values for these invariants: this is to test different configurations of task durations with a variable number of operators, on linear or U type layouts.

Please insert here Table 1. Description of the dataset: the 34 studied cases.

Experimentation results

Thereafter, only the result of case n°1 (figure 5) will be detailed. For x operators, the distribution of the 17 operations is given, for SALBP-2 (Lx) and SULBP-2 (Ux), as well as the cycle durations for each operator and the takt time of the associated cell. Each line also indicates the Chaku-Chaku of the operator concerned.

Please insert here Figure 5. Case n°1 results

The 34 studied cases are summarized in Table 2. For each Lx & Ux configuration, the cycle durations and the improvement percentage are given.

We observe that for one single operator, the results are the same in a straight or U-shaped cell. This is a trivial result. This experiment is thus not significant.

For 2 or 3 operators, the results are equivalent in rare cases: the advantage generally is in favour of U-Cells.

Please insert here Table 2. Summary of the experimentation results (in bold:

improvements greater than 10%)

Similarly, for a higher number of operators, optimization is no longer efficient.

This can be easily explained since the cycle duration is bounded by the bottleneck due to the longest task assigned to the cell: Lx and Ux results converge, either for 7 operators, or starting from 6 operators. These cases are not significant. Only case 32 achieves a slight performance advantage with 7 operators, advantage which is lost with 8 operators.

Between these two extreme cases (123 significant experiments), the results are always in favour of a U-shaped layout, very often with improvements in cycle time exceeding 10%.

Based on these significant experiments, 117 provide an improvement of the performance (95%). This result is substantial.

Analysis of experimentation results

Table 3 shows a summary of the previous results. An average 8% improvement of the maximum possible production rate is obtained, with peaks up to 25%. Figure 6 displays the scattering of the results.

Please insert here Table 3 We note that the scattering of the results is much greater for a linear layout than for a Ushaped cell, which demonstrates the robustness of the U-shaped layout approach in the case of a linear precedence graph.

Experimentations with a complex precedence graph

In many manufacturing systems, flows are not always fully ordered: tasks are partially organized by precedence relations. For instance, if there are jumps between operations, or if operations can be performed in parallel or independently, these anteriority constraints between tasks can be modelled by a precedence graph, in which tasks are the nodes and precedence constraints the arcs. In the previous models, the ANT table (Figure 4) makes it possible to present each precedence relationship as an elementary flow between two operations.

The whole is then a network of precedence constraints.

We will examine the constraints induced by these considerations when operating with a variable takt time.

Problem of implementation with a complex precedence graph

The optimization models presented in the 3 rd section for the resolution of SALBP-2 and SULBP-2 provide excellent solutions regardless of the precedence graph. Nevertheless, they cause a major problem in the case of an operation with takt time varying with regular frequency: for instance if we wish to adapt on a daily basis the operation to the customer's demand. Indeed, when the number of operators is changed the order in which the tasks are arranged on the production line can be changed. However, it is unconceivable to move resources and workstations at each takt time variation. This of course does not occur when the precedence graph is strictly linear, as seen in paragraph 4.

To overcome this difficulty, we must linearize the complex precedence graph in order to identify the most appropriate linear precedence graph that would be able to support the different chaku-chaku configurations. It is thus necessary to identify the linear precedence graph without any back-tracking, regardless of the chaku-chaku configuration applied. So as to realize this linearization while respecting this goal, we propose a twostep process. First of all, a breakdown into chronological levels uses the 'Activity on node' approach applied in project management. Then, the total moment of From-to chart is minimized at each level.

Methodology for linearizing a complex precedence graph

The complex 'Activity on node' diagram is first broken down into levels [START_REF] Elmaghraby | Activity Networks: Project Planning and Control by Network Models[END_REF]. The principle is the one used in 'precedence diagramming method' approaches [START_REF] Wiest | Precedence diagramming method: some unusual characteristics and their implications for project managers[END_REF], [START_REF] Tavares | A multi-stage non-deterministic model for project scheduling under resources constraints[END_REF]. The first level contains activities with no antecedents.

These activities are removed from the activity list, and we start the same process again with the remaining activities to get to the next level [START_REF] Crandall | Project planning with precedence lead/lag factors[END_REF]. This allows us to structure the precedence graph with a global order.

Then, an order has to be found within activities on the same level. To do so, the notion of total moment Z of From-to chart is applied [START_REF] Esgin | A Hybrid Approach for Process Mining: Using From-to Chart Arranged by Genetic Algorithms[END_REF]:

𝑀𝑖𝑛 𝑍 = ∑ ∑ 𝑓 𝑖,𝑗 |𝑗 -𝑖| 𝑁 𝑗=1 𝑁 𝑖=1 (17)
where fi,j = 1 if a precedence relation exists between the 2 activities, and fi,j = 0 otherwise.

From-to chart is the relationship flow matrix. It gives the dominance of the flows between resources and indicates the closeness of resources in a facility layout design [START_REF] Mahdavi | A flow matrix-based heuristic algorithm for cell formation and layout design in cellular manufacturing system[END_REF]. From-to chart could be classified as loose and tight [START_REF] Moore | The zone of compromise for evaluating lay-out arrangements[END_REF], depending on whether it contains respectively a few and a large number of strong precedence constraints. In our case, the From-to chart contains many entries with zero values: this corresponds to the case of flow line layout [START_REF] Hassan | Machine layout problem in modern manufacturing facilities[END_REF].

To minimize this moment, we proceed to activity permutations at the same level. This allows us to find a local order which remains compatible with the previously established global order.

The linear sequence which results from this treatment can be used as a reference sequence to apply SALBP-2 and SULPB-2. Other sophisticated mathematical techniques may also be applied, but what interests us here is the possibility to compare SALBP-2 and SULBP-2 on equal and logical datasets. Thanks to this approach, we can consider that all cases with complex precedence graphs can be transformed in cases with simple linear precedence graphs. We can then use the approach mentioned in sections 3 & 4.

Experimentations with well-known benchmarks

More than 250 benchmark instances are presented in (Scholl and Klein, 1999a). A complete description of the datasets can be downloaded from the web at: http://alb.mansci.de/files/uploads/UALBP-1 data sets.zip It seemed interesting to us to test some of these instances to verify our results. We considered 12 instances containing complex precedence graphs: 4 small-sized and 8 large-sized datasets.

Small-sized benchmarks

The first 4 instances can be considered small since the number of tasks as well as the number of required operators is around 10. A description of these instances is given in Table 4.

Please insert here Table 4. Description of small-sized instances Table 4 gives an idea on the distribution of task durations and indicates the linear sequence equivalent to the initial complex precedence graph.

The results are shown in Table 5.

Please insert here Table 5. Results for small-sized instances The performance is improved in 14 significant experiments among 16 (87.5%).

Large-sized benchmarks

The same principle was applied to large-sized datasets. Table 6 describes the datasets.

Please insert here Table 6. Description of large-sized datasets.

The results are given in Table 7 Please insert here Table 7. Results for large-sized datasets.

Analysis of the results

This series of tests confirms the previous observations while leading to new observations. For small-sized instances, the results confirm the performances achieved, similar to those presented in section 4. We will observe 32% performance improvement for the Jackson instance with 3 operators.

For large-sized instances, improvements are small for low rate productions (number of operators equal to two or three). On the other hand, for higher production rates (lower takt time), improvements are substantial, easily exceeding 10%. We also note that the maximum production rate is more quickly reached with a U layout, since the maximum rate limit, defined by the longest task duration, is reached with fewer operators: out of the 8 cases tested, we note that for ['Heskia'], the highest rate is reached for a U layout with 11 operators while 12 would be needed in a straight line, for ['Gunther ', Lutzl' & 'Mitchell'], the maximum rate is reached with 2 fewer operators, and for 3 others ['Buxey', 'Kilbrid' & 'Sauvey30'], with 3 fewer operators, i.e. almost 1/5 of the payroll. This is a meaningful result.

Conclusion

All the experiments presented here show unequivocally that, with the same number of operators, a cell in a U layout is always at least as efficient, in terms of reachable production rates, than an equivalent cell in a linear layout. In summary, in order to produce more with similar loads, U-Cell layouts are needed. This was shown through our experiments and also with academic benchmarks.

Based on all tested cases (230 significant experiments), 96% of cases give an improvement of performance. This remains of course to be validated on industrial situations.

The productivity gains observed are substantial. Furthermore, the performance robustness obtained makes it possible to validate this approach in the case of a manufacturing production subject to wide variations. Indeed, regardless of the rate requested by the downstream production, the ad hoc configuration of the U layout will be the most suitable for obtaining an optimal performance. This shows the great interest of the U-shaped cell layout approach.

We conclude that U cells constitute a solution for implementing any relevant and effective production system with variable structure. It is an essential tool for a lean manufacturing approach. ,2,3,4,5,7,6,8,9,10, C 3,1,2,4,5,6,7,8,9,10,11 1,2,7,3,6,9,12,25,26,4,27,10,5,14,15,8,13,19,11,16,21,17,18,22,20,23,24,28,29 Gunther 35 40 13,8 1 483 12, 09 C 17,1,2,5,10,12,3,6,4,7,8,9,11,14,18,13,15,19,16,20,21,22,25,30, 23,26,31,24,32,27,33,28,34,29,35 Heskia 28 108 36,57 1 1024 33, 09 C 2,1,3,4,5,6,8,17,19,21,22,23,24,26,7,9,20,25,27,10,18,11,12,13, 15,14,16,28 Kilbrid 45 55 12,27 3 552 9, 65 C 1,2,11,12,39,3,4,7,8,13,37,5,6,14,15,43,9,10,16,17,18,23,24,25, 29,30,31,32,19,26,27,20,33,21,34,35,36,22,28,38,40,41,42,44,4 5 Lutz1 32 1400 441,9 100 14140 244,9 C 3,4,2,1,5,6,9,8,7,10,11,12,13,15,14,16,17,18,19,20,21,22,23,24, 25,26,27,28,29,30,31,32 Mitchell 21 13 5 1 105 2,97 C 1,2,3,4,21,5,6,7,8,14,9,10,11,12,13,15,16,18,19,17,20 Roszieg 25 13 5 1 125 2,81 C

1 ,2,3,4,5,8,6,9,10,7,11,12,13,15,17,14,16,19,20,23,18,21,22,24, 25 Sauyer30 30 25 10,8 1 324 6,07 C

1 ,10,3,2,4,5,11,12,16,17,6,13,18,7,14,19,8,15,20,9,21,24,22,25, 23,26,27,28,29,30

  , which consists in moving k operators (out of phase by 2/k) along the entire cell (figure 2(a)), with an identical Chaku-Chaku, is not taken into account. We will focus on the balancing of the line between different Chaku-Chaku (one operator corresponding to one Chaku-Chaku), involving a manual transfer of parts between neighbouring Chaku-Chaku operators (figure 2(b)). Please insert here Figure 2. Chaku-Chaku in 'caravan' system (a), Chaku-Chaku with manual transfer (b)

  , called 'Standard Operation Chart' and all the chronograms show the cyclic activity of operator Ok in a Gantt diagram, called 'Standard Operation Combination Chart', If Tk is the cycle time of operator k, the following relationships must be fulfilled:

Figure 1 .Figure 2 .Figure 6 .

 126 Figure 1. Straight production cell (a), U-shaped production cell (b)

  . Average and standard deviation of cycle durations of the

	various configurations, for Lx & Ux (x  [2, 6])
	Please insert here Figure 6. Distribution of the results, for Lx & Ux (x  [2, 6])

Table list

 list

	Case n° P1 P2 configurations, for Lx & Ux (x  [2, 6]) P3 P4 P5 P6 P7 P8	P9	P10 P11 P12 P13 P14 P15 P16 P17
	1	5	23	12	10		15		5		30	8	12	10		20	32	9	24	13	60	15
	2	1	9	17	60		5		14		22	50	11	29		10	20	7	10	6	9	14
	3	1	60	15	5		13		22		30	40	4	16		20	20	14	7	10	8	9
	4 5 6 7 8 9	0 0 5 5 5 1 1 1 1 Instance 8 10 15 5 60 21	5 20 20 12 22 25 number of tasks	11 18 60 34 15 5 Task duration max	10 25 35 24 9 22	Average duration	12 20 10 18 6 27	Task duration min	12 52 8 14 56 14	15 6 9 42 10 24 ∑tasks durations	13 60 50 11 16 60 Standard deviation	17 11 12 60 6 11 Lin/complex	14 8 41 16 7 9	19 22 5 8 11 13 obtained sequence 20 24 11 9 7 10 10 15 13 8 5 9	23 5 6 6 20 13	35 10 5 8 5 20	60 11 5 9 21 8
	10 11	8 5 Bowman8 9 7 1 5 Jackson	16 8 60 11	15 7	17 7	13 16	60 9,38 24 4,18	16 3 11 1	5 75 31 46	50 4,31 12 1,94	34 9	C C	10 1,2,4,3,6,5,8,7 15 11 28 1	6 19	9 5	14 39	8 7	18 9
	12	1	21	22	60		5		19		27	17	23	12		16	8	10	18	11	7	13
	13	0 5	60	51	18		39		6		12	18	11	21		6	6	7	14	7	14	8
	14	4 6	8	14	36		13		18		7	25	5	5		17	16	12	30	6	10	21
	15	2	11	8	5		60		14		30	6	12	8		25	15	36	27	7	9	10
	16	0 1	8	15	22		18		5		60	25	29	6		18	26	9	13	19	7	13
	17	0 5	20	20	5		5		30		60	20	8	5		30	40	10	5	20	10	10
	18	0 1	35	5	60		23		8		20	11	8	17		15	24	6	16	14	22	9
	19	6	5	40	30		20		10		5	7	10	22		20	39	8	5	7	10	5
	20	0 6	5	25	10		15		7		30	5	20	6		23	10	9	49	60	7	16
	21	0 5	23	14	25		11		9		7	30	13	29		8	30	9	5	60	17	8
	22	3	37	22	9		34		11		12	10	7	8		60	5	9	20	6	8	15
	23	5	12	29	35		8		7		10	30	14	9		8	60	10	7	15	23	21
	24	0 8	17	60	6		12		15		9	17	20	22		26	5	21	24	13	22	6
	25	1	6	15	17		10		60		30	24	5	20		10	16	9	12	10	18	29
	26	9	6	12	19		5		7		31	24	17	60		11	15	10	18	24	7	28
	27	2 1	25	5	32		12		15		11	60	12	10		11	35	11	5	22	15	7
	28	1	15	5	30		23		8		10	15	12	14		16	60	13	20	13	32	5
	29	5 5	20	12	9		16		5		27	8	15	10		20	28	9	28	16	60	15
	30	2 7	9	31	21		9		11		10	8	32	60		8	9	25	33	14	11	5
	31	6	5	11	52		12		32		5	7	8	9		9	9	40	20	6	6	12
	32	1	6	20	5		26		18		30	15	33	7		21	10	12	14	6	60	9
	33	0 6	17	60	5		5		51		9	13	30	13		21	20	5	11	5	25	7
	34	1 2	60	10	5		13		16		33	19	35	11		8	15	20	9	12	7	5
	5 Table 1. Description of the dataset: the 34 studied cases
	.

Table 4 .

 4 Description of small-sized instances

	Instance	x	1	2	3	4	5	6	7	8
	Bowman8	Lx	75	42	28	22	20	17
		Ux	75	38	26	21	18	17
		%	0%	10%	7%	5%	10%	0%
	Jackson	Lx	46	24	19	13	12	11	9	7
		Ux	46	24	13	12	10	9	8	7
		%	0%	0%	32%	8%	17% 18% 11%	0%
	Jaeschke	Lx	37	20	14	10
		Ux	37	19	13	10
		%	0%	5%	7%	0%
	Mansoor	Lx	185	98	68	50	45
		Ux	185	93	65	49	45
		%	0%	5%	4%	2%	0%

Table 5 .

 5 Results for small-sized instances

		2 operators		3 operators		4 operators		5 operators		6 operators
		L2	U2	%	L3	U3	%	L4	U4	%	L5	U5	%	L6	U6	%
	average	160	153	5%	112	103	8%	86.1	78.0	9%	72.4	63	13%	62.7	60	7,7%
	standard deviation	6.44	0.71		5.77	0.84		4.31	1.10		4.92	1.13		2.90	0

Table 6 .

 6 Description of large-sized datasets

	Instance	x 16 % 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 0% 3% 4% 3% 8% 10% 10% 4% 14% 10% 14% 2%
	Mitchell	Lx	105	56	37	28	23	21	18	16	16	15	13
		Ux	105	53	35	27	22	18	16	15	13	13	13
		%	0%	5%	5%	4%	4%	14%	11%	6%	19%	13%	0%
	Roszieg	Lx	125	64	43	34	29	23	21	18	16
		Ux	125	63	42	32	26	21	19	17	16
		%	0%	2%	2%	6%	10%	9%	10%	6%	0%
	Sauyer30	Lx	324	166	111	87	73	57	53	48	44	38	35	33	31	28	26	26	25
		Ux	324	162	109	82	66	55	48	42	38	35	31	29	27	25	25	25	25
		%	0%	2%	2%	6%	10%	4%	9%	13%	14%	8%	11%	12%	13%	11%	4%	4%	0%

Table 7 .

 7 Results for large-sized datasets.

Acknowledgements.

The authors thank FICO, which, through its Academic Partner Program, has kindly put the FICO TM Xpress Optimization Suite at our disposal.


Enrico Lovisari
email: enrico.lovisari@gipsa-lab.fr.

Carlos Canudas De Wit
email: carlos.canudas-de-wit@gipsa-lab.fr.

Alain Y Kibangou
email: alain.kibangou@gipsa-lab.fr.

Data fusion algorithms for Density Reconstruction in Road Transportation Networks

Keywords: Road Transportation systems, Dynamical flow network, Density reconstruction, Floating Car Data

This paper addresses the problem of density reconstruction in traffic networks with heterogeneous information sources. The network is partitioned in cells in which vehicles flow from their origin to their destination. The state of the network is represented by the densities of vehicles in each cell. Density estimation is of crucial importance in future Intelligent Transportation Systems for monitoring, control, and navigation purposes. However, deploying fixed sensors for this purpose can be very expensive. Therefore, most of fixed sensors networks are rather sparse. On the contrary, recent technologies have enormously increased the availability of relatively inexpensive Floating Car Data. A data fusion algorithm is then proposed to incorporate the two sources of information into a single observer of density of vehicles. The efficiency of the proposed algorithm is shown in a real scenario using data from the Grenoble Traffic Lab fixed sensor network and INRIX Floating Car Data on the Rocade Sud in Grenoble.

I. INTRODUCTION

The increase of the number of vehicles observed in the last decades has steered crucial roads towards a state of near saturation, exhibiting on a daily basis periods of congested traffic [START_REF] Papageorgiou | ITS and traffic management[END_REF], in turn causing severe economical and environmental losses. Built upon recent technological and theoretical advancements, Intelligent Transportation Systems (ITSs) are expected to provide robust techniques for real-time monitoring, prediction and actuation of traffic networks.

This paper studies the problem of estimating the density of vehicles in traffic networks. Such an information allows to forecast travel time and traffic evolution, therefore providing navigation suggestion to drivers and allowing traffic managers to control the network through traffic lights, ramp metering and speed limits [START_REF] Papageorgiou | Review of road traffic control strategies[END_REF], [START_REF] Papageorgiou | ALINEA: A local feedback control law for on-ramp metering[END_REF], [START_REF] Pisarski | Optimal balancing of road traffic density distributions for the cell transmission model[END_REF], [START_REF] Como | Throughput optimal distributed routing in dynamical flow networks[END_REF].

The main tool to estimate the state of the network are fixed sensors, such as induction loops or magnetometers. These devices are able to 1) count the number, 2) estimate the density, and 3) measure the average speed of the vehicles that cross a certain section of road in a period of time.

In addition, the spread of wireless devices allows new sensing and communication capabilities. In particular for the traffic application, any vehicle equipped with GPS devices can act as a probe in the traffic and provide Floating Car Data (FCD). If a non negligible fraction of vehicles acts as probe, the collected data provides an estimate of the evolution of speed in the network. Due to privacy reasons, single vehicles traces are usually not directly used, but rather aggregated as average speed of vehicles in segments of road. Advanced methodologies, such as the one employed by INRIX, ensure a very fine spatial partition of the network, with segments as short as 250 meters ([6]). Compared to fixed sensors, this technology is less precise, but since it exploits existing communication systems it is relatively less expensive and already covers all major traffic networks.

In this paper we propose an algorithm that aims to reconstruct the traffic density by fusing fixed sensors measurements and Floating Car Data. We employ a macroscopic model, partitioning the network in cells and assigning to each cell a density of vehicles. The latter evolves dynamically according to a first order mass-conservation law.

The most celebrated macroscopic traffic model is the PDE based Lighthill-Whitham and Richards (LWR) model [START_REF] Lighthill | On kinematic waves. ii. a theory of traffic flow on long crowded roads[END_REF], which is able to reproduce crucial phenomena such as traffic shock waves, and a well known discretization of the LWR-PDE is the Cell Transmission Model (CTM) [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF]. Huge efforts have been put in the last 15 years to calibrate the CTM [START_REF] Muñoz | A piecewiselinearized cell transmission model and parameter calibration methodology[END_REF] and to unveil its system-theoretical properties [START_REF] Morbidi | A new robust approach for highway traffic density estimation[END_REF]. Fusion of flow, density and speed measurements, mostly considering single vehicles traces, has also been addressed. Approaches range from signal processing techniques such as the generalized Treiber-Helbing filter [START_REF] Van Lint | A robust and efficient method for fusing heterogeneous data from traffic sensors on freeways[END_REF], to nonlinear versions of the Kalman filter [START_REF] Work | Lagrangian sensing: traffic estimation with mobile devices[END_REF], to stochastic versions of the three-detector model [START_REF] Deng | Traffic state estimation and uncertainty quantification based on heterogeneous data sources: A three detector approach[END_REF].

Our approach inherits from the CTM the cell-based topology, but we do not directly employ the resulting dynamical model. Instead, inflows and outflows are estimated on the basis of the available flow measurements only, and speed measurements are employed to compute a pseudomeasurement of the density. These quantities are the inputs for the density observer. In addition, we propose a gradient descent method to calibrate the Fundamental Diagram, and we implement the proposed solution using real fixed sensor measurements from the Grenoble Traffic Lab [START_REF] Canudas De Wit | Grenoble traffic lab: An experimental platform for advanced traffic monitoring and forecasting[END_REF], a sensing network deployed along the Rocade Sud in Grenoble, France, and speed FCD measurements provided by INRIX, one of the most well known traffic solutions companies.

The remainder of the paper is organized as follows: Section II formulates the problem and details the available measurements. Section III describes the proposed solution, while Section IV presents our numerical experiments. Finally, Section V draws the conclusions and presents several future research directions.

A. Notation

The symbols R n , R n + and R n×m denote the sets of real and positive real valued vectors of dimension n, and of real valued matrices of dimension n×m, respectively. For a finite set A, the symbol R A denotes the set of vectors indexed by elements of A, and |A| its cardinality. The transpose of A ∈ R n×m is A T . For a vector x ∈ R n , x ≥ 0 is meant component-wise. The symbol I denotes the identity matrix of suitable dimension.

A graph G is a pair (V, E) where V is called the set of nodes and E the set of edges. Edges are equipped with two functions t : E → V and h : E → V, the head and tail functions, respectively, so that e is the directed edge between t(e) and h(e). We allow for parallel edges, namely, edges having the same head and tail, but not for loops, namely, edges whose head and tail coincide. A path of length n ≥ 2 is a sequence of edges e 1 , . . . , e n that are consecutive, namely, such that h(e i ) = t ( e i+1 ) for i = 1, . . . , n -1. A path of length 1 is a path made of a single link.

II. PROBLEM FORMULATION

A Transportation Network is a graph G = (V, E), where junctions v ∈ V are interfaces between cells e ∈ E, multi-lane roads are partitioned in sequences of parallel cells. Cells that carry vehicles in and out from the network are called onramps and offramps, respectively, and each cell e lies on the path from an onramp to an offramp. Time is discrete and slotted in intervals of length T > 0. On each cell e ∈ E we denote by ρ e (t) the density of vehicles, in vehicles per km (veh/km), during the t-th time slot. The state of the network ρ(t) = [ρ e (t)] e∈E changes dynamically in time according to the first-order model

ρ e (t + 1) = ρ e (t) + 1 e (f in e (t) -f out e (t)) (1)
where e is the length of cell e, and f in e (t) and f out e (t) are the inflow and the outflow at cell e during the t-th time slot.

The fraction of vehicles on cell e that turns into the cell j, the splitting ratio from e to j, is denoted R ej ≥ 0. If j does not follow e, then R je = 0, and j∈E R ej ≤ 1, with strict inequality at offramps. By definition of splitting ratios, one has f in e = j∈E R je f out j , which, stacking inflows and outflows into vectors f in and f out , is written in matrix form

f in = R T f out (2)
where the matrix R = [R ej ] is the matrix of splitting ratios.

In this paper we assume that R is fixed and predetermined. As a Godunov scheme for the LWR-PDE model, the CTM and derived models postulate that the outflow f out e on cell e depends on the density of vehicles on e and on cells immediately following e. We adopt a different point of view and renounce to explicitly determine the flow-density relation. The sole assumption that we make is the standard

f out e = ρ e v e (ρ), ∀e ∈ E (3)
namely that the outflow from a cell is the product of density and average speed of vehicles on the same cell. While in

e j k R ej R ek FCD segment s(e) = s(j) = s(k)
Fig. 1. A stretch of road partitioned in cells and FCD segments. Splitting ratios are shown from a cell e to following cells, j and k. A FCD segment including, among others, cells e, j and k, is also shown.

the CTM the speed is a function of the density of vehicles around cell e, we will leave it unmodelled. The reason is that, as it will be clarified in the following, we assume having measurements of speed in each cell, so there is no need to model its dependence on the state of the network. We shall consider from now on the dynamics of the real system to be dictated by Eqs. ( 1)-( 3), where v e (ρ), e ∈ E, is an unmodelled function of the local state of the network.

For data fusion and estimation purposes, we make a further simplifying assumption, writing

ϕ e = ϕ e (ρ e ), ∀e ∈ E (4)
where ϕ e is the flow of vehicles at the sensor locations.

The graph of ϕ e (•) is called the Fundamental Diagram on cell e, and it is usually supposed to be a concave function with ϕ(0) = ϕ(ρ jam ) = 0, where ρ jam is the jam density (standard values for ρ jam varying from 200 to 300 veh/km).

A. Available measurements

As already mentioned, in this paper we consider a scenario in which heterogeneous measurements are available.

1) Flow and density measurements: Measurements of flows and of density are obtained via fixed sensors, which count the number (flow) and estimate the density of vehicles that, during a time slot, cross the section of road where they are placed. Without loss of generality and for sake of simplicity, we assume that a new set of flow and density measurements is available at every time slot, so we can write

ϕ m e (t) = f out e (t) + ω ϕ e (t), e ∈ E m ρ m e (t) = ρ e (t) + ω ρ e (t), e ∈ E m (5)
where E m ⊆ E is the set of cells equipped with sensors, ϕ m e (t) and ρ m e (t) are flow and density measures at time t, and ω ϕ e (t) and ω ρ e (t) are measurement errors whose stochastic properties depend on the performance of the sensor as well as on road and weather conditions, with sources of noise ranging from temporary inability to detect changes of the magnetic field, too fast or too slow vehicles, etc.

2) Speed measurements: As mentioned, vehicles equipped with GPS devices can act as proves and communicate to ITSs Floating Car Data including their position and speed (traces). Since for privacy reasons traces of private vehicles cannot be used, roads are divided into segments, let S denote the set of all segments, and data are aggregated as average speed per segment s ∈ S. Notice that the resulting measurement is an average of the speed of a fraction of vehicles in the road.

Floating Car Data are obtained through existing communication architecture, making them less expensive than fixed sensors and available on great part of the major roads. However, their precision is related to their penetration rate and therefore is not as high as fixed sensors, and moreover they cannot distinguish vehicles on different lanes. Finally, FCD are usually averaged over a relatively long period of time. For example, while the GTL sensor network provides measurement every 15 seconds, FCDs are aggregated by INRIX every minute, standard practice being around 5/10 minutes. We will assume that new speed aggregate data is available every N time instants, i.e., at times N , 2N , 3N , . . . , corresponding to the average speed in the periods [0, N -1], [N, 2N -1], etc., respectively. As such, speed measurements can be formally written as

1 v FCD e (t) = v ff e , t ∈ [0, N -1] v FCD s(e) (k), t ∈ [kN, (k + 1)N -1] (6)
where

• v ff e > 0 is the freeflow speed on cell e, namely, the speed of vehicles in low density regime;

• v FCD s(e) (k) is given by v FCD s(e) (k) = 1 N |s(e)| j∈s(e),τ ∈It v j (τ ) + ω FCD s(e) (k)
where s(e) denotes the segment of which e is one of the cells (see Fig. 1), ω FCD s(e) (k) is a measurement error whose stochastic properties depend on the performance of the sensor as well as on road and weather conditions, and

I t = {τ : t N -1 ≤ τ N < t N }.

III. A NONLINEAR OBSERVER FOR TRAFFIC NETWORKS

Eqs. ( 1)-( 2) cannot be used as they are to observe the system, as errors in flow measurement, unavoidable in real systems, are integrated in the dynamics, therefore possibly producing unbounded estimation errors. To overcome this difficulty, we consider the following Luenberger-like observer, in which ϕ m (t) = [ϕ m e (t)] e∈Em and v FCD (t) = [v FCD e (t)] e∈E are the stacked versions of the measurements,

               ρe (t + 1) = ρe (t) + 1 e ( f in e (t) -f out e (t)) +κ(ρ e (t) -ρe (t)) f in e (t) = f in e (ϕ m (t)) f out e (t) = f out e (ϕ m (t))) ρe (t) = ρe (ϕ m (t), v m (t))
∀e ∈ E (7) where

• ρe (t) is the estimate of the density on cell e at time t;

• f in e (t), f out e (t) are estimates, based on the flow measurements, of inflow and outflow in cell e at time t;

• ρe (t) is an instantaneous estimate, based on flow and speed measurements, of the density on cell e at time t; • κ is a tunable gain trading off between flow and instantaneous density estimates. Notice that it should not be too large, as to avoid meaningless negative densities.

Our aim is to design the maps

f in = { f in e } e∈E : R E m + → R E f out = { f out e } e∈E : R E m + → R E ρ = {ρ e } e∈E : R E m + × R E + → R E in
such a way that the observer provides a good estimate of the real density of the system.

A. Proposed solution

The procedure consists of two steps, namely offline calibration and online update. This choice is driven by the empirical observation that the standard triangular Diagram tends to overestimate the flow in congestion. Other solutions include the inverted-λ Fundamental Diagram [START_REF] Hall | Empirical analysis of freeway flow-density relationships[END_REF], but the resulting model involves hysteresis and becomes more complex.

Let e ∈ E m be a cell equipped with sensors, and denote by {(ρ k , ϕ k )} k∈K , K = 1, . . . , K, the set of K density and flow measurements used as learning set. Calibration of the Fundamental Diagram involves two steps (we write ρ c and C instead of ρ c e and C e for sake of notation) • Gradient descent algorithm for estimation of ρ c and C = v ff ρ c : in the first step, we estimate the critical density ρ c and the capacity C, namely, the nominal maximum flow, by solving the non-linear and nonconvex minimization problem

min (ρ c ,C) V (ρ c ,C) = k (ϕ k -ϕ (ρ c ,C) (ρ k )) 2 s.t. 0 < ρ c < ρ jam C > 0 ϕ (ρ c ,C) (x) = C ρ c x, x ≤ ρ c C(ρ jam -x) ρ jam -ρ c , x > ρ c (8)
To solve [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF] we propose the following gradient descent with diminishing stepsize algorithm -Basic step: initialize ρ c 0 , C 0 . A choice is ρ c 0 = 20 veh/km (vehicles start influencing each other when their relative distance is less than 50 meters), and C 0 = v limit e ρ c 0 , where v limit e is the speed limit on cell e normalized by the sampling time T ; -n-th step: let (ρ c n , C n ) descend along the gradient of the cost, namely

ρ c n+1 = ρ c n + δ n ∇ ρ c V (ρ c ,C) C n+1 = C n + δ n ∇ C V (ρ c ,C)
where the gradients

∇ ρ c V (ρ c ,C) and ∇ C V (ρ c ,C) are computed at (ρ c , C) = (ρ c n-1 , C n-1
), and δ > 0 is a fixed initial step size;

-Stopping criterion: stop if || ρ c n C n - ρ c n-1 C n-1 || < ε
for some small threshold ε > 0.

• Calibration of the congested part: given (ρ c , C), calibration of the quadratic function for the congested region can be cast as the following quadratic problem

min (a,b,c) k∈I FF (ρ c ) (ϕ k -(aρ 2 k + bρ k + c)) 2 s.t. a (ρ c ) 2 + bρ c + c e = C a ρ jam 2 + bρ jam + c = 0 a ≥ 0 (9)
which can be readily solved using off-the-shelf tools. The Fundamental Diagrams can be extended on cells e ∈ E \ E m by convex interpolation of the parameters ρ c , C, a, b and c, using those calibrated in the closest cells e ∈ E m , with weights proportional to the relative distance.

2) Online density reconstruction algorithm: The online algorithm proceeds as follows

• at the beginning of the t-th time slot receive the measurements {ϕ m e (t)} e∈E m . The vector of estimate of the outflows f out (t) is computed as solution of the following minimization problem

min f out ||(I -R T ) f out || 2 +γ e∈E m ( f out e -ϕ m e (t)) 2 s.t.
f out ≥ 0 (10) which aims both to match outflows and measurements and to balance the outflows according to the splitting ratios, as if the network were at steady state. The parameter γ trades off between the two: the higher, the more flows on cells in E m will be forced to match the measurements, the lower, the more the solution will be as if the system were at steady state. Finally, f in (t) = R T f out (t).

receive the measurements {v m e (t)} e∈E when available, or hold the previous measurements; -For each cell e, assume that the outflow can be ; -Set ρe (t) = arg min i=1,2 {|v e (ρ i e ) -v m e (t)|} as a pseudo-measurement of the density, Since it is solely based on latest measurements of flows and speed, it is very noisy, especially in congestion regime. Therefore, we don't use it directly; -For each cell e ∈ E, let the density estimate evolve according to the observer equation [START_REF] Lighthill | On kinematic waves. ii. a theory of traffic flow on long crowded roads[END_REF].

IV. DENSITY RECONSTRUCTION -EXPERIMENTAL

RESULTS

Our experimental setting is the Grenoble Traffic Lab (GTL) [START_REF] Canudas De Wit | Grenoble traffic lab: An experimental platform for advanced traffic monitoring and forecasting[END_REF], a network of sensors deployed for monitoring and research purposes along the "Rocade Sud", a 12 km long freeway encircling the town of Grenoble, France (Fig. 2, upper panel). The network consists of 135 magnetometers buried in the asphalt on both lanes of the main line and on each onramp and offramp, totalling 68 sensing locations. Fig. 2, upper panel, shows the position of the 22 sections of the main line equipped with sensors. With few exceptions, each sensing location consists of two magnetometers deployed in pairs at a known distance. On each sensing location and every T = 15 seconds, the system provides flow, occupancy and speed measurements. Occupancy o e is defined as the percentage of the last period of T = 15 seconds a vehicle was sitting over the sensor, and is approximatively related to the density as ρ e ≈ oe 100 ave , where ave is the average length of a vehicle (in km). Therefore, we assume that density is measured directly. Fig. 2, lower panel, shows a stylized representation of the Rocade Sud, the positions of the 68 sensing locations, and the distance between consecutive measurement sections along the main line. For a detailed report on the GTL, we refer to [START_REF] Canudas De Wit | Grenoble traffic lab: An experimental platform for advanced traffic monitoring and forecasting[END_REF].

We partition the network in cells in such a way that each cell comprises the space between two consecutive sensing locations. As such, the numbered circles in Fig. 2 also represent cells. We use Floating Car Data provided by INRIX, one of the main navigation and traffic monitoring companies. The Rocade has been partitioned into FCD segments, and average speed is computed in each FCD segment every 1 minute. FCD segments cover the whole main line of the Rocade and most onramps and offramps, but lanes are not distinguished along the main line. FCD segments are represented in Fig. 2 as rectangles encircling several sensing locations/cells.

For our experiments, and to prove that the method performs well even with sparse equipment, we consider a very limited subset of he available sensors. In particular, we only use the sensors on the 9 sections shown in light blue in Fig. 2, of which 8 correspond to positions in which loops inductors have been deployed by the Government Agency Centre national d'information routière (CNIR) [START_REF]Centre national d'information routière (CNIR)[END_REF] for monitoring purposes. No information from ramps is used.

We calibrate the Fundamental Diagrams employing data from April 10th, 2014, a working day (a Thursday) exhibiting standard traffic pattern:

• very limited night time traffic; • a peak of congestion in the morning (8:00 -10:00), triggered by vehicles exiting the Rocade at the offramp of Eybens (cells 37/38) and spilling back until Meylan; • a second congestion triggered by vehicles entering in A480 at Rondeau (unable to do so due to the high traffic on A480), and spilling back until around Libération; • a third, smaller, congestion triggered at Eybens around 14:00-15:00; • medium/heavy but fluid traffic from 10:00 to 16:00 • a second peak of congestion in the afternoon, again triggered by congestion at Rondeau at around 16:00, spilling back on the whole freeway in around 60 minutes, and lasting around two hours.

The matrix of splitting ratios is set according to the following rule of thumb:

• let e be a fast lane cell, and j and k be the following fast and slow lane cells. Then R ej = R ek = 0.5; • let e be a slow lane cell, and j and k be the following fast and slow lane cells. If among the cells that follow e there is not an offramp, then R ej = R ek = 0.5. Otherwise, R ej = R ek = 0.4 and R er = 0.2, where r is the offramp cell that follows e;

• if e is an onramp cell and j is the following slow ramp cell, then R ej = 1.

A. Validation

Our method was validated over a another Thursday working day, April 24th 2014. Calibration of the Fundamental Diagram at sensing locations is illustrated in Fig. 3, which shows in thick black the linear-quadratic Fundamental Diagram, in dashed thick black the standard linear Fundamental Digram in congestion regime, and as crosses the pairs (density, flow) measured on April 24th, 2014. As standard, data in freeflow are in good accordance with the linear part, while data in congested regime are much more scattered. As it can be observed, a standard triangular Fundamental Diagram would overestimate the flows in congested regime (dashed think line), while the convex quadratic curve better captures the average flow-density relation.

The algorithm was implemented in Matlab on a commercial laptop with 2.1 GHz i7-4600U CPU and 8 GB RAM. The optimization problems were solved using Matlab functions and the optimization system CVX [START_REF] Grant | Cvx: Matlab software for disciplined convex programming[END_REF], [START_REF] Grant | Graph implementations for nonsmooth convex programs[END_REF].

The results are reported in Fig. 4. For validation purposes only, density measurements from all GTL fixed sensors are considered ground truth. As such, the upper panel shows the evolution of the "true" measured density in all the cells on the main line, over the whole day. On the x-axis, the 46 sensing locations along the main line (numbers correspond to the labels in Fig. 2), on the y-axis, the 5760 time slots over the whole validation day (one slot every T = 15 seconds). Colors vary from green to yellow to red as density increases, with a minimum of 0 vehicles per km (green) to a maximum of ρ jam = 200 vehicles per km (red). In the lower panel, we show, in the same format, the results of the density reconstruction. The estimation algorithm captures the four congestions described above in a reasonably good way, given the limited amount of information employed; in particular, the two small congestions at Rondeau during the morning and at Eybens during early afternoon are both detected. The biggest discrepancy between estimates and measured data is around Eybens (cells 37-38) and can be explained as follows: The offramp of Eybens is a critical point at which many vehicles exit the Rocade, thus causing congestion before it. However, that point belongs to a long FCD segment running from Saint-Martin-d'Hères (cells 32-33) to Eybens entrance (cells 41-42), which thus provides a unique, rather low speed measurement, in turn yielding a high density estimate over the whole segment. A second obvious discrepancy is the smoothness of the reconstructed densities in congested regime compared with more scattered density measurements. While the former is related to the smoothing effect of the first-order observer, the latter is due to the high measurement rate of sensors, which well captures stop-and-go phenomena on the Rocade.

V. CONCLUSIONS

This formalizes the problem of data fusion of heterogeneous sources of information for density reconstruction, and proposes an easily implementable solution that employs sparse fixed sensor measurements and ubiquitous Floating Car Data. Calibration algorithms for the Fundamental Diagram are discussed. Future research directions include and are not limited to estimation of statistical properties of measurement noises, development of theoretical guarantees on maximum mean-square reconstruction error, calibration of the matrix of splitting ratios, and extension of the considered numerical scenario to part of the town of Grenoble.

1 ) 2 + b e ρ c e + c e = v ff e ρ c e a e ρ jam e 2 +

 122 Offline calibration: In this section we propose a simple solution for calibrating the Fundamental Diagram, namely, for estimating the function ϕ e (•) on the cells e ∈ E m . We adopt a linear-quadratic Fundamental Diagram ϕ e (ρ) = v ff e ρ, ρ ≤ ρ c e a e ρ 2 + b e ρ + c e , ρ > ρ c e where • ρ c e is the critic density, which separates the freeflow lowdensity region [0, ρ c e ), in which vehicles do not influence one each other, from the high-density congested region (ρ c e , ρ jam e ], in which speed decreases with density; • ρ jam e is the jam density, at which vehicles stop; • v ff e > 0 is the freeflow speed on cell e; • the flow in congested region is a convex quadratic function of the density, where for for consistency b e ρ jam e + c e = 0 a e > 0

-

  approximated by the local flow ϕ e = ϕ e (ρ e ) determined by the Fundamental Diagram. Compute the two possible densities ρ 1 e (freeflow) and ρ 2 e (congested) corresponding to flow f out e (t); For each cell e, compute the two corresponding velocities v e (ρ 1 e

Fig. 2 .

 2 Fig. 2. The experimental setting: the town of Grenoble and the Rocade Sud (upper panel), and a stylized version of the freeway (lower panel). The positions along the main line of the 22 sections of the main line in which sensors have been placed is shown as red pins. The positions of the 68 fixed sensing locations are shown in the stylized map. Each sensing location also corresponds to a cell. Light blue circles denote fixed sensors that are used in our implementation. Each rectangle represents one FCD segment, often providing average speed measurements over more than one cell.

Fig. 3 .

 3 Fig. 3. Calibration of the Fundamental Diagram on the cell Eybens exit -slow lane. The linear-convex Fundamental Diagram, calibrated using data from April 10th, 2014, is shown in thick line. Dashed thick line represents the corresponding linear Fundamental Diagram in congested regime. Each cross is a (flow, density) pair measured on April 24th, 2014, one for each time slot of T = 15 seconds during the whole day.

Fig. 4 .

 4 Fig. 4. Numerical results of density estimation on all cells: measured densities (upper panel) and estimated densities (lower panel).

The choice to set the speed measurement to freeflow speed for the first N samples, when no information is available, is done uniquely for the sake of simplicity, and can be replaced by any other educated guess.

VI. ACKNOWLEDGEMENTS

The authors wish to gratefully thank INRIX Europe, whose collaboration and permission to use data from the Grenoble area led to the results presented in this work.

This work as been partially

supported by the EU project SPEEDD FP7 n. 619435.