Cyril Furtlehner

Jean-Marc Lasgouttes

Arnaud De

La Fortelle

A Belief Propagation Approach to Traffic Prediction using Probe Vehicles

This paper deals with real-time prediction of traffic conditions in a setting where the only available information is floating car data (FCD) sent by probe vehicles. Starting from the Ising model of statistical physics, we use a discretized space-time traffic description, on which we define and study an inference method based on the Belief Propagation (BP) algorithm. The idea is to encode into a graph the a priori information derived from historical data (marginal probabilities of pairs of variables), and to use BP to estimate the actual state from the latest FCD. The behavior of the algorithm is illustrated by numerical studies on a simple simulated traffic network. The generalization to the superposition of many traffic patterns is discussed.

I. INTRODUCTION

With an estimated 1% GDP cost in the European Union (i.e. more than 100 billions euros), congestion is not only a time waste for drivers and an environmental challenge, but also an economic issue. Today, some urban and interurban areas have traffic management and advice systems that collect data from stationary sensors, analyze them, and post notices about road conditions ahead and recommended speed limits on display signs located at various points along specific routes. However, these systems are not available everywhere and they are virtually non-existent on rural areas.

In this context, the EU-funded REACT project developed new traffic prediction models to be used to inform the public and possibly to regulate the traffic, on all roads. The REACT project combines a traditional traffic prediction approach on equipped motorways with an innovative approach on nonequipped roads. The idea is to obtain floating car data from a fleet of probe vehicles and reconstruct the traffic conditions from this partial information. Two types of approaches are usually distinguished for traffic prediction, namely data driven (application of statistical models to a large amount of data, for example regression analysis) and model based (simulation or mathematical models explaining the traffic patterns). Models (see e.g. [START_REF] Chowdhury | Statistical physics of vehicular traffic and some related systems[END_REF], [START_REF] Klar | Mathematical models for vehicular traffic[END_REF] for a review) may range from microscopic description encoding drivers behaviors, with many parameters to be calibrated, to macroscopic ones, based on fluid dynamics, mainly adapted to highway traffic and subject to controversy [START_REF] Daganzo | Requiem for second order fluid approximation of traffic flow[END_REF], [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF]. Intermediate kinetic C. Furtlehner is with Project Team TAO at INRIA Futurs, Université de Paris-Sud 11 -Bâtiment 490, 91405 Orsay Cedex, France. (cyril.furtlehner@inria.fr)

J.-M. Lasgouttes is with Project Team IMARA at INRIA Paris-Rocquencourt, Domaine de Voluceau -BP 105, 78153 Rocquencourt Cedex, France. (jean-marc.lasgouttes@inria.fr) A. de La Fortelle is with INRIA Paris-Rocquencourt and École des Mines de Paris, CAOR Research Centre, 60 boulevard Saint-Michel, 75272 Paris Cedex 06, France. (arnaud.de la fortelle@ensmp.fr) description including cellular automata [START_REF] Nagel | A cellular automaton model for freeway traffic[END_REF] are instrumental for powerful simulation and prediction systems in equipped road networks [START_REF] Chrobok | Traffic forecast using simulations of large scale networks[END_REF]. On the other hand, the statistical approach mainly focuses on time series analysis on single road links, with various machine learning techniques [START_REF] Guozhen | Traffic flow prediction based on generalized neural network[END_REF], [START_REF] Chun-Hsin | Travel-time prediction with support vector regression[END_REF], while global prediction systems on a network combine data analysis and model simulations [START_REF] Chrobok | Traffic forecast using simulations of large scale networks[END_REF], [START_REF] Kanoh | Short-term traffic prediction using fuzzy c-means and cellular automata in a wide-area road network[END_REF]. For more information about traffic prediction methods, we refer also the reader to [START_REF] Benz | Information supply for intelligent routing servicesthe INVENT traffic network equalizer approach[END_REF], [START_REF] Versteegt | PredicTime -state of the art and functional architecture[END_REF].

We propose here a hybrid approach, by taking full advantage of the statistical nature of the information, in combination with a stochastic modeling of traffic patterns and a powerful message-passing inference algorithm. The beliefpropagation algorithm, originally designed for bayesian inference on tree-like graphs [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF], is widely used in a variety of inference problem (e.g. computer vision, coding theory. . . ) but to our knowledge has not yet been applied in the context of traffic prediction. The purpose of this paper is to give the first principles of such an approach, able to exploit both space and time correlation on a traffic network. The main focus is on finding a good way to encode some coarse information (typically whether traffic on a segment is fluid or congested), and to decode it in the form of real-time traffic reconstruction and prediction. In order to reconstruct the traffic and make predictions, we propose to use the socalled Bethe approximation of an underlying disordered Ising model (see e.g. [START_REF] Mézard | Spin Glass Theory and Beyond[END_REF]), to encode the statistical fluctuations and stochastic evolution of the traffic and the belief propagation (BP) algorithm, to decode the information. Those concepts are familiar to the computer science and statistical physics communities since it was shown [START_REF] Yedidia | Constructing free-energy approximations and generalized belief propagation algorithms[END_REF] that the output of BP is in general the Bethe approximation.

The paper is organized as follows: Section II describes the model and its relationship to the Ising model and the Bethe approximation. The inference problem and our strategy to tackle it using the belief propagation approach are stated in Section III. Section IV is devoted to a more practical description of the algorithm, and to numerical results illustrating the method. Finally, some new research directions are outlined in Section V.

II. TRAFFIC DESCRIPTION AND STATISTICAL PHYSICS

We consider a road network for which we want to both reconstruct the current traffic conditions and predict future evolutions. To this end, all we have is partial information in the form of floating car data arriving every minute or so. The difficulty is of course that the up-to-date data only covers part of the network and that the rest has to be inferred from this. In order to take into account both spatial and temporal relationships, the graph on which our model is defined is made of space-time vertices that encode both a location (road link) and a time (discretized on a few minutes scale). More precisely, the set of vertices is V = L ⊗ Z + , where L corresponds to the links of the network and Z + to the time discretization. To each point α = ( , t) ∈ V, we attach an information τ α ∈ {0, 1} indicating the state of the traffic (1 if congested, 0 otherwise). On such a model, the problems of prediction and reconstruction are equivalent, since they both amount to estimating the value of a subset of the nodes of the graph. The difference, which is obvious for the practitioner, lies mainly in the nature (space or time) of the correlations which are most exploited to perform the two tasks.

