A DISTRIBUTED MODEL OF TRAFFIC FLOWS ON EXTENDED REGIONS

. This work deals with the modelling of traﬃc ﬂows in complex networks, spanning two-dimensional regions whose size ( macroscale ) is much greater than the characteristic size of the network arcs ( microscale ). A typical example is the modelling of traﬃc ﬂow in large urbanized areas with diameter of hundreds of kilometers, where standard models of traﬃc ﬂows on networks resolving all the streets are computationally too expensive. Starting from a stochastic lattice gas model with simple constitutive laws, we derive a distributed two-dimensional model of traﬃc ﬂow, in the form of a non-linear diﬀusion-advection equation for the particle density. The equation is formally equivalent to a (non-linear) Darcy’s ﬁltration law. In particular, it contains two parameters that can be seen as the porosity and the permeability tensor of the network. We provide suitable algorithms to extract these parameters starting from the geometry of the network and a given microscale model of traﬃc ﬂow (for instance based on cellular automata). Finally, we compare the fully microscopic simulation with the ﬁnite element solution of our upscaled model in realistic cases, showing that our model is able to capture the large-scale feature of the ﬂow.


Introduction.
Modeling traffic dynamics is the subject of a large literature. Since the early seminal works in the 1990s (we refer the reader to the review [10] for a detailed chronological presentation), different approaches were developed, focusing on different scales. All models are based on a graph or network of streets, in which vehicles can move. In particular, macroscopic models were developed in which the main variables are average quantities such as the mean vehicle density 526 FABIO DELLA ROSSA, CARLO D'ANGELO AND ALFIO QUARTERONI and the mean vehicle velocity at each point of each street of the network. Macroscopic models are expressed by partial differential equations (PDE) representing conservation laws. The typical length scale for such models is usually by far greater than the typical size of a single vehicle. On the other hand, microscopic models were developed in which the interactions between vehicles are modelled individually, based on either cellular automata (CA) or analogous discrete evolutionary models. In spite of the microscopic models being more accurate, their main drawback is the huge computational cost when applied to complex networks.
Unfortunately, when traffic flows have to be evaluated in very large regions with complex road networks, the computational cost of traffic simulations can be prohibitive even for macroscopic models. The problem is related to the length scale associated to single roads, which is very small compared to the size of the considered region. Hence, a natural question is whether the details of the geometry of the single roads could be taken into account by averaged quantities distributed over the region, and whether a suitable evolutionary model for the relevant variables such as the average vehicle density and the mean vehicle velocity could be developed using such information. Notice that such model would be a macroscopic model of traffic flow over a 2D region, in which we only keep trace of the underlying road network by means of some distributed quantities related to the original "microscopic" geometry.
Similar problems in geophysics have motivated the use of "homogenized" models, based on upscaling techniques, to study fluid flow in heterogeneous media. In this work we want to focus on upscaling techniques for traffic flows in complex networks. In particular, we develop and study homogenized models of traffic flow and compare the results with those obtained by more standard approaches, such as CA techniques, on complex networks. Our model is derived on the basis of simple assumptions and standard asymptotic expansion techniques. Its formulation is similar to classic models of fluid flow in porous media, with a non-linear advection term. Moreover, in the simple case of one road, it reduces to well-known macroscopic LWR model of traffic flow, introduced indipendently by Lighthill, Whitham [13] and Richards [19] in the 50's.
There are several motivations that make such homogenized approach attractive. First, they allow for a fast computational analysis of traffic flows. Second, comparing the macroscopic model with more accurate microscopic models, we will show that the relevant feature of the flow dynamics are captured by the homogenized version. This is possible since some information about the microscopic network geometry is contained in the homogenized quantities, entering the model as parameters. In particular, such quantities have the (suggestive) meaning of "porosity" and "permeability" of the road network thought as a continuum. In several cases, those quantities are interesting per se, as they carry the essential local information about connectivity and anisotropy [14] of the network. Hence, our analysis may be helpful to better understand the relevant feature of an urban tissue. Indeed, the methodologies presented in this work have been used in collaboration with Studio08 (prof. B. Secchi and P. Viganò, IUAV Unversity, Venice) to analyze the urban morphology of the greater Paris area on a length scale of 100 Km, in the framework of the international research consulting programme Le Grand Pari(s) de l'agglomération parisienne [11].
We will derive our distributed model starting from a discrete one. We consider a population of particles (vehicles) evolving on a two-dimensional regular lattice (lattice gas model ). The lattice is assumed uniformly spaced in the x and y direction in the plane. The coordinates of the lattice nodes will be denoted by (x m , y n ) = (mh, nh), being m, n ∈ Z, where h > 0 is the characteristic lattice size (see fig. 1). Each node is representing a site that can be occupied by a particle. The lattice gas evolves according to simple rules: a) Each site can be occupied (by one particle) or free. b) Each particle can move from a certain (occupied) site only to the free neighboring ones, with some prescribed probabilities. c) Particle moves can happen only at discrete times t k = k∆t, k ∈ N, ∆t > 0 being a characteristic time.
In the one-dimensional case, Brieger and Bonomi [1] highlighted the equivalence between the classic LWR model and this lattice gas model. They took a collection of sites x m = mh on a line, assuming that if a particle can move, it moves forward with probability α + ∈ [0, 1], and backward with probability α − = 0 (no reverse flow). Then, they wrote an evolution equation for the occupation probability p k m of site x m at time t k based on rules a), b) and c). Next, they interpreted the discrete evolution equation as a finite difference method for the approximation of a suitable continuous evolution equation. Let c = h ∆t , and define v = α + c the free velocity (the maximum statistical particle speed); the continuous evolution equation reads where p(t, x) is now a function of time t ∈ R + and space x ∈ R, derivatives are denoted by subscripts, and κ is a dissipative term, given by κ = 1 2 hv, which can be interpreted as the upwind artificial viscosity of the finite difference scheme.
Equation (1) can be reformulated as where u = −κp x + vp(1 − p) is the flux. However, the variable of interest here is the actual particle density ρ, defined as the number of particles per unit length. Since the number of sites per unit length ρ max = 1 h is the maximum particle density, the actual particle density is given by ρ(t, x) = ρ max p(t, x). Hence, we can reformulate (2) in terms of the density as follows, where If κ = 0, we recognize in (3), (4) the well-known Lighthill-Whitham-Richards (LWR) model of traffic flow on a single road [13], [19].
In this work we first derive a two-dimensional traffic flow model based on the same arguments; the main difference is that in 2D we cannot assume a preferred direction for the motion. The model is a (non-linear) advection-diffusion equation for ρ.
Then, we discuss its homogenization via asymptotic expansion. We obtain explicit formulae for the porosity and the permeability of the network, which appear as parameters in our equations, and of the advection field. Basing our arguments on the availability of a microscale model of traffic flow, we propose a method to efficiently evaluate the corresponding macroscopic parameters. Finally, we show how the distributed model is able to capture the main features of the microscopic flow.
2. Distributed traffic flow model on a 2D regular network. Let us consider a full two-dimensional lattice ( fig. 1, left). In order to proceed with an analysis similar to the one considered for the 1D case, let us assume that a particle at a given site ( fig. 1, on the right) can move to the east, west, north and south nearest sites with transition probabilities given by (respectively) α + , α − , β + , β − . We shall consider these probabilities as being constants. Using conditional probabilities for expressing the traffic rules a)-c), now gives the following evolution equation: Notice that the probabilities are subject to the constraint α + + α − + β + + β − ≤ 1, and 1 − (α + + α − + β + + β − ) is the probability that a particle in a free lattice does not move. Define α, δα, and β, δβ respectively the mean value and the deviation from mean value of transition probabilities: We have: +α − (p k m+1,n − p k m,n ) + β − (p k m,n+1 − p k m,n ) +α + p k m,n (p k m+1,n − p k m−1,n ) + β + p k m,n (p k m,n+1 − p k m,n−1 ) +α − p k m,n (p k m−1,n − p k m+1,n ) + β − p k m,n (p k m,n+1 − p k m,n−1 ), Since in our lattice gas model transition probabilities are constant, this can be regarded as a finite difference approximation of the following problem: where K is a diffusion tensor and v the bias in the motion trend of particles, given by We assume that h, ∆t have finite characteristic values, that determine K and v according to (9). As a particular case, notice that if we only allow forward horizontal displacements, setting the problem reduces to the 1D LWR (dissipative) traffic model (1). We point out that problem (8) can be recast in a conservative form: Again, we can rewrite our model in term of the particle density ρ(t, x) = ρ max p(t, x), where ρ max = 1 h 2 is the density of available sites (per unit area). We get Eq. (11) constitutes our distributed traffic flow model on a regular lattice with constant transition probabilities, which extends LWR model to the bi-dimensional case. It is a (nonlinear) evolution equation for the particle density, featuring a (linear) diffusive term, and a (nonlinear) advection term. We point out that, since K is proportional to h 2 ∆t , a purely hyperbolic conservation law would be recovered in the limit (h, ∆t) → 0 with constant ratio c = h/∆t. (i.e. K = 0 in (11)). In practice, h and ∆t have characteristic values (for instance, h could be of the order of magnitude of the lenght of a car, and ∆t = h/V max , where V max is the maximum velocity of cars). For this reason, we will assume that (9) are explicit formulae to be used also in the continuous model (8). Moreover, when two-dimensional traffic flows are considered, the unbiased situation in which transition probabilities to opposite directions are almost equal (i.e., δα ≃ δβ ≃ 0), is more likely to occur than in the 1D case of a single road, where the transition probabilities are biased according to the fixed direction and the orientation of the flow. In this case K = 0 even in the limit (h, ∆t) → 0. This is what happens for instance if we assume the scaling δα = h δα, δβ = h δβ (meaning that for h → 0 the transition probabilities will converge to α and β), for h, ∆t → 0 with constant ratio ∆t . In this case, the discrete lattice equations converge to the solution of the continuum problem (10) where In general, we expect the role of the diffusion term K to be more important in 2D.
3. The homogenized equations. In section 2 we have considered a regular lattice ( fig. 1) and constant transition probabilities. If we allow transition probabilities to vary on lattice nodes, by eq. (7) we can model networks different from being regular grids. As depicted in fig. 2, a given road network ( fig. 2(a)) can be represented by a lattice in which the only active nodes lie on roads. The other inactive nodes are such that the probability to jump into them is very low, possibly zero. In this case, particles move on the sub-lattice of active nodes only. Following a similar ideas to extend the continuous model to treat heterogeneous microstructures, let Ω be the considered 2D region, and Ω ǫ p ⊂ Ω be the subregion occupied by roads (the shaded area in fig 2(c)), where ǫ > 0 is the characteristic length scale of the road network. In this case, the quantities K and v computed by (9) are very low or even negligible in the "impervious" region Ω\Ω ǫ p , at least with respect to their values in the "pore" region Ω ǫ p . These "impervious" and "pore" (roads) regions, are easily detected by thresholding algorithms applied to images. For instance, under the hypothesis that vehicles move on each street only parallel to the street uniformly in both directions, one can easily provide the values of transition probabilities for any point of a given roadmap (in particular, transition probabilities will vanish outside the roads).
However, if the length scale ǫ of the single roads is too small with respect to the size of Ω, resolving the small scale is computationally too expensive. In such cases, the parameters K and v are rapidly oscillating according to the underlying road network geometry, which calls for suitable homogenized models, by means of averaging techniques.
As customary in homogenization theory, let us consider the domain Ω as ǫYperiodic, Y being a unit cell. The unit cell is thus representative of the microstructure of the network, for instance featuring highly permeable regions (the pores being the roads) included in low-permeable regions, associated to active and inactive nodes in a corresponding lattice gas framework, as discussed before. An example is given in fig. 2(a), where Y is the unit square [0, 1] 2 and the shaded subregion represented corresponds to the active nodes in the lattice gas framework, and by the ǫY -periodic network in 2(c), where the shaded region is Ω ǫ p . We are interested in deriving asymptotic flow models for small values of ǫ. For the sake of simplicity, assume that flow in the pore region is governed by problem (10) and consider the steady case, with homogeneous boundary conditions: Here b(x) is a datum (for instance due to the lifting of non-homogeneous Dirichlet boundary conditions). Next, consider the two-scale expansion: with u k and p k being Y -periodic in the second argument. By substituting these expansions in (12), applying the chain rule d dx = ∂ ∂x + ǫ −1 ∂ ∂y and collecting the terms of the same order in ǫ, we have and, in Ω ǫ p × Y , where we used p 0 (x, y) = p 0 (x) as a consequence of (13) 1 . Notice that (14) is a one parameter (x) family of Darcy problems in Y . For i = 1, 2, define w i (y), z(y) as the Y -periodic solutions of the cell problems where e i the i − th unit vector in R 2 , and let K 0 be the matrix whose columns are w i . In the same way, define v 0 as the Y -periodic solution of the cell problem Then, because of the linearity of problem (14), we have the following expression for u 0 : Notice that by (13) due to the Y -periodicity of u 1 . Hence, by averaging eq. (14) over the pore region and using (13) 2 we end up with the following homogenized model: where, denoting by |Y | the area of Y , we have set Equipped with the boundary conditions n ·ū 0 = 0 on ∂Ω, problem (17) represents the asymptotic homogenized flow model on Ω for ǫ → 0. At this stage we can draw some important observations. Remark 1. In spite of the presence of the nonlinear term, the two-scale analysis reveals that problem (12) admits a simple homogenization which has exactly the same structure of the non-homogenized model (11). Moreover, the homogenized model allows to upscale the original model using equations (15), (16) and (18) that provide the correct "averaged" valuesK 0 andv 0 of the permeability tensor and advection velocity.
Remark 2. The extension to the time dependent case is easy to carry out. In the same way, the model can be reformulated in terms of vehicle density. As previously discussed, we usually have a subregion Y p ⊂ Y occupied by roads and the complementary "impervious" region. Notice that according to this partition, φ = |Yp| |Y | will be referred to as the porosity of the network. Let us define the particle density A h is the maximal surface density of vehicles inside the cell Y , A h = h 2 being the occupational area associated to a single vehicle (and h is the lattice spacing parameter, the minimal length between two possible sites). Let us assume that ρ max is constant (or at least that its relative variations are small); then the time dependent homogenized model takes the following form, which has exactly the same structure as (11). network into elementary cells by means of a finite element triangulation ( fig. 3(b)), and consider the homogenized model (19) where the parametersK 0 andv 0 are computed in each cell by suitable techniques. We propose the following alternative methods to accomplish this task.
M1. In each element, we solve the cell problems (15), (16). Since the network is no more periodic, we consider no-flux boundary conditions instead of periodic boundary conditions. These cell problems can be solved by means of finite differences in which the computational grid is given by the matrix of pixels. In this case, each pixel corresponds to a site, so that the image resolution is dictated by the lattice parameter h. The thresholded gray image can be used to detect the pore region Y p (roads) where transition probabilities α ± m,n and β ± m,n at each pixel (m, n) have to be assigned; in the "impervious regions" Y \Y p , the transition probabilities will be set to zero. At this stage, the values of K and v at each pixel are known. Then, the associated homogenized quantities K 0 andv 0 can be computed by (15)- (18). In fig. 3(c) we see a representation ofK 0 , for a very simple situation in which the transition probabilities were perfectly isotropic. Here and in the sequel, we represent permeability tensors by means of ellipses, whose principal axes are given by unit eigenvectors scaled by the corresponding eigenvalues. The calculation of fig. 3(d) clearly shows that, even if the transition probabilities are isotropic in the pore region, the microscopic geometry of the network may affect the homogenized perme-abilityK 0 , which can be anisotropic. This behavior is emphasized in 3(d) by the stretched ellipses aligned with the Seine river which acts as a barrier for the flow. In general, we observe that the ellipses are stretched according with the underlying microstructure of the road network. This shows how the microscopic features of the network can affect the macroscopic flow, even for unbiased particles.
M2. Assume that the average values K p and v p of K 0 and v 0 are known on the pore region Y p . For instance, if the micro-model is a lattice gas model, Y p is the subset of active sites and such values may be provided by (9). Then, since K 0 and v 0 vanish on Y \Y p , thanks to (18) we can approximatē on each cell.
Method M2 allows to make use of external micro-models to compute the parameters of the upscaled equations. This is helpful since M1 has a practical drawback: it requires a detailed knowledge of transition probabilities at the microscale. The latter is a rather strong assumption, and even if some statistical data are available at a macroscopic level, the same situation is unlikely to occur for the microscale.
With little information about the transition probabilities, the computed parameters may not be representative of the underlying phenomena.
In the next section, we will present a general statistical algorithm to compute K p and v p on the pore region Y p (i.e. the local microscopic road network) of each cell Y using a microscopic flow model, such as those based on cellular automata (CA), in which the behavior of a population of particles is simulated accounting for individual interactions. Hence, by (20) the method M2 allows to plug the computed values of K p and v p into the distributed model (19).

