Numerical approximation of continuous traffic congestion equilibria

Starting from a continuous congested traffic framework recently introduced in [Carlier, Jimenez, Santambrogio, 2008], we present a consistent numerical scheme to compute equilibrium metrics. We show that equilibrium metric is the solution of a variational problem involving geodesic distances. Our discretization scheme is based on the Fast Marching Method. Convergence is proved via a $\Gamma$-convergence result and numerical results are given.


Submitted on 19 May 2009
HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt età la diffusion de documents scientifiques de niveau recherche, publiés ou non, emanant desétablissements d'enseignement et de recherche français ouétrangers, des laboratoires publics ou privés. 1. Introduction. Congestion is an important issue in real-life applications such as road or communication networks. In the early 50's, Wardrop [17] considered the situation where a large number of vehicles have to go from one location to another, connected by a finite number of different roads. Each vehicle has to choose a pathway along the roads to minimize some transportation cost which depends not only on the chosen pathway but also on the number of vehicles (the traffic) along it. Naturally, shorter or wider roads are preferred by the vehicles. Nevertheless, the transportation "cost" on a given road depends increasingly on the traffic. To avoid the congestion some vehicles may choose a "longer" pathway to minimize the cost. Wardrop gave a minimal stability requirement for transportation strategies: the cost of every actually used pathway should be equal or less than that which would be experienced by a single vehicle on any other roads. In particular there is an equilibrium concept: all actually used pathways have the same cost, i.e. they compensate in term of congestion their differences given by length and other conditions. Beside the equilibrium concept there is a minimality concept as well: those pathways are minimal among all the possible ones. In other words, a Wardrop equilibrium is a situation where every traveler uses only shortest paths from his source to his destination, given the overall congestion pattern resulting from the individual strategies of all the users of the road network. This natural equilibrium concept, somehow similar to Nash's equilibrium, has been very popular since its introduction.
Few years later, Beckmann et al. [3] discovered that Wardrop equilibria can be obtained as minimizers of a certain convex optimization problem. This result is very useful for proving theoretical existence, uniqueness and stability results. The convex minimization formulation (or some dual one) can also be used to compute numerically traffic equilibria as long as the complexity of the problem is tractable which, in practice, often excludes the case of large networks. On the contrary, in the present paper, we consider a continuous case, where there is no network and every path is allowed. Our work is based on a continuous traffic equilibrium formulation recently introduced in [8]. We will see how this problem can be dealt without huge computational costs.
The study of Wardrop equilibria has mainly been restricted to the finite-dimensional case where the road network is a given finite graph. In [8], there is no network and congestion effects are captured by a metric that depends increasingly on the traffic intensity which gives a natural framework to study optimal transportation problems with congestion effects. It turns out that the problem addressed in [8] can also be interpreted as a continuous traffic equilibriumà la Wardrop. The equilibrium problem is then formulated in terms of the traffic intensity or the metric itself. As shown, in section 4, the equilibrium metric is obtained by minimizing some convex functional that is the difference of two terms. The first one is a standard convex integral functional. The second one is a weighted sum of geodesic distances between a given source and a given destination as a function of the metric. Building upon this continuous minimization problem, we introduce a discretized minimization scheme based on Sethian's FMM [13] (Fast Marching Method). We prove convergence of the discretized problems to the continuous one, obtaining a Γ−convergence result which is interesting by itself. It is based on the convergence of the FMM scheme to the viscosity solution of the Eikonal Equation |∇U| = ξ in the case ξ ∈ C 0 and on some semicontinuity properties coming from the comparison between different concepts of this equation for ξ ∈ L p . We then solve numerically the discretized problem by an iterative subgradient descent method. Each iteration of the discretized problem requires the computation of an element of the subgradient of the geodesic distance with respect to the current metric. In our companion paper [4], we show how to compute such a subgradient in the same routine of the FMM algorithm and develop other applications of this method.
The continuous traffic congestion model is introduced in section 2. In section 3, traffic equilibria are defined and characterized by means of some convex optimization problem. A more tractable dual formulation is given in Section 4. In section 5, we consider a discretization of this dual problem, prove that this discretization Γconverges to the continuous problem and give a consistent numerical scheme to solve the discrete problem. Numerical results are presented in section 6. 2. The model.