Each vertex is correlated to its neighbors (in time and space) and the evaluation of this local correlation determines the model. In other words, we assume that the joint probability distribution of τ

V def = {τ α , α ∈ V} ∈ {0, 1} V is of the form p({τ α , α ∈ V}) = α∈V φ α (τ α ) (α,β)∈E ψ αβ (τ α , τ β ) (1)
where E ⊂ V 2 is the set of edges, and the local correlations are encoded in the functions ψ and φ. V together with E describe the space-time graph G and V(α) ⊂ V denotes the set of neighbors of vertex α.

The model described by [START_REF] Chowdhury | Statistical physics of vehicular traffic and some related systems[END_REF] is actually equivalent to an Ising model on G, with arbitrary coupling between adjacent spins s α = 2τ α -1 ∈ {-1, 1}, the up or down orientation of each spin indicating the status of the corresponding link (Fig. 1). The homogeneous Ising model (uniform coupling constants) is a well-studied model of ferro (positive coupling) or anti-ferro (negative coupling) material in statistical physics. It displays a phase transition phenomenon with respect to the value of the coupling. At weak coupling, only one disordered state occurs, where spins are randomly distributed around a mean-zero value. Conversely, when the coupling is strong, there are two equally probable states that correspond to the onset of a macroscopic magnetization either in the up or down direction: each spin has a larger probability to be oriented in the privileged direction than in the opposite one. From the point of view of a traffic network, this means that such a model is able to describe three possible traffic regimes: fluid (most of the spins up), congested (most of the spins down) and dense (roughly half of the links are congested). For real situations, we expect other types of congestion patterns, and we seek to associate them to the possible states of an inhomogeneous Ising model with possibly negative coupling parameters, referred to as spin glasses in statistical physics [START_REF] Mézard | Spin Glass Theory and Beyond[END_REF].

The practical information of interest which one wishes to extract from ( 1) is in the form of local marginal distributions p α (τ α ), once a certain number of variables have been fixed by probe vehicles observations. They give the probability for a given node to be saturated at a given time and in turn can be the basis of a travel time estimation. From a computational viewpoint, the extraction cost of such an information from an Ising model on a multiply connected graph is known to scale exponentially with the size of the graph, so one has to resort to some approximate procedure. As we explain now such an approximation exists for dilute graphs (graphs with a tree-like local structure).

On a simply connected graph, the knowledge of p α (τ α ) the one-vertex and p αβ (τ α , τ β ) the two-vertices marginal probabilities is sufficient [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF] to describe the measure [START_REF] Chowdhury | Statistical physics of vehicular traffic and some related systems[END_REF].

p(τ V ) = (α,β)∈E p αβ (τ α , τ β ) α∈V p α (τ α ) qα-1 = α∈V p α (τ α ) (α,β)∈E p αβ (τ α , τ β ) p α (τ α )p β (τ β ) , (2)
where q α denotes the number of neighbors of α. Since our space time graph G is multi-connected, this relationship between local marginals and the full joint probability measure can only be an approximation, which in the context of statistical physics is referred to as the Bethe approximation. This approximation is provided by the minimum of the socalled Bethe free energy, which, based on the form (2), is an approximate form of the Kullback-Leibler distance,

D(b p) def = τV b(τ V ) ln b(τ V ) p(τ V ) ,
between the reference measure p and an approximate one b. This rewrites in terms of a free energy as

D(b p) = F(b) -F(p),
where

F(b) def = U(b) -S(b), (3)
with the respective definitions of the energy U and of the entropy S U(b

) def = - (α,β)∈E τα,τ β b αβ (τ α , τ β ) log ψ αβ (τ α , τ β ) - α∈V τα b α (τ α ) log φ α (τ α ), S(b) def = - (α,β)∈E τα,τ β b αβ (τ α , τ β ) log b αβ (τ α , τ β ) + α∈V τα (q α -1)b α (τ α ) log b α (τ α ).
The set of single vertex b α and and two-vertices b αβ marginal probabilities that minimize (3) form the Bethe approximation of the Ising model. For reasons that will become evident in Section III-B, these will also be called beliefs. It is known that the quality of the approximation may deteriorate in the presence of short loops. In our case, the fact that the nodes are replicated along the time axis alleviates this problem.

In practice, what we retain from an inhomogeneous Ising description is the possibility to encode a certain number of traffic patterns in a statistical physics model. This property is also shared by its Bethe approximation (BA) and it is actually easier to encode the traffic patterns in this simplified model rather than the original one. Indeed, it will be shown in Section III-C that the computation of the BA from the marginal probabilities is immediate.

The data collected from the probe vehicles is used in two different ways. The most evident one is that the data of the current day directly influences the prediction. In parallel, this data is collected over long periods (weeks or months) in order to estimate the model [START_REF] Chowdhury | Statistical physics of vehicular traffic and some related systems[END_REF]. Typical historical data that is accumulated is

• pα (τ α ): the probability that vertex α is congested (τ α = 1)
or not (τ α = 0); • pαβ (τ α , τ β ): the probability that a probe vehicle going from α to β ∈ V(α) finds α with state τ α and β with state τ β . The computation of pα and pαβ requires a proper congestion state indicator τ α that we assume to be the result of the practitioner's pretreatment of the FCD. The definition of this indicator is a problem of its own and is outside of the scope of this article. A relevant FCD variable is instantaneous speed. An empirical threshold may be attached to each link in order to discriminate (in a binary or in a continuous manner) between a fluid and a congested situation. Another approach is to convert the instantaneous speed in a probability distribution of the local car density, when an empirical fundamental diagram is known for a given link. Aggregation over a long period of these estimators yields then the desired historical data. The edges (α, β) of the space time graph G are constructed based on the presence of a measured mutual information between α and β, which is the case when pαβ (τ α , τ β ) = pα (τ α )p β (τ β ).

III. THE RECONSTRUCTION AND PREDICTION ALGORITHM

A. Statement of the inference problem