Computation of averaged parameters from CA microscopic models.
In this section we propose a general strategy to extrapolate the mean values of pore permeability K p and drift velocity v p using a microscopic model M Y defined over the active or "pore" region Y p of each cell Y of a coarse tessellation of the domain Ω. The "pore" region Y p is actually a network: we will refer to the sub-network contained in the cell as the micro-network. Since the global network is no more periodic, we cannot rely on homogenization formulae given in section 3. Moreover, we are interested in providing a general method that allows the use of the generic traffic model M Y on the micro-network to extrapolateK 0 ,v 0 , i.e. we aim to follow the approach M2 of section 3.1.
Let us observe that the pore permeability K p and drift vector v p depend on the motion of the particles in the micro-network, but not on the traffic intensity. The idea is thus to run several instances of the microscopic model M Y for a single particle, then averaging the outcomes. This approach is based on the linearity of the lattice gas model for a single particle. In fact, let us rewrite eq. (11) considering a single particle. In this case, the evolution equation of the occupation probability becomes that is which is the discrete finite differences counterpart of the linear equation with the same K and v of eq. (11). In other words, the nonlinearity in (11) is only due to interactions. Now, looking at K and v respectively as the covariance matrix and the mean value of a random walk, we propose to use a Monte Carlo method [23,20] to estimate the K p , v p respectively with the sample covariance and sample average of a collection of simulated random paths by the micro-model M Y . Possible instances of micro-models are CA algorithms. In our experiments, we considered for M Y an extension of the well known Nagel model [16,17], adopting the following simple laws for a single particle. L1 At every time step the particle accelerates until a given maximal velocity is reached (there is no interaction with other particles). L2 The particle always decelerates when approaching road junctions and gets into the crosses at a given minimal velocity. L3 Suitable "normalized fluxes" F i are defined at road junctions, assigning the probability of choosing the i-th outgoing arc.
Obviously, the trajectories obtained by this extended Nagel algorithm are not pure bi-dimensional random walks. In order to reduce the error, we use a multistart Monte Carlo algorithm based on the following steps: S1 N trajectories x k i , i = 1, . . . , N , k = 1, . . . , T i + 1, are created in Y p using the micro-model M Y . T i + 1 is the exit time of the i-th trajectory. The starting point of each path is uniformly distributed on all the possible entries of the map Y p . S2 The sample mean m i and the covariance S i of the increments of each trajectory are computed:

FABIO DELLA ROSSA, CARLO D'ANGELO AND ALFIO QUARTERONI
S3 The drift velocity v p and the permeability tensor K p of the element are computed as the weighted sample means over the trajectories,

S4
The parametersv 0 andK 0 of the distributed model (19) are computed on each cell Y by (20) using the cell porosity φ.
The choice of the fluxes F i inserts a bias in the motion of a particle. In particular, we distinguish the case in which particles have no preferred direction from the case in which a preferred direction e exists, meaning that particles will statistically prefer directions that promote motion along e. Notice that the local network geometry may have an impact on the drift vector, giving rise to a macroscopic bias different than e. In the first case, at each cross we require F i to assign the same probabilities to all the possible outlets. In the second case, those probabilities will not be uniform and have to be computed. To set up a statistical law describing the preference of drivers for direction e, we proceed as follows. Consider a given node, and let θ i be the angle between the line passing through the initial and the final point of the i-th outgoing arc and e. Then, we associate to this arc the weight where the constant γ is a tuning parameter of the biasing. Finally, we assign the probability F i of choosing the i-th arc by normalizing each weight w(θ i ). We point out that suitable boundary conditions are needed in step S1, to be used when the particle hits the boundary of the cell micro-network. We consider free-flow boundary conditions, i.e. if a particle reaches the boundary, the simulation stops. Correspondingly, we stop the algorithm when all the particles have reached the border of the map. Indeed, at step S3 the upscaled data are computed as weighted sample means, using the time of permanence of the car in the map (i.e. the number of iterations before reaching the border of the map) as weight.