Notations and definitions.
Given a separable and complete metric space X, M + (X) and M 1 + (X) denote respectively the set of positive and finite Borel measures on X and the set of Borel probability measures on X. If X and Y are separable metric spaces, µ ∈ M 1 + (X) and f : X → Y is a Borel map, f µ denotes the push forward of µ through f i.e. the element of M 1 + (Y ) defined by f µ(B) = µ(f −1 (B)) for every Borel subset B of Y .
In the sequel, L d denotes the d-dimensional Lebesgue measure. If µ and ν are in M + (R d ) then dµ dν denotes the Radon-Nikodym derivative of µ with respect to ν. We write µ ν to express that µ is absolutely continuous with respect to ν, in which case, slightly abusing notations, µ is identified with the Radon-Nikodym derivative dµ dν . The data of our problem are Ω (modelling a city, say), an open bounded connected subset of R 2 with a Lipschitz boundary, two probability measures, µ 0 and µ 1 in M 1 + (Ω), giving respectively the distribution of residents and services in the city Ω and a transport plan γ i.e. a probability measure on Ω × Ω having µ 0 and µ 1 as marginals. The transport plan γ models the travelers' everyday movements, more precisely for A and B Borel subsets of Ω, γ(A × B) represents the fraction of the travelers' population that commutes from a location in A to a location in B.
Introducing congestion naturally leads to consider spaces of paths, lengths of such paths and sets of probability measures on sets of paths. From now on, we shall denote: Roughly speaking, a transportation strategy is a probability Q over paths (Q(Σ) represents the number of travellers that use a path in Σ) which is compatible with the transport plan γ: If ξ ∈ C 0 (Ω, R + ), it is proven in Lemma 2.6 in [8] that the map σ → L ξ (σ) is l.s.c. on C hence Borel for the uniform topology. Since 0 ≤ L ξ (σ) ≤ ξ ∞ l(σ), L ξ ∈ L 1 (C, Q) as soon as: which means that the average (with respect to Q) length is finite. If (2.1) holds then L ξ ∈ L 1 (C, Q) for every ξ ∈ C 0 (Ω, R + ). If ξ ∈ C 0 (Ω, R), let us write ξ = ξ + − ξ − with ξ + and ξ − the positive and negative part of ξ. Since L ξ = L ξ+ − L ξ− then L ξ is Borel and in L 1 (C, Q) provided (2.1) holds. Any Q ∈ M 1 + (C) induces a traffic intensity which is a nonnegative measure on Ω, intuitively the traffic intensity of a subregion of Ω is the cumulated traffic of this subregion (or the total amount of mass travelling through the subregion). For simplicity, in the next definition, to avoid the case of measures taking the value +∞, we restrict to the case where Q satisfies (2.1).
Let us remark that the total mass of i Q is the average length with respect to Q.
In the sequel, we shall need to extend the definitions of c ξ and L ξ to the case where ξ only belongs to some L p space. This is possible when p > 2 (we refer to [8] for details), by proceeding as follows. Let us assume that p > 2 and ξ ≥ 0, ξ ∈ L p (Ω), then the following is proved in [8].
This allows to extend the definition of the geodesic distance c ξ to the case ξ ≥ 0 and ξ ∈ L p (with p > 2) as follows: This definition really extends that of geodesic distance since when ξ is continuous and non-negative, then c ξ = c ξ (see [8] for details). Section 5.3 shows that this definition coincides with a different one, given in terms of pointwise a.e. subsolutions of an Eikonal Equation.
It remains to extend the definition of L ξ . To do that, let us define for every q ∈ [1, +∞], the set: where in the previous definition we have slightly abused notations to intend that i Q L 2 and that the density of i Q with respect to L 2 is in L q . The requirement will naturally appear as a consequence of our traffic congestion modelling detailed in the next paragraph (roughly speaking, in our model, congestion effects are so strong that when (2.5) is not satisfied then the total transportation cost is +∞). If p > 2, ξ ∈ L p , ξ ≥ 0, (ξ n ) n is a sequence of non-negative continuous functions that converges to ξ in L p and Q ∈ Q p * (γ) (p * being the conjugate exponent of p), then the following holds (see [8]): • (L ξn ) n converges strongly in L 1 (C, Q) to some limit which is independent of the approximating sequence (ξ n ) n and which is, again, denoted as L ξ .