We turn now to our present work concerning an inference problem, which we set in general terms as follows: a set of observables τ V = {τ α , α ∈ V}, which are stochastic variables are attached to the set V of vertices of a graph. For each edge (α, β) ∈ E of the graph, an accumulation of repetitive observations allows to build the empirical marginal probabilities {p αβ }. The question is then: given the values of a subset τ V * = {τ α , α ∈ V * }, what prediction can be made concerning V * , the complementary set of V * in V?

There are two main issues:

• how to encode the historical observations (inverse problem) in an Ising model, such that its marginal probabilities on the edges coincide with the pαβ ? • how to decode in the most efficient manner-typically in real time-this information, in terms of conditional probabilities P (τ α |τ V * )? The answer to the second question will somehow give a hint to the first one.

B. The belief propagation algorithm

BP is a message passing procedure [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Network of Plausible Inference[END_REF], which output is a set of estimated marginal probabilities (the beliefs b αβ and b α ) for the measure (1). The name "belief" reflects the artificial intelligence roots of the algorithm. The idea is to factor the marginal probability at a given site in a product of contributions coming from neighboring sites, which are the messages. The messages sent by a vertex α to β ∈ V(α) depends on the messages it received previously from other vertices:

m α→β (τ β ) ← τα∈{0,1} n α→β (τ α )φ α (τ α )ψ αβ (τ α , τ β ), (4)
where

n α→β (τ α ) def = γ∈V(α)\{β} m γ→α (τ α ). (5)
The messages are iteratvely propagated into the network with a parallel, sequential or random policy. If they converge to a fixed point, the beliefs b α are then reconstructed according to

b α (τ α ) ∝ φ α (τ α ) β∈V(α) m β→α (τ α ), (6)
and, similarly, the belief b αβ of the joint probability of (τ α , τ β ) is given by

b αβ (τ α , τ β ) ∝ n α→β (τ α )n β→α (τ β ) × φ α (τ α )φ β (τ β )ψ αβ (τ α , τ β ). (7)
In the formulas above and in the remainder of this paper, the symbol ∝ indicates that one must normalize the beliefs so that they sum to 1. A simple computation shows that equations ( 6) and ( 7) are compatible, since (4)-( 5) imply that

τα∈{0,1} b αβ (τ α , τ β ) = b β (τ β ).
In most studies, it is assumed that the messages are normalized so that

τ β ∈{0,1} m α→β (τ β ) = 1.
holds. The update rule (4) indeed indicates that there is an important risk to see the messages converge to 0 or diverge to infinity. It is however not immediate to check that the normalized version of the algorithm has the same fixed points as the original one (and therefore the Bethe approximation). This point has been analyzed in [START_REF] Furtlehner | Belief propagation and Bethe approximation for traffic prediction[END_REF] and the conclusion is that the fixed points of both version of the algorithm coincide, except possibly when the graph has a unique cycle. We can safely expect that it is not the case in practical situations.

It has been realized a few years ago [START_REF] Yedidia | Generalized belief propagation[END_REF] that the fixed points of the BP algorithm coincide with stable points of the Bethe free energy (3), and that moreover stable fixed points correspond to local minima of (3) [START_REF] Heskes | Stable fixed points of loopy belief propagation are minima of the Bethe free energy[END_REF]. BP is therefore a simple and efficient way to compute the Bethe approximation of our inhomogeneous Ising model. We propose to use the BP algorithm for two purposes: estimation of the model parameters (the functions ψ αβ and φ α ) from historical data and reconstruction of traffic from current data.

C. Setting the model with belief propagation

The fixed points of the BP algorithm (and therefore the Bethe approximation) allow to approximate the joint marginal probability p αβ when the functions ψ αβ and φ α are known. Conversely, it can provide good candidates for ψ αβ and φ α from the historical values pαβ and pα .

To set up our model, we are looking for a fixed point of the BP algorithm satisfying ( 4)-( 5) and such that b αβ (τ α , τ β ) = pαβ (τ α , τ β ) and therefore b α (τ α ) = pα (τ α ).

It is easy to check that the following choice of φ and ψ,

ψ αβ (τ α , τ β ) = pαβ (τ α , τ β ) pα (τ α )p β (τ β ) , (8)
φ α (τ α ) = pα (τ α ), (9)
leads (1) to coincide with [START_REF] Klar | Mathematical models for vehicular traffic[END_REF]. They correspond to a normalized BP fixed point for which all messages are equal to 1/2. It has been shown in [START_REF] Furtlehner | Belief propagation and Bethe approximation for traffic prediction[END_REF] that this form of φ and ψ is in some sense canonical: any other set of functions yielding the same beliefs is equal to ( 8)-( 9), up to a change of variable. This equivalence holds for the beliefs at the other fixed points of the algorithm and their stability properties. This crucial result means that it is not needed to learn the parameters of the Ising model, but that they can be readily recovered from the Bethe approximation. The message update scheme (4) of previous section can therefore be recast as

m α→β (τ β ) ← τα∈{0,1} n α→β (τ α )p αβ (τ α |τ β ) (10)
and the beliefs are now expressed as

b α (τ α ) ∝ pα (τ α ) γ∈V(α) m γ→α (τ α ), (11)
b αβ (τ α , τ β ) ∝ pαβ (τ α , τ β )n α→β (τ α )n β→α (τ β ). ( 12
)
There is no guarantee that the trivial constant fixed point is stable. However, the following theorem, proved in [START_REF] Furtlehner | Belief propagation and Bethe approximation for traffic prediction[END_REF], shows that this can be decided from the mere knowledge of the marginal which we want to model.

Theorem 1: The fixed point {p} is stable if, and only if, the matrix defined, for any pair of oriented edges (α, β) ∈ E, (α , β ) ∈ E, by the elements J α β αβ = pαβ (1|1) -pαβ (1|0) 1 1 {α ∈V(α)\{β}, β =α} , has a spectral radius (largest eigenvalue in norm) smaller than 1.

A sufficient condition for this stability is therefore

pαβ (1|1) -pαβ (1|0) < 1 q α -1
, for all α ∈ V, β ∈ V(α).

In addition, on a dilute graph, the knowledge of the Jacobian coefficient distribution and the connectivity distribution of the graph is enough to determine the stability property by a mean field argument [START_REF] Furtlehner | Belief propagation and Bethe approximation for traffic prediction[END_REF].

D. Traffic reconstruction and prediction

Let V * be the set of vertices that have been visited by probe vehicles. Reconstructing traffic from the data gathered by those vehicles is equivalent to evaluating the conditional probability

p α (τ α |τ V * ) = p α,V * (τ α , τ V * ) p V * (τ V * ) , (13)
where τ V * is a shorthand notation for the set {τ α } α∈V * . The BP algorithm applies to this case if a specific rule is defined for vertices α ∈ V * : since the value of τ α is known, there is no need to sum over possible values and ( 10) becomes

m α→β (τ β ) ← n α→β (τ α )p αβ (τ α |τ β ). (14)

IV. PRACTICAL CONSIDERATIONS AND SIMULATION

The algorithm outlined in Section III can be summarized by the flowchart of Fig. 2. It is supposed to be run in real time, over a graph which corresponds to a time window (typically a few hours) centered around present time, with probe vehicle data added as it is available. In this perspective, the reconstruction and prediction operations are done simultaneously on an equal footing, the distinction being simply the time-stamp (past for reconstruction or future for prediction) of a given computed belief. The output of the previous run can be used as initial messages for a new run, in order to speedup convergence. Full re-initialization (typically a random set of initial messages) has to be performed within a time interval of the order but smaller than the time-scale of typical traffic fluctuations. We have tested the algorithm on the artificial traffic network shown on the program's screenshot of Fig. 3. To this end, we used a simulated traffic system which has the advantage to yield exact empirical data correlations. For real data, problems may arise because of noise in the historical information used to build the model; this additional difficulty will be treated in a separate work. The simulator implements a queueing network, where each queue represents a link of the traffic network (a single-way lane) and has a finite capacity. To each link, we attach a variable ρ ∈ [0, 1], the car density, which is represented by a color code on the user interface snapshot.

As already stated in Section II, the physical traffic network is replicated, to form a space time graph, in which each reconstruction algorithm

∀α ∈ V, β ∈ V(α) m α→β ← 1 it ← 0 ∀α ∈ V, β ∈ V(α) m old α→β ← m α→β it ← it + 1 Choose randomly α 0 ∈ V Q Q Q Q Q Q Q Q Q Q n y α 0 ∈ V * ?
Compute m α 0 →β for all β ∈ V(α 0 ) using [START_REF] Benz | Information supply for intelligent routing servicesthe INVENT traffic network equalizer approach[END_REF] Compute m α 0 →β for all β ∈ V(α 0 ) using ( 14)

Q Q Q Q Q Q Q Q Q Q n y it = itmax or m old -m <
Compute beliefs using [START_REF] Chrobok | Traffic forecast using simulations of large scale networks[END_REF] and [START_REF] Guozhen | Traffic flow prediction based on generalized neural network[END_REF] end Fig. 2. The traffic reconstruction algorithm. The parameters, besides those that have already been defined elsewhere, are the total number of iterations itmax, the maximal error > 0 and a norm • on the set of messages, which choice is up to the implementer. Likewise, the update policy described here is "random", but parallel or sequential updates can also be used.

vertex α = ( , t) corresponds to a link at a given time t of the traffic graph. To any space-time vertex α, we associate a binary congestion variable τ α ∈ {0, 1}. The statistical physics description amounts to relating the probability of saturation P (τ α = 1) to the density ρ α . For the sake of simplicity, we consider a linear relation and build our historical p according to some time averaging procedure. In practice, the car density would not be available from the FCD and a preprocessing of information would be necessary.

In our oversimplified setting, the single-vertex beliefs yield directly an estimation of the car density. Nevertheless, more realistic data collection and modeling would be completely transparent w.r.t. the algorithm.

To estimate the quality of the traffic restoration we use the following estimator:

reconstruction rate def = 1 |V| α∈V 1 1 {|bα-ρα|<0.2} ,
which computes the fraction of space-time nodes α for which the belief b α does not differ by more than an arbitrary threshold of 0.2 from ρ α . A typical prediction time series is shown in Fig. 4. The overall traffic level, characterized by some real number comprised between 0 and 1, oscillates between dense and fluid conditions with a certain amount of noise superimposed. In this setting, we observe first that BP has three fixed points, among which the reference b = p (see Section III-C), which is is in fact unstable because it is a superposition of distinct measures. The two additional fixed points represent actually the dense and fluid traffic conditions. These additional states appear spontaneously and some fine tuning is required to control saturation effects [START_REF] Furtlehner | Belief propagation and Bethe approximation for traffic prediction[END_REF]. This is reflected in the sudden drops of the reconstruction rate when the algorithm jumps from one state to the other during transient phases. A selection criteria based on free energy measurements may be used to choose the most relevant fixed point [START_REF] Furtlehner | Belief propagation and Bethe approximation for traffic prediction[END_REF].

Concerning the dependence of the reconstruction rate with respect to the number of probe vehicles, Fig. 5 indicates that the knowledge on less than 10% of the links (10 vehicles for 122 links) is sufficient in this setting to identify most of the time correctly the traffic regime. However, when the number of probes is increased, the reconstruction rate given by our algorithm saturates around 80%. In addition to the fact that time correlations are not incorporated in our model, the main reason for this saturation is that in our practical implementation [START_REF] Furtlehner | Belief propagation and Bethe approximation for traffic prediction[END_REF] correlations between probe vehicles are neglected when imposing (13).

V. CONCLUSION AND PERSPECTIVES

We have presented a novel methodology for reconstruction and prediction of traffic using the belief propagation algorithm on floating car data. We have shown how the underlying Ising model can be determined in a straightforward manner and that it is unique up to some change of variables. The algorithm has been implemented and illustrated using an artificial traffic model. While our main focus is currently on testing on more realistic simulations, several generalizations are considered for future work, using the extension of the results of Section III to a more general factor graph setting done in [START_REF] Furtlehner | Belief propagation inference with a prescribed fixed point[END_REF]:

Firstly, the binary description corresponding to the underlying Ising model is arbitrary. Traffic patterns could be represented in terms of s different inference states. A Potts model with s-states variables would leave the belief propagation algorithm and its stability properties structurally unchanged. However, since the number of correlations to evaluate for each link is s 2 -1, this number of states should be subject to an optimization procedure.

Secondly, our way of encoding traffic network information might need to be augmented to cope with real world situations. This would simply amount to use a factor-graph used to propagate this information. In particular it is likely that a great deal of information is contained in the correlations of local congestion with aggregate traffic indexes, corresponding to sub-regions of the traffic network. Taking these correlations into account would result in the introduction of specific variables and function nodes associated to these aggregate traffic indexes. These aggregate variables would naturally lead to a hierarchical representation of the factor graph, which is necessary for inferring the traffic on large scale network. Additionally, time dependent correlations which are needed for the description of traffic, which by essence is an out of equilibrium phenomenon, could be conveniently encoded in these traffic index variables.

Ultimately, for the elaboration of a powerful prediction system, the structure of the information content of a trafficroad network has to be elucidated through a specific statistical analysis. The use of probe vehicles, based on modern communications devices, combined with a belief propagation approach, is in this respect a very promising approach.

Fig. 1 .

 1 Fig. 1. Traffic network (a) and Ising model (b) on a random graph. Up (resp. down) arrows correspond to fluid (resp. congested) traffic.

Fig. 3 .Fig. 4 .Fig. 5 .

 345 Fig. 3. Traffic network as produced by the simulator. The continuous color code represents the traffic index from 0 (green/light) to 1 (red/dark). There are 35 physical nodes and 122 physical links (road segments), simulated on 40 time steps, which yields a time-space graph G with 4880 nodes.


Eliane Giraud
email: eliane.giraud@ulg.ac.be

Michel Suery

Michel Coret

M Suéry

High temperature compression behavior of the solid phase resulting from drained compression of a semi-solid 6061 alloy

High temperature compression behavior of the solid phase resulting from drained compression of a semi-solid 6061 alloy

Introduction

The rheological behavior of alloys in the solidification range is now studied more and more extensively [1][2][3][4][5] owing to the importance of processes during which solid and liquid phases are coexisting. This coexistence occurs obviously during conventional solidification of castings, ingots or billets but also during liquid phase sintering, welding and forming such as rheocasting or thixocasting [6][7][8][9][10]. The knowledge of this behavior is indeed important for modeling purposes to avoid numerous trial and error experiments. For example, during semi-continuous casting of Al billets, defects like hot tears or macrosegregations [11][12][13] can form depending on alloy composition and process parameters. Their prediction requires modeling the whole process by considering the behavior of the semi-solid alloy taking into account both the deformation of the solid and the flow of the liquid. Criteria for hot tear formation are then introduced based on critical solid deformation or cavitation pressure in the liquid [14,15].

To determine the rheological behavior of the solid, it is generally assumed that its composition is not far from that of the alloy. Experiments are then carried out at temperatures close to but below the solidus temperature of the alloy to determine the various parameters of the constitutive equation. In this temperature range, viscoplastic behavior is generally a good approximation so that the constitutive equation is mainly determined by the strain rate sensitivity parameter and by the activation energy [16][17][18]. This procedure therefore does not take into account the fact that the composition of the solid phase is not that of the alloy and also that it is changing with temperature and thus with the solid volume fraction present in the alloy.

Indeed, during solidification or partial melting of a binary alloy, solid and liquid are coexisting with compositions given by the phase diagram and proportions by the lever rule in equilibrium conditions. The composition and the proportion of these two phases thus continuously change with temperature. To determine the constitutive equation of the solid phase, it is therefore necessary to test alloys with various compositions corresponding to those of the solid phase at various temperatures below the solidus. Extrapolation of the results at the solidus temperature for each composition allows determining the behavior of the solid phase in a semi-solid alloy at various temperatures. This procedure can be quite easily considered in the case of a binary alloy but it is hardly possible in a multi-constituent alloy. In this case indeed, the composition of the solid phase is not simply given by the phase diagram and even the solid can be constituted of several phases. Specific softwares are thus required which are able to give the composition and the proportion of the various phases as a function of temperature for various solidification or partial melting conditions. The next step of the procedure would be to prepare alloys with the composition of the solid phase at various temperatures and to test them.

Another simpler procedure can be considered, i.e. drainage of the liquid from the semi-solid alloy and testing of the remaining solid. If all the liquid present at a given temperature can be drained out of a sample, the remaining solid has the exact composition of the solid phase at this temperature. Repeating this procedure at various temperatures would then lead to various specimens having the composition of the solid phase at these temperatures. One way to carry out this drainage of the liquid is to use drained oedometric compression, already performed on aluminum alloys [2,19]. It consists in compressing a semi-solid sample placed in a container with a filter on top of it by applying the compression using a hollow piston to drain the liquid through the filter. This procedure allows one to measure also the compressibility of the solid skeleton provided that the pressure required to drain the liquid is low compared to that required to compress the solid.

In the present work, the feasibility of the procedure exposed above has been studied in the case of a 6061 alloy. Drained compressive tests have been carried out at various temperatures within the solidification range. These tests allow measuring the compressibility of the solid for various temperatures and obtaining specimens representative of the solid phase found in the semi-solid state. Then compressive tests at high temperature (below the solidus temperature) have been carried out on these drained specimens in order to determine the constitutive equation of the solid. Finally, a comparison with the behavior of the non-drained 6061 alloy within the same temperature range has been performed.

Experimental procedure

The alloy used for this investigation is a 6061 alloy provided by Almet (France) in the form of a rolled plate of 50 mm thickness in the heat treated T6 condition. In order to drain the liquid from a semi-solid sample, a drained compression apparatus has been designed. It consists in a container of 35 mm diameter in which the alloy is placed while being still solid (Fig. 1). During the test, the alloy is initially melted, and then partially solidified at a cooling rate of 20 K/min until a given solid fraction is reached. At this stage, the temperature is kept constant, and a downward vertical displacement is imposed to the hollow piston. Two displacement rates of the piston have been studied: 0.008 and 0.015 mm/s. A system of filtration allows the liquid to flow out of the sample, and, consequently, the solid fraction increases. In fact this system is constituted of two different filters as shown in Fig. 1: a rigid stainless steel filter with quite large holes (about 2 mm diameter) associated with a much thinner stainless steel filter with very small holes (about 200 m diameter). The first filter allows transmitting the load from the hollow piston whereas the second is required for the filtration of the liquid. This test can be seen as a means to impose solidification mechanically at constant temperature. Assuming that the solid fraction evolves only by liquid drainage, the following equation can be used to relate the imposed axial strain ε z to the solid fraction g s :

g s = g s0 exp(-ε z ) (1)
where g s0 is the initial solid fraction before any strain is imposed. Two values of this initial solid fraction (0.6 and 0.8) have been selected by using two temperatures (910 K and 897 K, respectively) at which oedometric compression has been carried out. These two temperatures have been determined by using the ProPhase software (from ALCAN CRV, France). The compression of the sample is pursued until g s is theoretically equal to 1, which means that all the liquid should have been drained from the semi-solid sample. The test was thus stopped when ε z satisfied this condition. After drained compression, the solid remaining in the container has been cut in two pieces through the axis of the cylinder to observe the microstructure of the alloy after drainage of the liquid.

These two pieces have been thereafter machined in order to get compression specimens of 8 mm in diameter and 6 mm in height. These specimens have been used for compression tests carried out at high temperature by using the strain rate jump procedure. Four temperatures have been investigated: 723, 773, 803 and 823 K as well as five strain rates: 10 -4 , 2.5 × 10 -4 , 10 -3 , 2.5 × 10 -3 and 6 × 10 -3 s -1 .

Similar compression tests have been also carried out for comparison on samples machined in the non-drained 6061 rolled plate.

Experimental results

Drained compression tests

Influence of initial solid fraction

Fig. 2 shows the variation of the applied stress as a function of the theoretical solid fraction (given by Eq. ( 1)) for the two tests carried out with different initial solid fractions at a displacement rate of the piston of 0.015 mm/s.

The applied stress obviously increases with increasing strain i.e. with increasing solid fraction. It increases slowly initially when the initial solid fraction is small (0.60) and much more rapidly when the initial solid fraction is larger (0.80). It should be noted that the two curves are coming close to each other at large solid fractions but the curve corresponding to the larger initial solid fraction is always below that for the lower initial solid fraction. In addition, it is surprising to observe that the stress does not increase very sharply when the solid fraction is close to 1 although the solid is not compressible.

Influence of displacement rate

Fig. 3 shows the variation of the applied stress as a function of the solid fraction for the two tests carried out with different displacement rates of the piston and for an initial solid fraction of 0.8.

The displacement rate has an influence on the stress required to drain the liquid from the specimen: a high strain rate leads to a higher measured stress at a given solid fraction. This strain rate sensitivity is correlated to the viscoplastic behavior of the solid network in this solid fraction range as already observed for tensile [20] and shear [21] experiments.

3.2.

Compressive tests at high temperature (below the solidus temperature)

Compressive tests at constant strain rate on drained samples

Since most of the solidification defects form in the last stage of solidification (i.e. g s > 0.8), only samples machined in the remaining solid obtained after a drained experiment at a high initial solid fraction (i.e. g s0 = 0.8) have been tested in order to determine the high temperature behavior of the solid phase.

Fig. 4 shows the stress-strain curves obtained at the various temperatures and for the different strain rates applied successively from 10 -4 s -1 to 6 × 10 -3 s -1 . For a given strain rate, stress increases until a plateau where it remains relatively constant, thus showing negligible strain hardening during deformation of the solid material. Moreover, stress obviously increases with increasing strain rate and decreasing temperature, which highlights the viscoplastic behavior of the solid material.

Compressive tests at constant strain rate on non-drained 6061 samples

The same type of compression experiments at various temperatures in the solid state as for drained samples has been performed on non-drained 6061 samples. Fig. 5 shows the stress-strain curves obtained at the various temperatures and for the different strain rates applied successively from 2.5 × 10 -4 s -1 to 6 × 10 -3 s -1 . Negligible strain hardening and viscoplastic behavior of the solid material are again found.

Discussion

The behavior of the semi-solid alloy during the drained experiments will be discussed first. Then the behavior observed during compression at high temperature of the solid phase will be examined in order to determine the appropriate rheological law. Finally, a comparison with the behavior at high temperature of the nondrained 6061 alloy will be performed.

Behavior of the semi-solid alloy during drained compressive tests

The stress required for the drained oedometric compression of the partially solidified specimen increases with increasing strain or solid fraction (Figs. 2 and3). This stress is due to both the densification of the solid skeleton and the liquid flow out of the sample which becomes more and more difficult as the remaining liquid volume fraction decreases. Indeed, a decrease of the liquid volume fraction allows increasing the number of contacts between dendrites and decreasing the interdendritic spaces thus decreasing the permeability of the solid skeleton. Adopting Darcy's law to describe liquid flow and assuming homogeneous deformation in the specimen, its contribution to the stress measured during the drained compression test can be evaluated. The maximal interstitial pressure P L , which corresponds to the pressure at the bottom of the container is then given by [22]:

P L = V p • Á 2 • K • H 0 • h 2 + P 0 ( 2
)
where V p is the velocity of the piston, Á is the liquid viscosity, K is the solid skeleton permeability, H 0 and h are the initial and current positions of the piston from the bottom of the container and P 0 is the pressure due to the drained liquid located above the filters. The variation of the position of the piston during the test is given by: h = H 0 • g S0 /1g L where g S0 is the initial solid fraction and g L is the current liquid fraction. To determine the solid skeleton permeability, the Kozeny-Carman equation can be used [23,24]:

K = g 3 L • 3 1 /5
where 1 is the primary dendrite arm spacing. 1 is assumed equal to that given by Flemings for pure aluminum in [START_REF] Flemings | Solidification Processing[END_REF] since, to measure this parameter during our experiments, it would have been necessary to cool with a very high rate the specimen when a solid fraction of 0.8 is reached, which is not possible. P 0 is given by: P 0 = (H 0h) • L • g with L the liquid density and g the gravity acceleration.

The variation of interstitial pressure with liquid fraction, calculated using Eq. ( 2), is shown in Fig. 6: it increases slowly up to a solid fraction of 0.97 and then very sharply when the solid fraction exceeds 0.97 which corresponds to the coalescence solid fraction [20,21]. This sharp increase is due to the fact that, for solid fractions larger than that for coalescence (>0.97), the number of solid bridges increases drastically thus leading to a drop of the solid skeleton permeability and then to more and more difficult liquid flow. However, the interstitial pressure remains much lower (0.07 MPa) than the experimental measured stress (12 MPa as shown in Figs. 2 and3) so that it can be neglected during the deformation of the semisolid alloy. The mushy material therefore behaves during drained compression like a solid without liquid.

Metallographic observation of the part of the specimen which has been drained through the filter (Fig. 7) and of the part which remained in the container (Fig. 8) shows that liquid has been effectively drained from the specimen.

The solidified drained liquid (Fig. 7) contains much more eutectic and intermetallics than the non-drained 6061 alloy [20]. Since this liquid was drained from the specimen when the solid fraction was equal to 0.8, it is possible to determine its composition by using the ProPhase software for the solidification conditions of the experiment (20 K/min) (Table 1). The liquid is theoretically enriched in Mg, Si, Fe and Cu which are the main secondary elements of the 6061 alloy. An analysis by X-ray diffraction of the drained liquid, by using the Cu K␣ radiation, allows detecting various phases which contain these elements such as: Mg 2 Si, CuAl 2 , Al 15 (FeMn) 3 Si 2 and Al 4 Cu 2 Mg 8 Si 7 in addition to the Al-rich matrix (Fig. 9). This result confirms that liquid has been drained from the specimen. The solid part which remains in the container after drainage contains some intermetallics homogeneously distributed everywhere in the specimen (Fig. 8). Their concentration is not high enough to allow detection by X-ray diffraction: as shown in Fig. 9, only the Al matrix is detected. Since these intermetallics are observed after cooling, it is necessary to wonder whether they were already present in the material before liquid drainage or they formed upon cooling after liquid drainage. Prophase calculation indicates that the first intermetallics form in the 6061 alloy when the volume fraction of the primary phase is 0.84. Therefore, the intermetallics present in the solid part are more likely to result from the solidification of the liquid which was not drained from the specimen. Thus, the drainage of the liquid has not been complete which seems to be in contradiction with the stress-solid fraction curves (Figs. 2 and3): these curves indicate that the drained compression tests were stopped when the solid fraction was theoretically very close to 1, so that almost no liquid should have remained in the specimen. The fact that some liquid remains can be explained by two factors: the filter has been deformed slightly at the end of the compression test so that the real axial strain applied to the specimen was smaller than expected, and/or some solid was extruded through the filter. These assumptions allow us to explain the behavior observed when the apparent solid fraction reaches 1. Indeed, since the solid is not compressible, the stress should have increased asymptotically to infinite, which is not observed experimentally.

In addition to the viscoplastic behavior of the solid phase shown in Fig. 3, the influence of the initial solid fraction and of the accumulated strain on the mechanical behavior of the semi-solid material has been investigated. Fig. 2 also shows that stress increases more or less rapidly with increasing strain depending on the initial solid fraction. This is due to the initial morphology of the solid phase: when the initial solid fraction is high, the solid grains are more connected so that drained compression involves extensive deformation of the solid with less rearrangement thus leading to an important increase of stress with increasing strain. Fig. 2 shows also that, at a given solid fraction, the measured stress is higher when the initial solid fraction is low. This result can be explained by the level of accumulated strain required to reach this solid fraction, which is larger when the initial solid fraction is low. A larger accumulated strain obviously leads to a more deformed microstructure and consequently to a larger number of solid-solid contacts. However, when the strain is such that the solid fraction approaches 1, the microstructure is largely deformed with almost no liquid remaining, but, as strain hardening is not observed in this temperature range, the stress no longer depends on the initial solid fraction.

Behavior of the drained solid at high temperature

Fig. 4 has shown that the solid phase exhibits a viscoplastic behavior with negligible strain hardening. The most common approximation used to describe the mechanical behavior of such material consists in using a classical creep law as follows [16,18]:

ε = A • n • exp - Q RT ( 3
)
where ε is the strain rate, is the measured stress, n is the stress sensitivity parameter, Q is the activation energy, A is a material constant, T is the temperature and R is the ideal gas constant. The behavior of the solid material is thus described by two main parameters: n and Q. The slope of the curves showing the variation of stress as a function of strain rate in logarithmic scales (Fig. 10) allows determining the strain rate sensitivity parameter m which is close to 0.08 whatever the temperature. The stress sensitivity parameter n is then equal to 12 (n = 1/m). From the curves showing the variation of ln( ε/ n ) as a function of -1/RT (Fig. 11), it is possible to deduce the activation energy Q which is equal to 347 kJ/mol.

The values obtained for n and Q are very high. Indeed, as a general rule, the stress sensitivity parameter for materials deformed at high temperature is between 3 and 5. An exponent of 3 is generally attributed to dislocation glide controlled by viscous drag whereas a value of 5 corresponds to climb controlled dislocation glide [START_REF] Kloc | [END_REF]. A value of 8 can be sometimes observed when the substructure does not evolve with stress. The activation energy is usually around 130-140 kJ/mol which is the activation energy for self-diffusion in 4) for specimens resulting from drained compression of the 6061 alloy at an initial solid fraction of 0.8.