5.
Simulations. In this section we consider an application of our upscaling techniques to real urban networks. First of all, the estimation algorithm will be validated on simple networks. Then, it will be employed to compute the parameter of the distributed model of traffic flow over an extended region, and the results of the distributed model will be compared with a fully microscopic simulation over the extended region.
5.1. Sensitivity analysis with respect to γ. As a first case study, let us consider the influence of the parameter γ on the upscaled parameters. To do that, we performed 20000 single-particle Monte Carlo simulations on the simple network in The straight lines in fig. 4(a) are the computed drift velocity v p , the axes of the ellipses represent the principal directions of the computed permeability K p . From fig. 4(c), we see how the biasing parameter γ affects the effective drift velocity. On the other hand, from fig. 4(b), we conclude that γ has little influence on the permeability tensor (about 10%). The estimation of the permeability tensor is thus almost independent of the bias in the particles behavior.

5.2.
Sensitivity analysis with respect to N and convergence of the algorithm. As a second benchmark, we study the number N of required Monte Carlo simulations for the estimation to be good. We consider the standard deviation of the computed quantity as a performance index, in agreement with the literature (see [21] or [9]). In particular, we performed 100 runs using increasing N , and evaluated the standard deviation of v 0 and K 0 . For the estimation of the drift velocity v p we have used bias parameters given by γ = 0 and γ = 20, whereas for the permeability tensor K p we have considered, according to the previous observations, only γ = 0. The results for the estimation are shown in fig. 5. From these results, we can state that the convergence of the algorithm is satisfactory using N ≥ 10000.