Traffic congestion.
Traffic congestion modelling has to capture the fact that travellers spend more time (or more effort) in regions where the traffic is dense i.e. where the traffic intensity is high. A natural way to do so (assuming for a moment that i Q L 2 with a density with respect to L 2 still denoted i Q ), is to associate to i Q a metric of the form g(., i Q (.)) where g is a given congestion function. We shall always assume the following on the congestion function g: • there exists strictly positive constants a and b and a constant α ∈ (0, 1) such that Let us set q := 1 + α and q * = 1 + 1/α its conjugate exponent (note that q * > 2 so that one can define c ξ as in the previous paragraph whenever ξ ∈ L q * and ξ ≥ 0) and define Q q (γ) by (2.4). Let us also remark that every Q ∈ Q q (γ) satisfies (2.1). For Q ∈ Q q (γ), the congested metric associated to Q is well-defined and given by: Note that our assumptions on g imply that ξ Q ∈ L q * whenever Q ∈ Q q (γ). Since q * > 2, by the results recalled in paragraph 2.2, c ξ Q is a well-defined continuous function, L ξ Q ∈ L 1 (C, Q) and (2.6) and (2.7) hold for ξ = ξ Q .

3.1.
Equilibria. An equilibrium is a transportation strategy Q such that the congested metric ξ Q is well defined (i.e. Q ∈ Q q (γ) with q = α + 1 ∈ (1, 2)) and such that, given the congestion resulting from ξ Q , travellers only use shortest paths i.e. geodesics. This stability requirement on the transportation strategy (all used routes connecting x to y have the same commuting time and a shorter commuting time than unused routes) was first introduced by Wardrop [17] in the case of a finite number of roads. The generalization to our continuous setting then reads as: In other words, Q ∈ Q q (γ) is a Wardrop equilibrium if and only if Q gives full mass to the set of geodesics for the metric ξ Q defined by (2.9).
Of course, a necessary condition for the existence of a Wardrop equilibrium, in our model where congestion effects are very strong, is that Q q (γ) is nonempty. As shown in section 3.2, this condition is also sufficient for existence of equilibria. At first glance however, the condition Q q (γ) = ∅ may seem difficult to check a priori if we only know γ. Yet, in the case when there is a finite number of sources and destinations (which is basically the situation considered in our numerical simulations), this condition is automatically fulfilled.
Proof. It is of course enough to prove the claim when γ is a single Dirac mass γ = δ (x0,y0) and up to a change of coordinates we may assume x 0 = (0, 0) and y 0 = (1, 0). Now let β 0 > 0 be such that for all β ∈ [0, β 0 ], the curve σ β defined below lies in Ω (this is where the convexity assumption simplifies the argument eventhough it is not essential): Denoting by T 0 the triangle with vertices (0, 0), (1/2, 0), (1/2, β 0 /2) and performing the change of variables (β, t) → σ β (t), the first integral in the expression above equals 1 β 0 T0 ϕ(x, y) 1 x 1 + y 2 x 2 dxdy and a similar expression gives the value of the second integral. Now remarking that 1 + y 2 x 2 ≤ 1 + β 2 0 on T 0 and that 1/x ∈ L p (T 0 , L 2 ) for all p ∈ [1, 2) we get the desired result. Remark 1. Some techniques from incompressible fluid mechanics, as Yann Brenier pointed out to us, may be used to prove the existence of finite-energy Q s in the case where the transport plan γ is fixed but the marginal measures are not discrete. For instance, if µ = ν = L 2 and γ is a transport plan from µ to ν, one can find a measure Q which is compatible with γ (i.e. (e 0 , e 1 ) # Q = γ) and incompressible, in the sense that it realizes at any time t the same density (e t ) # Q = µ, and this measure may be chosen concentrated on K−Lipschitz curves. In this case one can see that i Q ≤ KL 2 ∈ L ∞ . Composing with diffeomorphisms (possibly time-dependent) one can generalize the same idea to other pairs of measures.

3.2.
Existence and characterization of equilibria. Our study of equilibria relies on the following convex optimization problem: where the function H is defined by By the same arguments as in [8], one can prove the following.
then Q is an equilibrium if and only if Q solves (P).
Actually, the problem addressed in [8] is slightly different from the one considered here, but all the proofs can be straightforwardly adapted. Indeed, in [8], the transport plan γ is not given a priori, only its marginals are, a situation slightly more general than the present one. Also, in [8], the congestion function g only depends on i but it is easy to check that the results of Theorems 3.2 and 3.3 extend to the case where g also depends on x under the assumptions of the present paper.
From now on, we assume that Q q (γ) = ∅ (q = 1 + α ∈ (1, 2)). In this case, equilibria exist and finding them amounts to solve (P). Since H(x, .) is strictly convex, if Q 1 and Q 2 are equilibria, then necessarily i Q1 = i Q2 . In other words, equilibria are not necessarily unique but they all induce the same intensity or equivalently the same metric ξ. The next section shows that finding the equilibrium metric amounts to solve another optimization problem which is dual to (P). This dual problem turns out to be easier than (P) to be numerically handled. By our assumptions on g, one has H * (x, ξ) = 0 for every x ∈ Ω and ξ ≤ ξ 0 (x), and Let us recall Young's inequality: and that inequality (4.1) is strict unless ξ = g(x, i) ≥ ξ 0 (x). In particular, for Q ∈ Q q (γ), we have the identity and 3) Let us now define the functional and consider: Theorem 4.1. Let us assume that Q q (γ) = ∅ (q = 1 + α ∈ (1, 2)), then min(P) = max(P * ) (4.6) and ξ ∈ L q * solves (P * ) if and only if ξ = ξ Q for some Q ∈ Q q (γ) solving (P).
In the sequel, we numerically approximate the unique equilibrium metric ξ Q by a descent method on (P * ). One can recover the corresponding equilibrium intensity i Q by inverting the relation ξ Q (x) = g(x, i Q (x)).
5. Discrete Algorithm. The traffic intensity that we are looking for, which is both the optimal one and the unique equilibrium, is determined by solving the dual problem in the variable ξ = H (i). This section describes and justifies a numerical approximation of the optimal ξ in the dual problem. This procedure (approximating the solution of the dual problem instead of the primal one) is classical as far as discrete Wardrop problems on networks are concerned (see [1] and the references therein). Obviously our numerical procedure passes itself through a discretization, but this has nothing to do with the network formulation of discrete Wardrop problem. The discretization of ξ is indeed done on a square lattice and we prove the convergence of our scheme to the solution of the continuous problem.

5.1.
A discrete functional via Fast Marching. This section is concerned with the discrete formulation of the dual problem. After setting the notations in this discrete setting, the convexity of the discrete functional, with respect to the discrete metric, is proven. The convergence of this discretization to its continuous counterpart is proven in the next section.
The metric is discretized as a vector ξ = (ξ i,j ) i,j ∈ R N where ξ i,j represents the value of the metric at a point x i,j = (ih, jh) ∈ Ω of a square lattice. The number of grid points is N = h −2 . The first term of the energy functional J defined in equation (4.4) is easily discretized as h 2 i,j H * (x i,j , ξ i,j ). The second term of the energy J requires the discretization of the geodesic length c ξ . This can be performed by computing the solution of a discretized Eikonal equation. For a fixed source S 0 , the geodesic distance to S 0 is denoted as U ξ (x) = c ξ (S 0 , x). The function U ξ is the viscosity solution of the Eikonal equation The vector (U i,j ) i,j is a discretization of the values U ξ (x i,j ) of the geodesic distance on the lattice vertices, where the dependence on ξ is dropped to ease notations. A discrete version of (5.1) is solved in order to compute U = (U i,j ). For the Eikonal equation, classic finite difference schemes tend to overshoot and are unstable. Rouy and Tourin [11] showed that a correct scheme for approximating the viscosity solution for of (5.1) is given by the following first order accurate scheme:

The Fast Marching Method (FMM) is a numerical method introduced by Tsitsiklis
[15] and Sethian [12], see also [13]. This method efficiently solves the isotropic Eikonal equation (5.2) on a cartesian grid. Values of U may be regarded as the arrival times of wavefronts propagating from the source point S 0 with velocity 1/ξ. The central idea behind the FMM is to visit grid points in an order consistent with the way wavefronts propagate, i.e. with the Huygens principle. The numerical complexity of the FMM is O (N log(N )) operations for a grid with N points.
Let us call c h ξ (S 0 , T ) = U(T ) the value of this solution, computed with FMM, at the target point T when the vanishing boundary datum is fixed at S 0 and the metric is ξ. The dependence with respect to the discretization parameter h will be sometimes omitted in this section, since the lattice is fixed.
The functional to be minimized is the following where the weights γ α,β represent the coupling on the set of pairs sources S α -targets T β and α,β γ α,β = 1.
Lemma 5.1. The functional J h is convex.
Proof. The only non trivial point is proving the concavity of the FMM term. It is sufficient to prove that for fixed S 0 and T the function ξ → U ξ (T ) is concave in ξ and, thanks to the homogeneity, it is sufficient to prove super-additivity. We want to prove the inequality U ξ1+ξ2 ≥ U ξ1 + U ξ2 .
Thanks to the comparison principle of Lemma 5.2 below, it is sufficient to prove that ξ 1 + ξ 2 ≥ D(U ξ1 + U ξ2 ). This is easily done if we notice that the operator D is convex (as it is a composition of the function (x, y) → x 2 + y 2 , which is convex and increasing in both x and y, and the operator D x and D y , which are convex since they are produced as a maximum of linear operators) and 1−homogeneous, and hence it is subadditive.
Proof. Let us suppose at first a strict inequality ξ < η. Take a minimum point for U η − U ξ and suppose it is not the fixed point S 0 . Computing D and using the previously noticed sub-additivity, we have Yet, at minimum points we should have D(U η − U ξ ) = 0 and this proves that the minimum is realized at S 0 , i.e. that U η − U ξ ≥ 0.
To handle the case ξ ≤ η without a strict inequality, just replace η by (1 + ε)η and notice that the application η → U η is continuous.

5.2.
Description of the algorithm and convergence of its iterations. The discrete function J h , defined in equation (5.3) is convex and can be minimized using standard methods from convex optimization. The first term of the functional J h is differentiable, while the second may not be. Hence a subgradient descent algorithm can be used for the discrete optimization problem.
In the following, the dependence of H * on the position is omitted to ease notations and we write H * (x, ξ) = f (ξ). The functional to optimize is written in the form for a regular function f : R + → R + and a convex, but in general not C 1 , function K which is not separable w.r.t. to the variables ξ i,j .
The subgradient method corresponds (taking for simplicity ξ 0 = g(., 0) = 0) to the following iterative scheme where v (k) ∈ ∂K(ξ (k) ) is a vector in the subdifferential of K at the previous point ξ (k) and ρ k is a suitable sequence of steps.
A well-known result on subgradient algorithms (see for instance [5]) states that the sequence (ξ (k) ) k converges to the unique minimizer of J (uniqueness comes from strict convexity) provided the following two constraints are enforced • The step sizes ρ k satisfy k ρ k = +∞ and This condition is satisfied for a sequence of steps such as ρ k = 1/k. • The sequence (w (k) ) k stays bounded, which can be enforced by slightly modifying the penalty function f . Since K is 1−homogeneous, it is Lipschitz continuous and hence its subgradients are bounded. However, the original penalty f is not necessarily Lipschitz, but, one can consider instead the following Lipschitz functionf defined as It happens that, if t 0 is sufficiently large, the minimizers ξ of the modified energyJ and of the original one are the same and satisfy ξ i,j ≤ t 0 . This choice of t 0 may always be performed as far as f is superlinear at infinity. A computation of an element v (k) of the subdifferential of K is required for the gradient descent. Indeed, in [4] a recursive method that computes such a vector in parallel to the iteration of the FMM is detailed. This methods produces a vector which is the gradient of U w.r.t. ξ at differentiability points and a vector in the subdifferential at all other points. The complexity of this algorithm depends on the complexity of the original FMM, as it visits the points in the same order. Yet, since it must compute a N −dimensional gradient at each iteration, the result has complexity O(N 2 log(N )) for a grid of N points.
Notice that, in the continuous framework, the gradient of U ξ (T ) with respect to ξ is concentrated along geodesics joining T to S 0 . The discrete gradient which is computed with the method of [4] is supported on the grid points explored by the front propagation of the FMM in order to reach the point T (see figure 1, which gives sort of fattened geodesics).
The method of [4] allows us to compute with a fast algorithm subgradients of a single term ξ → U ξ (T ) computed by FFM. Since the functional K is a conic combination of terms of this kind, it requires to use this algorithm several times and combines the resulting subgradient in a subgradient of K (we recall that the sum of subdifferentials is always included in the subdifferential of the sum and, here, since every term in K is continuous, we actually have the converse inclusion).
The convergence of the subgradient descent (5.4) ensures that one can find with a fast iterative algorithm the minimizer of the discrete functional J h . Section 5.3 relates this discrete solution to the original continuous solution and shows that they are close in a suitable sense. Section 6 shows numerical simulations that use the descent algorithm 5.4 for various examples of penalty functions H * (x, ξ).

Γ−convergence of the discrete functionals to the continuous one.
This section considers the problem of the convergence of minima and minimizers of the discretized problems (5.3) to the original minimal value and minimizer of the functional (4.4). This allows to say that, as far as both the discretization step of the lattice is small and the number of iterations of the subgradient descent algorithm is large, one gets a good approximation of the desired optimal ξ.
The result is presented through a Γ−convergence proof, which implies the convergences we want (see [10]). Both the discrete problems and the continuous one  Figure 1. Examples of the subgradient computation using our algorithm presented in [4]. On the left, an element of the subdifferential of the geodesic with respect to the metric. The metric is constant, and the geodesic is a straight line. In the middle, we chose a non constant metric ξ(x) = 1/(1.5 − exp(−d (O, x))), where O is the center of the image. On the right, an element of the subdifferential of the geodesic with respect to the metric given on the second figure. are embedded in the L q * space. To do this, a vector ξ ∈ R N , with N = h −2 , is identified to the function ξ ∈ L q * (Ω) that takes the constant value ξ i,j on the axis-aligned square of width h centered around x i,j . This actually means that the corresponding function are defined as piecewise constant on the dual lattice of the discretization grid.
The discrete optimization problem we have on R N corresponds to solving where the space L q * N (Ω) is the space of L q * functions which are constant on every square of the N −th lattice. The measure γ N is a discretization of γ on points of the lattice, and satisfies γ N γ. The functionc N ξ is the piecewise affine extension of c h ξ on a triangular grid refining the square lattice (each square being divided into two triangles). One can notice that in order to compute integration, any extension with the same values on the support of γ N would have given the same result. The same energy J N may be extended to the whole L q * (Ω) by +∞ outside L q * N (Ω). Before presenting the Γ−convergence result, we need to add some notions and equivalences concerning the meaning of c ξ in the continuous case, when ξ ∈ L q * . Section 2.2 introduces the quantity c ξ by approximation with continuous functions. We are here concerned with the functions T → c ξ (S 0 , T ) that is denoted U ξ for simplicity when the dependence on S 0 is not ambiguous. Let us consider also the function v ξ given by which is the maximal a.e. subsolution of |∇v| = ξ. What we are looking for is very much linked to viscosity solution of the Hamilton-Jacobi equation |∇v| = ξ, when ξ is only L q * (see [6] for a definition via W 1,q test functions and local a.e. inequality and [7] for some other notions and equivalences in the case ξ ∈ L ∞ ). Yet, we are here interested only in the following two notions: the one presented in [8] (and in Section 2.2) and the maximal a.e. one.
We can now go through the essential part of the Γ−convergence proof. We denote byŨ N ξ N the functions T →c N ξ N (S 0 , T ).
. Then for each S 0 the sequence (Ũ N ξ N ) N is bounded in W 1,q * and any weak limit u satisifes u ≤ v ξ . Proof. Let us prove boundedness. We consider at first derivatives with respect to the component x of the variable: sinceŨ N ξ N is an affine extension, its ∂/∂x derivative on a triangular cell coincides with the same derivative on the horizontal side of the cell. The slope on the side is bounded by the value of ξ N at one of the vertices. If T is a triangle of the partition (its area being A), P 1 and P 2 are the two vertices we consider and Q 1 and Q 2 the two squares where ξ N is constant (with area 2A each), respectively centered around P 1 and P 2 , we get Summing up over all triangles any square is counted at most four times, and in the end we get This proves boundedness (since (ξ N ) N is bounded in L q * and the same argument may be performed for ∂/∂y) and the existence of weak limits u (which are uniform limits as well, and hence the condition u(S 0 ) = 0 is conserved).
Let us now define Since ξ N is bounded in L q * , the same is true for ξ N,x and ξ N,y and we denote by ξ x and ξ y their weak limits (up to subsequences). Let us fix an arbitrary rectangle R whose vertices belong to the N −th lattice. Let us now estimate the positive part of ∂Ũ N ξ N /∂x on R: in this case it is sufficient to estimate all ascending slopes. Take an horizontal line in R: all these slopes coincide with the slopes on the horizontal boundaries of the triangles and are hence bounded by the values of ξ N,x at the right extremal points of the segments bounding the triangles. Fixing the ascending direction allows us to avoid superpositions in the choice of the extremal point. The same can be done for descending slopes and for vertical slopes. In the case of ascending horizontal slopes we get where R N is a rectangle a little bit larger than R, since it includes all the squares in the dual lattice around the points of ∂R. Yet, |R N \ R| → 0 as N → ∞. When N → ∞ we can pass to the limit and obtain R ∂u ∂x which, R being arbitrary, implies ∂u ∂x and, taking the maximum over the positive and the negative part, This proves |∇u| ≤ ξ 2 x + ξ 2 y . It is sufficient to prove ξ 2 x + ξ 2 y ≤ ξ to obtain that u is an a.e. subsolution of the Hamilton-Jacobi equation, and hence u ≤ v ξ . The inequality we want is a consequence of standard properties of weak convergence: if a sequence of vectorvalued maps z N (in this case, the pairs (ξ N,x , ξ N,y )) weakly converges to a map z (in this case z = (ξ x , ξ y )), then, for any convex function f (in this case f (x, y) = x 2 + y 2 ), f (z N ) h implies f (z) ≤ h (to prove it, it is sufficient to express f as a supremum of affine functions). The thesis in hence proven.
Theorem 5.5. The sequence of functionals J N previously defined Γ−converges with respect to the weak L q * convergence to the limit functional J. Moreover, as the sequence (J N ) N is equi-coercive and every functional, J included, is strictly convex, convergence of the unique minimizers and of the values of the minima is guaranteed.
Proof. To prove Γ−convergence it is sufficient to prove that ξ N ξ implies and that for every ξ there exists a sequence ξ N ξ such that J N (ξ N ) → J(ξ). For proving the first part (Γ−liminf inequality), it is sufficient to notice that the first addend in the functionals is lower semicontinuous w.r.t. weak convergence, and that by Lemma 5.4, the second behaves the same way (but we have to use the equality provided by Lemma 5.3).
For the second part (Γ−limsup inequality), we know by [10] that it is sufficient to prove it for ξ in a class which is sufficiently dense ("dense in energy", i.e. so that we can guarantee approximation with convergence of the limit energy J). The class of regular functions satisfies this density property. Indeed, by convolution we know that both terms of the functional J converge (for the first it is dominated convergence and for the second use Lemma 5.3 again). For regular functions it is a well known fact that we can choose ξ N by taking the values of ξ at the points of the lattice and we have uniform convergence of the solutions computed by the Fast Marching Method to the true solution of the Eikonal equation (this fact is mainly proved in [11], with techniques coming from [2]; for other convergence results on monotone discrete schemes for Hamilton-Jacobi equations see also [9] and [14]). 6. Numerical results. In this section some numerical results are shown. Algorithm described in section 5.2 is used. The equilibrium states are visualized using the metric ξ computed by this algorithm.
The first example consists in a source-target pair (S, T ) and a homogeneous functional H, which does not depend on the spatial location. As previously seen, ξ = H (i), and H is strictly increasing function, thus (at least in the case where H does not depend on x), ξ and i are visualized as very similar. In the first example, we took H * (x, ξ) = 1 3 ξ 3 , so that H(x, i) = 2 3 i 3 2 . On Figure 2 the equilibrium metric on a "desert" square is shown. Note that the equilibrium metric is symmetric as the initial configuration of the source and the target.
Next, we consider the case where the functional H is spatially varying. We take H(x, i) = 2 3 ξ 0 (x)i Points x i for i ∈ {1, 2, 3, 4} are chosen in a symmetric manner (see figure 3 (a)). We took A = 10, B = 4, σ 1 = 15 and σ 2 = 7 on a square grid of size 100 × 100. For figures 2 and 3, we used a single source-target pair, thus the traffic weight is γ = 1.
Let us now introduce obstacles and multiple sources and targets. In a symmetric configuration of two sources S 1 and S 2 , and two targets T 1 and T 2 ; we consider a river where there is no traffic and a bridge linking the two sides of the river (see figure 4 (a)). We chose the traffic weights such that γ 1,1 + γ 1,2 = 2(γ 2,1 + γ 2,2 ) and  Figure 2. Traffic congestion equilibrium metric on a desert square between two points. Top: the equilibrium metric ξ obtained via the subgradient descent is shown as a 3D surface. Bottom, left: flat view of the same metric. Bottom, right: the minimal action map U associated to ξ, together with some isolevels (black curve) and some geodesics (white curves). On the last image, one can see that the Wardrop conditions are satisfied. γ2,2 γ2,1 = γ1,1 γ1,2 = 2. The traffic intensity going out from S 1 is twice S 2 's. One can note the two hollows on each side of the river appearing because of the inter-sides and intra-sides crossed traffics (see figure 4).
Finally, we solve with our algorithm an asymmetric and general problem with several sources and targets. In figure 5 (a)