Table 2

Parameters for the rheological law (4) and for specimens resulting from drained compression of the 6061 alloy at an initial solid fraction of 0.8. aluminum [27,28]. Thus, the classical creep law (3) is not adapted to describe the behavior of the solid phase of the 6061 alloy since the values of the various parameters have no physical meaning. Another possibility consists in using a hyperbolic sine law. Indeed, some alloys can exhibit a stronger stress increase than that predicted by an exponential law when exceeding a given stress or strain rate level [18,[29][30][31][32]. This leads to the presence of two regimes: a power-law region and a power-law breakdown region, presumably governed by different physical processes not already completely known. For intermediate strain rates and stresses, the alloy exhibits a behavior governed by a power law, whereas, for high strain rates and stresses, a hyperbolic sine law intervenes. Therefore, due to the stress level reached during the experiments, it can be assumed that the experimental conditions lead to the typical behavior of the power-law breakdown regime. In this case, the following equation can be used to relate the strain rate ε to the experimental stress [18,29,30]:

ε = A • G • b k • T • exp -Q RT • [sinh(˛• )] n (4)
where b is the burgers vector (2.86 × 10 -10 m), k is the Boltzman constant (1.38 × 10 -23 J/K), A and ˛are material constants and G is the shear modulus. It is assumed that the shear modulus depends on temperature according to the equation suggested in [33]: G = 3.022 × 10 4 -16 × T (in MPa).

The standard procedure to determine the parameters values can be divided into two main steps: (i) determination of the activation energy Q and the stress sensitivity parameter n by applying Eq. (3) on experimental data within the power-law regime and (ii) determination of the material constants A and ˛by applying Eq. (4) on experimental data within the power-law breakdown regime and by using the values of Q and n calculated in (i). This procedure assumes that sufficient data for each regime are available, which is not the case in this study. Therefore, the values of Q and n have been fixed to classical values for hot deformation of aluminum alloys [34,35] i.e. 142 kJ/mol and 5, respectively. Parameter's and material's constants A and ˛have then been determined so that the variation of ε • k • T • exp(Q/RT )/G • b as a function of sinh(˛• ) in a logarithmic scale exhibits a slope equal to 5, as shown in Fig. 12. The values of the different parameters of Eq. ( 4) are summarized in Table 2. It can be noticed that the material constant ˛is in agreement with classical values of the literature for aluminum alloys: from 0.01 to 0.08 [29][30][31]36], which tends to show that a hyperbolic sine law is relatively appropriate to describe the behavior of the solid phase of the 6061 alloy. The same study as for the drained alloy has been performed on the 6061 alloy in order to determine the most appropriate rheological law describing its behavior at high temperature.