5.3.
Comparison between upscaled and microscopic models. We are now ready to apply our techniques to a more complex and extended network. The area of application is the urban tissue of Le Havre. This is today considered the "Paris' harbour", even though it is 150 Km far away from the French capital. We mention that the role of "peripheral" urban areas at length scales of 100Km, either poorly integrated banlieues or external well-developed and independent centers, was a keynote of the Le Grand Pari(s) programme [11]. The considered road network is shown in fig. 6: it spans over a 30Km×30Km square domain Ω. In order to validate our upscaling techniques, we propose a comparison between an accurate microscopic CA model, based on the extension of the Nagel model to road networks, and the upscaled model. For the sake of simplicity, we only describe the modifications with respect to the standard Nagel model and to rules L1-L3 of section 4. Consider the one-dimensional lattice with spacing h between sites, let ∆t be the time step, let v k i and g k i be respectively the velocity (number of spanned sites per ∆t) of the ith vehicle and the gap (number of sites) between the ith vehicle and the next one at time k. The extended Nagel model evolves according to L2, L3 and the following modified L1 rule (to account for interactions): (towards Paris). The time step of the CA algorithm was 2 seconds. For the sake of simplicity, the sites of the Nagel model were identified to the pixels of the image, each pixel corresponding to 6 meters. As boundary conditions, we let car exiting the network disappear; moreover, we assign a flux of 10 cars per second entering the domain by the nodes belonging to the upper half edge on the west side of the map, Γ * in . As initial conditions, we distribute 62500 cars on a surface of 45 Km 2 on the north-western 3 Km×15Km subregion of map, corresponding to the half-saturation of roads in that area (50% free sites on all roads).
Next, we computed the resulting cell average density ρ CA , as a piecewise constant function computed by counting the vehicles in each cell Y divided by |Y | at each time step. ρ CA will be our "reference" density distribution to test the upscaled model. Now let us consider the upscaled model. The map was subdivided in 10×10 cells as in fig. 6; each cell Y was then processed following steps S1-S4 to extract cell values ofK 0 andv 0 to be used in the upscaled model. The biasing parameters were the same as in the global microscopic simulation. According to the findings of sections 5.2 and 5.1, we choose N = 20000. Results are presented in fig. 6, where we reported on the right the magnitude of the permeability fieldK 0 and the drift velocityv 0 computed in each cell. Again, the geometry of the street pattern influence the macroscopic drift vector. This fact is emphasized by the cells near the river, where the computed drift vector may significantly deviate from e. The porosity φ was computed as the ratio between road pixels and non-road pixels in each cell Y . Since pixels are also the sites of the Nagel model, we set ρ max = φ in each cell: this is the maximal car density at which all available sites on each road are occupied so that vehicles can not move at all.
OnceK 0 andv 0 have been computed, it is possible to set up the upscaled distributed model (19) of traffic flows in the considered region. In our specific case, we seek the vehicle density ρ(t, x), t ∈ R + , x ∈ Ω, satisfying where Ω is the considered square domain, Γ in := {x ∈ ∂Ω :v 0 (x) · n(x) < 0} is the inflow boundary, Γ out = ∂Ω\Γ in is the outflow boundary, n is the outward unit normal vector, f is the total flux, and ρ 0 = ρ 0 (x) is the initial condition, f in = f in (x) is the car flow rate at the inflow. We will assume that a positive lower bound for the eigenvalues of the symmetric tensorK 0 exists. In this case, problem (25) is parabolic, and we may resort to backward Euler or Cranck-Nicolson time advancing scheme with continuous linear finite elements in space. Equation (25) will be generally dominated by advection (as results from section 2) so that the numerical scheme will need stabilization with respect to transport. For a survey of different stabilization techniques we refer to [18]. In our application we choose an interior penalty stabilization of the convective terms. For more advanced numerical schemes based on interior penalty techniques, which are also robust with respect to the vanishing permeability case, we refer to [7], [24]. We also notice that for purely hyperbolic problems, the solution can exhibit shocks, so that alternative schemes have to be chosen, such as finite volumes. We consider a collection T h of disjoint elements K with characteristic size h > 0, such that Ω = K∈τ h K, diam(K) < h. Of course, in our applications, T h may coincide with the tessellation introduced to upscale the micro-model. However, it can also be finer for computational purposes. Let V h = {v h ∈ H 1 (Ω) : v h | K ∈ P 1 (K), ∀K ∈ T h }. Consider the following bilinear and linear forms, where c 0 = 1 ∆t −∇·β is assumed to be positive, F h is the set of all internal edges, · is used to denote the jump of functions across edges, and γ β > 0 is the stabilization parameter.
Then, the approximate solution ρ n h at time t n = n∆t is computed as follows: For k = 0, 1, . . . until convergence, compute ρ n,k+1 For n = 0, 1, . . ., set ρ n+1 Here, Π V h is the L 2 orthogonal projection on the finite element space V h , andv 0,h is the orthogonal H div projection on [V h ] 2 . Projection is needed since the considered finite element scheme requires β to be continuous across edges.
The results of computations and the comparison between the microscopic CA model and the upscaled model are reported in fig. 7. In each column, we have displayed some snapshots of ρ CA on the left, and the corresponding distribution ρ n h on the right. For the simulations, we considered the tessellation T h obtained by splitting each square cell Y (used for the upscaling algorithm) in two triangles. In fig. 7(b), we have highlighted the inflow Γ in (solid and dashed red lines); moreover, according to the CA simulation, the entering flux f in is non-zero only on the subset Γ * in represented by the solid line. First, let us comment on the average cell density ρ CA computed by the accurate microscopic model. The evolution is characterised by the migration of cars according to the biasing direction e. Correspondingly, after a certain time, the car distribution reaches an equilibrium in which the flux given by the entering cars form Γ * in is balanced by the flux of cars exiting the map. The "barrier action" of the river, which is not easily crossed due to lack of connecting roads (bridges), is clearly seen. In fact, the density front advances in the upper part of the map rather than in the lower part.
The simulation with the upscaled model provides similar results. We observe that the high-density regions are the same for both models. The front expansion is similar, and the information about the low permeability of cells containing the river clearly affects the distributed model: in fact, the front does not invade the lower part of the map and still prefers to expand in the north-west side.
The main differences are related to small density regions. In fact, thanks to the log scale used in fig. 7(b,d,f), we can appreciate the few cars distributed in the lower part of the map according to the CA micromodel. The upscaled model is not able to describe such very fine details. However, the behavior of the upscaled model is satisfactory for the regions where the density is not very low. Moreover, the simulation of the upscaled model resulted to be 20 times faster than that of the CA microscopic model. Figure 7. Evolution of the particle density ρ CA computed with the accurate microscopic model and of the particle density ρ h computed with the upscaled distributed model. On each column, the left plot displays ρ CA and the right plot reports ρ h . Rows correspond respectively to the initial time, to t =4h30 and t =9h. Data are represented with interpolated surfaces and colors on the left column, and with colors on the right column. Colors are graded logarithmically, due to the wide range of density. In subfigure (b) Γ * in and Γ in \Γ * in are highlighted respectively by solid and dashed red lines. 6. Conclusions. In this work, we considered traffic flows model on extended twodimensional regions, representing broad metropolitan areas. Starting from lattice gas models, we derived an upscaled model to be used when the length scale of the considered region is huge with respect to the characteristic lattice spacing parameter. It is written as a semilinear Burgers equation, in which the viscous and transport terms depend on the transition probabilities of the lattice model.
Next, we extended the range of application of the upscaled model to arbitrary networks by two-scale averaging techniques, following the homogenization approach. The concepts of permeability and porosity of the network were naturally introduced. We also discussed different methods to compute the relevant parameters to be used in the upscaled model, starting from widely available data (e.g. images). Those methods are essentially based on a partition of the region in elementary cells; cell problems are then solved to compute the parameters to be used in the upscaled model.
In particular, a very general method to compute the upscaled parameters using an arbitrary micromodel for each cell and Monte Carlo techniques has been discussed and accurately tested on real cases.
Finally, computational results from the upscaled model have been compared with those of a fully microscopic simulation in the region of Le Havre, France. The comparison shows that the car densities found using the two models agree in those subregions where cars are mostly concentrated. On the other hand, fine details about low density subregions are lost. However, the most relevant features of vehicle distribution are captured. In particular, it seems that the upscaled parameters are able to correctly identify high-permeability and low-permeability regions, that affect traffic flows as expected. Potential application of the described model is the multiscale simulation of traffic flow on extended regions, where highways and principal roads could be represented by macroscopic traffic models on networks, and the matrix of smaller roads could be modelled by resorting to our averaging approach.