Since Eq. ( 3) still leads to high values for the activation energy Q (about 287 kJ/mol) and the stress sensitivity parameter (about 11), an identical procedure has been followed to estimate the values of the parameters for Eq. ( 4). These values are summarized in Table 3.

Although the behavior of the non-drained 6061 alloy is also governed by a hyperbolic sine law, it can be noticed, by comparing Tables 2 and3, that the values of the constitutive parameters are different depending on the type of tested material. Therefore, the deformation at high temperature of the non-drained 6061 alloy differs from the deformation of the drained alloy, which confirms the necessity to study the behavior of the material with the appropriate composition to determine the response of the solid phase in the semi-solid state.

Conclusion

In this study, an original technique to determine the behavior at high temperature of the solid phase of a multi-constituent alloy was presented. It consists in: (i) drainage of the liquid present at a given temperature in order to obtain a solid with the exact composition of the solid phase at this temperature and (ii) deformation of the solid at high temperature to determine the rheological law.

The densification behavior of the 6061 alloy in the mushy state has been investigated thanks to drained compressive tests. These tests consist in a mechanical solidification of the semi-solid alloy at a constant temperature resulting from liquid drainage through filters which leads to a nearly complete densification of the solid phase. The results have shown that the compression behavior of the semi-solid alloy depends on strain rate, on the initial morphology of the solid skeleton and on the accumulated strain. Indeed, at a given solid fraction achieved during compression, the pressure required to densify the solid increases with decreasing initial solid fraction and thus increasing accumulated strain. Since strain hardening is not present in this temperature range, this result is correlated to the liquid distribution in the specimen which depends on the accumulated strain.

Specimens resulting from drained compression have been further tested in simple compression at high temperature to determine the constitutive equation of the solid phase present in a semi-solid 6061 alloy at a given temperature and to compare it with that of the 6061 alloy. Compression tests with strain rate jumps at various temperatures have shown that the behavior of the two materials is governed by a hyperbolic sine law. The values of the parameters of the constitutive equation are however different. This thus shows that the behavior of the solid phase in the 6061 alloy differs from that of the alloy. It is therefore necessary to study the behavior of the material with the appropriate composition if the behavior of the solid phase must be known within the solidification range.

Fig. 1 .

 1 Fig. 1. Schematic view of the oedometric compression apparatus (a) with picture of the two filters (b).

Fig. 2 .

 2 Fig. 2. Variation of the applied stress during drained compression as a function of the solid fraction present in the specimen. The two curves correspond to two different initial solid fractions.

Fig. 3 .

 3 Fig. 3. Variation of the applied stress during drained compression as a function of the solid fraction present in the specimen. The two curves correspond to two different displacement rates of the piston.

Fig. 4 .

 4 Fig. 4. Stress-strain curves at various temperatures and various strain rates applied by compression on specimens resulting from drained compression of the 6061 alloy at an initial solid fraction of 0.8.

Fig. 5 .

 5 Fig. 5. Stress-strain curves at various temperatures and various strain rates applied by compression on non-drained 6061 specimens.

Fig. 6 .

 6 Fig.6. Pressure necessary to drain the liquid through a porous solid mediumcalculated using Eq. (2).

Fig. 7 .

 7 Fig. 7. Microstructure of the liquid part of the 6061 specimen which has been drained through the filters for an initial solid fraction equal to 0.8: overall view (a) and magnified view (b).

Fig. 8 .

 8 Fig. 8. Microstructure of the solid part of the 6061 specimen remaining in the container after a drained compression at an initial solid fraction of 0.8: in the vicinity of the filters (a) and in the bottom of the specimen (b).

Fig. 9 .

 9 Fig. 9. X-ray diffraction patterns on the liquid part of the 6061 specimen which has been drained through the filters (a) and on the solid part which remained in the container (b) for an initial solid fraction equal to 0.8.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Stress-strain curves in a logarithmic scale at various temperatures for specimens resulting from drained compression of the 6061 alloy at an initial solid fraction of 0.8.

Fig. 12 .

 12 Fig. 12. Determination of the material constants A and ˛from Eq. (4) for specimens resulting from drained compression of the 6061 alloy at an initial solid fraction of 0.8.

4. 3 .

 3 Comparison with the behavior at high temperature of the non-drained 6061 alloy

Table 1

 1 Composition (in wt%) of the liquid in a 6061 alloy when the solid fraction is equal to 0.8 obtained by ProPhase calculation.

	Si	Mg	Cu	Fe	Mn	Cr
	2.5	2.9	1.1	1.2	0.2	0.06

Table 3

 3 Parameters for the rheological law (4) and for non-drained 6061 specimens.

	n	Q (kJ/mol)	˛(MPa -1 )	A (m 2 s -1 )
	5	142	0.055	6 × 10 -15

Acknowledgements

One of the authors (EG) is grateful to CNRS (French National Center for Scientific Research) and AREVA for financial support through a scholarship. The authors thank Cédric Gasquères, ALCAN CRV (France), for providing the ProPhase calculations and Stéphane Coindeau, CMTC (France), for the analysis by X-ray diffraction.