Florent Berthelin
email: florent.berthelin@unice.fr

Thierry Goudon
email: thierry.goudon@inria.fr

Bastien Polizzi
email: bastien.polizzi@imft.fr

Magali Ribot
email: magali.ribot@univ-orleans.fr

Asymptotic problems and numerical schemes for traffic flows

Keywords: 2010 Mathematics subject classification. 65M08, 35L60, 90B20, 35Q91

come

Introduction

In order to describe traffic flows and to reproduce the formation of congestions, several models based either on Ordinary Differential Equations (ODE) or Partial Differential Equations (PDE) have been proposed. Starting from individual-based "Follow-the-Leader" models [START_REF] Gazis | Nonlinear follow-theleader models of traffic flow[END_REF], a very active stream in the traffic community considers now PDE models. A first example dates back to Lighthill and Whitham in the 50's [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF]: the evolution of the density of cars is described by means of a mass conservation equation, where the flux is defined by a prescribed function of the density. In these so-called first-order models, the relation between flux and density is referred to as the fundamental diagram in the traffic flows community. A more accurate description can be expected by considering second-order models where a system of PDE governs the evolution of the density and the speed of cars. A first attempt in this direction is due to Payne [START_REF] Payne | Freflo: A macroscopic simulation model of freeway traffic[END_REF], strongly inspired by the principles of fluid mechanics. However, Daganzo [START_REF] Daganzo | Requiem for second-order fluid approximations of traffic flow[END_REF] pointed out the drawbacks of this approach: the Payne-Whitham model may lead to inconsistent behaviors for the flow, such as vehicles going backwards. The model introduced independently by Aw and Rascle [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] and by Zhang [START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF], which still has the form of a 2 × 2 system of conservation laws, is intended to correct these inconsistencies. In [START_REF] Aw | Derivation of continuum traffic flow models from microscopic follow-the-leader models[END_REF], a derivation of the system is proposed from a Follow-the-Leader model. We can also mention that some kinetic models [START_REF] Bellomo | On the multiscale modeling of vehicular traffic: From kinetic to hydrodynamics[END_REF][START_REF] Nelson | The Prigogine-Herman kinetic model predicts widely scattered traffic flow data at high concentrations[END_REF][START_REF] Paveri-Fontana | On Boltzmann-like treatments for traffic flow: A critical review of the basic model and an alternative proposal for dilute traffic analysis[END_REF][START_REF] Puppo | Fundamental diagrams in traffic flow: the case of heterogeneous kinetic models[END_REF][START_REF] Wegener | A kinetic model for vehicular traffic derived from a stochastic microscopic model[END_REF] are under consideration, after the pioneering work [START_REF] Prigogine | Kinetic Theory of Vehicular Traffic[END_REF].

Further details and references can be found in the survey [START_REF] Bellomo | On the modeling of traffic and crowds:a survey of models, speculations, and perspectives[END_REF]. This work is concerned with the numerical simulation of certain variants of the Aw-Rascle-Zhang model. Let ρ(x, t) and v(x, t) be the density and the velocity of cars at position x ∈ R and time t > 0, respectively. The

   ∂ t ρ + ∂ x (ρv) = 0, ∂ t v + p(ρ) + v∂ x v + p(ρ) = 0, (1)
where ρ → p(ρ) plays the role of the pressure in the gas dynamic equations.

In fact, the quantity w = v+p(ρ) describes the desired velocity of the drivers, whereas v corresponds to the actual velocity of the cars. Therefore, the term p(ρ), the velocity offset, stands for the difference between these two velocities, reflecting the fact that the drivers slow down because of the density of cars.

It is convenient to rewrite the equations in the more convenient form of a conservative system; namely (1) is, at least formally, equivalent to

   ∂ t ρ + ∂ x (ρv) = 0, ∂ t ρ(v + p(ρ)) + ∂ x ρv(v + p(ρ)) = 0. (2)
Of course, a crucial modeling issue relies on the expression of the velocity offset p(ρ). At first glance, again inspired from gas dynamics, we can set p(ρ) = ρ γ for some γ > 1. However, such a model does not permit to impose a priori a limitation to the cars density. Consequently, Berthelin, Degond, Delitala and Rascle proposed in [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF] to define the velocity offset as follows:

given 0 < ρ < ∞, ρ ∈ [0, ρ ) -→ p(ρ) = ρ ρ ρ -ρ γ , γ > 1,
which can be rewritten, for ρ = 0, by p(ρ

) = 1 ρ - 1 ρ -γ
, where ρ denotes a maximal value for the density. The velocity offset tends to infinity when ρ → ρ while we get the classical expression p(ρ) ∼ ρ γ when ρ → 0. Such a pressure law also arises in gas dynamics, where it is referred to as the Bethe-

Weyl law [START_REF] Benzoni-Gavage | Multidimensional Hyperbolic Partial Differential Equations: First-Order Systems and Applications[END_REF]; for instance it is used to model close-packing effects in multifluid flows, see [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF] and the references therein. This expression for p has the role of enforcing that ρ satisfies the constraint 0 ≤ ρ(x, t) ≤ ρ for all x ∈ R and t > 0, as it can be seen from the bounds on the Riemann invariants of the system [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF]. Moreover, [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF] points out that drivers do not reduce significantly their speed unless they reach a congested region. Accordingly the velocity offset is appropriately rescaled with a small parameter 0 < ε 1, and we use

p ε (ρ) = ε ρ ρ ρ -ρ γ , for 0 ≤ ρ < ρ , (VO1)
in the Aw-Rascle-Zhang model. We are thus led to the Rescaled Modified Aw-Rascle (RMAR) system

       ∂ t ρ ε + ∂ x (ρ ε v ε ) = 0, ∂ t ρ ε (v ε + p ε (ρ ε )) + ∂ x ρ ε v ε (v ε + p ε (ρ ε )) = 0. (3)
In this model, the velocity offset is small unless the density is getting close to the threshold ρ . Finally, [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF] studies the limit when ε → 0 in [START_REF] Bellomo | On the multiscale modeling of vehicular traffic: From kinetic to hydrodynamics[END_REF]. In this regime we obtain (at least formally) the constrained system

           ∂ t ρ + ∂ x (ρv) = 0, ∂ t (ρ(v + π)) + ∂ x (ρv(v + π)) = 0, 0 ≤ ρ ≤ ρ , π ≥ 0, (ρ -ρ)π = 0. (4)
In ( 4), the limit "pressure" π = lim ε→0 p ε (ρ ε ) appears as the Lagrange multiplier associated to the unilateral constraint 0 ≤ ρ ≤ ρ . In particular, π becomes active only in the congested regions, where the density reaches the threshold ρ . Otherwise, in absence of congestion, the system is reduced to the pressure-less gas dynamic model [START_REF] Brenier | Sticky particles and scalar conservation laws[END_REF][START_REF] Grenier | Existence globale pour le système des gaz sans pression[END_REF]   

∂ t ρ + ∂ x (ρv) = 0, ∂ t (ρv) + ∂ x (ρv 2 ) = 0. (5)
The asymptotic model is further investigated in [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF], exhibiting the formation of clusters, and proving the existence of weak solutions to the system (4) through the stability analysis of "sticky blocks" dynamics. It is also worth pointing out the original numerical approach developed in [START_REF] Maury | Pressureless Euler equations with maximal density constraint : a time-splitting scheme[END_REF] for (4) which uses ideas from the modelling of crowd motion and includes a fine description of the non elastic collision processes.

The asymptotic system ( 4) is thus specifically intended to describe the formation and the dynamics of jams. In this paper, we are interested in the numerical simulations of the system (4), and in the asymptotic regime ε → 0 in [START_REF] Bellomo | On the multiscale modeling of vehicular traffic: From kinetic to hydrodynamics[END_REF]. The difficulty is two-fold. On the one hand, in the free flow case, it is well-known that the pressureless gas dynamics system (5) can lead to deltashocks formation, which makes it difficult to treat numerically [START_REF] Bouchut | Numerical approximations of pressureless and isothermal gas dynamics[END_REF]. On the other hand, with the formation of a congestion, there is no direct access to the limit velocity offset π which is defined in a quite abstract way. Therefore, in order to go beyond the simple particulate approach in [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF], we wish to develop numerical simulations of the RMAR model (3) with the velocity offset (VO1) for small values of ε. We are still facing several numerical challenges.

Firstly, the model prohibits that the density exceeds the threshold ρ . Secondly, (one of) the characteristic speeds of the system become very large in congested region, which makes the time step shrink: the smaller ε, the more severe the stability constraint. Therefore, we need to design a scheme which can preserve the natural estimates of the problem, in particular the density limitation. As already observed in [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF][START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF] standard schemes may fail this objective due to the very specific structure of the PDE system. We also refer the reader to [START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF] for further examples related to fluid mixtures. Moreover, we would like to relax the stability constraints on the time step. To this end, a first attempt would be to adapt the explicit-implicit method pro-posed in [START_REF] Degond | Numerical simulations of the Euler system with congestion constraint[END_REF] for the Euler system with congestion constraint. However, this method applied to (3) is not satisfactory: it produces excessively smoothed density profiles and it overestimates the velocity in congested regions. Thus we design a different splitting strategy, partly inspired from [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF]. Beyond the conception of a numerical scheme able to handle the stiffness of the problem, we will also discuss different asymptotic approaches of the constrained problem (4), based on different definitions of the scaled offset velocity in [START_REF] Bellomo | On the multiscale modeling of vehicular traffic: From kinetic to hydrodynamics[END_REF] which all lead asymptotically to (4). It is interesting to study how the shape of the pseudo-pressure affects the intermediate states (for not so extreme values of the scaling parameter), and the numerical costs.

The outline of this article is the following. In Section 2, we go back to some properties of the Aw-Rascle-Zhang system and we detail the numerical difficulties we face. Additionally, we propose different velocity offsets and scaling that can be used to recover asymptotically the constrained system (4). Then, in Section 3, we propose a new explicit-implicit scheme based on a splitting strategy. The splitting is constructed to reduce the characteristic speeds in the explicit part so that we can expect to use larger time steps.

Finally, in Section 4, we display some numerical simulations in order to prove the efficiency of the scheme and to compare the behavior of the system when using different velocity offsets.

2 Properties of the Aw-Rascle-Zhang model and numerical difficulties

We will describe in this Section the main numerical difficulties we have to deal with, when computing solutions of system (3).

Different velocity offsets

With the velocity offset (VO1), it is forbidden to produce numerical densities larger than the threshold ρ : if, due to any numerical error, the code returns a density larger than ρ , we cannot update further the system and the simulation breakdowns. To cope with this difficulty, we propose to slightly modify the law, replacing (VO1) by a function ρ → pε (ρ) which is defined for any positive entry, which behaves like p ε for ρ < ρ , and which blows up as ε → 0 for ρ ≥ ρ . For instance, we set

pε (ρ) =      ε ρ ρ ρ -ρ γ , if ρ ≤ ρ ε tr , c 0 ε + c 1 ε (ρ -ρ ε tr ) + c 2 ε (ρ -ρ ε tr ) 2 2 , if ρ > ρ ε tr . (VO2)
In this formula, ρ ε tr is a transition density, which has a modeling nature; it should satisfy ρ ε tr → ρ as ε tends to 0. Beyond the transition, pε is a second order polynomial, computed so that pε remains a C 2 function. We thus set

ρ ε tr = ρ -h(ε), c 0 ε = p ε (ρ ε tr ) = ε ρ (ρ -h(ε)) h(ε) γ , c 1 ε = (p ε ) (ρ ε tr ), c 2 ε = (p ε ) (ρ ε tr ).
The expected behavior holds for instance with h(ε) = ε, since it satisfies the two following properties when ε → 0: h(ε) → 0 and c 0 ε → +∞. We point out that ρ ε tr is purely a modeling parameter and the model (VO2) leads us to the same difficulties as (VO1) in terms of stability issues, as we will see in the next Section.

Another option is to use the following velocity offset

p γ (ρ) = V ref ρ ρ γ , γ > 1, (VO3)
for large values of the exponent γ. In this formula V ref > 0 is a reference velocity, to bear in mind the physical meaning of p (in the numerical simulations below we will simply set V ref = 1). This approach is used in fluid mechanics, for modeling certain free boundary problems where bubbles are immersed in a gas [START_REF] Lions | On a free boundary barotropic model[END_REF]. This function is defined on [0, ∞), it behaves proportionally to the gas-law ρ γ for ρ → 0 and it blows up as γ → +∞ for ρ ≥ ρ . Using (VO3) and the regime of large γ's in traffic flows modeling is quite new; we shall see that this simple law has certain advantages in the numerical simulations of congested situations.

In what follows, p refers to (VO1), (VO2) or (VO3). We will see that similar behaviors, corresponding to what can be expected for (4), are captured asymptotically (namely as ε → 0 or γ → ∞) by these velocity offsets. However, the intermediate behaviors can significantly differ and the definition of the velocity offset seriously impacts the numerical costs.

Stability issues

As long as the functions ρ and v are smooth enough, we can rewrite system (1) in the fully non-conservative form

∂ t   ρ v   + A(ρ, v)∂ x   ρ v   = 0, A(ρ, v) =   v ρ 0 v -ρp (ρ)   .
The two eigenvalues related to the system are therefore equal to

λ 1 = v -ρp (ρ) ≤ λ 2 = v (6)
with related eigenvectors

r 1 =   1 -p (ρ)   , r 2 =   1 0   .
The system is strictly hyperbolic, away from the regions where ρ = 0. Let us just note that the largest eigenvalue λ 2 is always linearly degenerate, leading to contact discontinuities and that λ 1 is genuinely non-linear, except for certain forms of the velocity offset we are not considering here. Therefore, the first eigenvalue will admit shocks or rarefaction waves. One of the difficulties of the computations is that vacuum regions may appear. Observe that the information does not travel faster than the actual cars speed v, and that the system preserves the natural properties ρ ≥ 0, v ≥ 0.

We are considering here some Finite Volume (FV) numerical schemes in order to compute the solutions of system (1). Let us denote by ∆t and ∆x the time step and the space step of the method, respectively. We consider the discrete times t n = n∆t, for n ∈ N and the discretization cells

C j = [x j-1/2 , x j+1/2
], j ∈ Z (neglecting for the time being the issue of the boundary conditions) where x j+1/2 = (j + 1/2)∆x. We go back to the conservative form (2) and we denote

U (x, t) =   ρ(x, t) y(x, t)   with y(x, t) = ρ(v + p(ρ))(x, t),
the conservative variables. In terms of the conservative variables ρ and y,

we simply have    ∂ t ρ + ∂ x (ρv) = 0, ∂ t y + ∂ x (yv) = 0,
which recasts as follows, using only the variables ρ and y,

     ∂ t ρ + ∂ x (y -ρp(ρ)) = 0, ∂ t y + ∂ x y 2 ρ -yp(ρ) = 0. ( 7
)
The numerical unknown U n j = (ρ n j , y n j ) is thought of as an approximation of the mean value of U (x, t) on the cell C j at time t n . The FV scheme has the following general form

U n+1 j = U n j - ∆t ∆x F n j+1/2 -F n j-1/2
which mimics what we obtain by integrating the continuous equation [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF] over t n , t n+1 × C j . For the simple schemes we wish to deal with, the numerical flux at the interface x j+1/2 is a function of the neighbouring cells

F n j+1/2 = F U n j+1 , U n j .
Without entering into the details of the schemes, the numerical stability of such a method relies (at least) on the following inequality, see for example [9, Sect. 2.3.3] or [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF],

∆t ≤ 1 2 ∆x max(|λ 1 |, |λ 2 |) , (8)
where λ 1 and λ 2 are the two eigenvalues given by ( 6). It imposes that the state U x j+1/2 , t n+1 on the interface x j+1/2 at time t n+1 = t n + ∆t only depends on the states of the unknown at time t n on the neighbouring cells C j and C j+1 . It can thus be obtained by solving the corresponding Riemann problem with data U n j and U n j+1 .

Let us first consider the case of the pressure (VO1). In case of a congestion formation, ρ ε → ρ but we expect that p ε (ρ ε ) remains bounded and admits the limit π as ε → 0; it leads to the ansatz

ρ = ρ -O ε 1/γ , when ε → 0.
Accordingly the behavior of the characteristic speeds is given by max λ

(ε) 1 , λ (ε) 2 = O ε -1/γ , when ε → 0,
since v ε should remain bounded when ε → 0. Hence, as ε goes to zero, the time step ∆t shrinks due to the condition [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF]. The same remarks apply to the velocity offset (VO2), which essentially behaves like (VO1) when ε → 0.

For the velocity offset (VO3), we find max λ

(γ) 1 , λ (γ) 2
= O (γ) , when γ → +∞, which again imposes tiny time steps. This observation motivates the design of a scheme based on splitting strategy so that the fast waves can be treated implicitly.

Invariant regions

Let us detail another difficulty which is very specific to the traffic flow system

(1). The Riemann invariants for the system (1) are given by, see [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF],

z 1 = v + p(ρ), z 2 = v.
Therefore, the domain

(z 1 , z 2 ) ∈ R 2 with z 1 ∈ [w m , w M ], z 2 ∈ [v m , v M ]
is an invariant region for (1): if the initial datum lies in such a region, the solution will still be contained in the same region for all times. However, numerical difficulties arise due to the fact that such domains are non-convex for the conserved quantities ρ, y. This point has been observed in [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF] for the traffic flows model, see also [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF][START_REF] Jung | Schémas numériques adaptés aux accélérateurs multicoeurs pour les écoulements bifluides[END_REF] for similar problems. At first sight, it would be tempting to define the numerical fluxes by using the Godunov scheme, which is a standard for systems of conservation laws. It works into two steps. Owing to the stability condition (8), we solve a set of uncoupled Riemann problems, centered at the interfaces x j+1/2 with data U n j , U n j+1 .

Then, we project the obtained piecewise constant solution to obtain the updated numerical unknown, constant on the cells C j . This projection step does not preserve the invariant region, since the latter is non-convex (for the role of the convexity of the invariant domain we refer the reader to [START_REF] Roux | Stability for some equations of gas dynamics[END_REF] for gas dynamics equations, and more generally to the textbook [9, Prop. 2.11]).

A counter-example is detailed in [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF] to explain why the Godunov scheme fails to satisfy the maximum principle for the Riemann invariants of (1), and especially for v. To solve this difficulty, [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF] proposes a hybrid scheme which mixes the Glimm scheme to compute the contact discontinuities, and the Godunov procedure to compute 1-shock or 1-rarefaction wave. This hybrid method is well-adapted to handle the specific velocity offsets dealt with in [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF][START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF], see Remark 2 below, which differ from the models we wish to consider here.

We bear in mind that, instead of using the mean of the solutions of the Riemann problems over the cells, the Glimm scheme uses a random sampling strategy in the reconstruction procedure. Hence, by construction, the Glimm scheme preserves the invariant regions, despite the defect of convexity. It is thus well adapted to the simulation of the system (1). Note however the final scheme is non conservative.

Remark 1 Note that, depending on the definition of the numerical fluxes, the stability condition can be even more constrained than [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF], for instance in order to fulfill the bound from above on the density with the "close-packinglike" velocity offset (VO1), see e.g. [START_REF] Berthelin | Multifluid flows: a kinetic approach[END_REF].

Remark 2 In [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF][START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF], the discussion focuses on the velocity offset

p(ρ) = V ref ln ρ ρ , (9)
which has some very specific features:

• First of all, ρ → p(ρ) is defined on (0, ∞) and non decreasing. The model cannot treat vacuum regions since the velocity offset is not defined for ρ = 0. With this model, the velocity offset is well-defined beyond ρ , and we note that p(ρ) < 0 for 0 < ρ < ρ , which might be physically questionable.

• Second of all, we have ρp (ρ) = V ref .

In particular, we get λ 1 = v -V ref so that the CFL condition depends only on v and it does not shrink as the density approaches ρ .

Description of the scheme

Let us now explain in more details the construction of the scheme that we wish to use for the simulation of (1), with a velocity offset ρ → p(ρ) that introduces some stiffness in order to reproduce the expected behavior of the constrained system [START_REF] Bellomo | On the modeling of traffic and crowds:a survey of models, speculations, and perspectives[END_REF]. In what follows, p thus refers to p ε in (VO1), to pε in (VO2) or to p γ in (VO3).

In order to get rid of the large characteristic speeds, the idea consists in splitting the velocity offset into two parts p = p exp + p imp , so that the system with p exp has a stability CFL condition (8) of order 1 with respect to the scaling parameters. Namely, the eigenvalue λ

(ε) exp,1 (resp. λ (γ)
exp,1 ) is bounded with respect to 0 < ε 1 (resp. γ 1). The corresponding system can thus be treated explicitly by means of the Glimm scheme, which preserves the invariant domains. Next, only the stiff part that involves p imp is treated implicitly. We expect also that the implicit part has a simple structure that can be handled with a not too complicated scheme.

Such splitting approach appeared in [START_REF] Degond | Self-organized hydrodynamics with congestion and path formation in crowds[END_REF][START_REF] Degond | Numerical simulations of the Euler system with congestion constraint[END_REF][START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF] for more standard fluid mechanics systems, for instance as an efficient strategy to handle low Mach number regimes. However, as explained above, the structure of the PDE system (1) significantly differs from the Euler equations, in particular with the lack of convexity of the invariant domains and these methods cannot be directly applied to (1).

Definition of the explicit velocity offset

In the first step of the splitting, we consider the system

   ∂ t ρ + ∂ x (ρv) = 0, ∂ t ρ(v + p exp (ρ)) + ∂ x ρv(v + p exp (ρ)) = 0. ( 10
)
It has the same structure as the Aw-Rascle-Zhang system (2), just replacing the full velocity offset p, that can be (VO1), (VO2) or (VO3), by p exp . The characteristic speeds are

λ 1 = v -ρp exp (ρ), λ 2 = v.
In order to relax the stability constraint, we define p exp so that the characteristic speed does not blow up as the scaling parameters ε goes to 0 or γ goes to ∞. It leads to require that p exp (ρ) is bounded uniformly with respect to ε (resp. γ for the law (VO3)), when ρ lies in a compact set of (0, ∞).

The definition depends on a truncation parameter 0 < ρ num < ρ . Following an idea of [START_REF] Degond | Self-organized hydrodynamics with congestion and path formation in crowds[END_REF][START_REF] Degond | Numerical simulations of the Euler system with congestion constraint[END_REF], we set

p exp (ρ) =            p(ρ), if 0 ≤ ρ ≤ ρ num , p(ρ num ) + p (ρ num )(ρ -ρ num ) + p (ρ num ) (ρ -ρ num ) 2 2 , if ρ > ρ num . (11)
We use a second order polynomial for ρ > ρ num to ensure that p exp is a C 2 function. We will see below that it is important to reach such a regularity.

The transition density ρ num is chosen so that p exp (ρ num ) is bounded when

ε → 0 (resp. γ → +∞).
If such a condition is satisfied, then the two eigenvalues are bounded with respect to ε (resp. γ) and therefore the stability condition ( 8) is independent of ε (resp. γ). Let us now explain how to choose ρ num for the velocity offsets under consideration.

a) For the laws (VO1) and (VO2), when ε → 0 and ρ ε → ρ , we expect

that ρ ε = ρ -O(ε 1/γ ), since p ε (ρ ε ) remains finite. We set ρ num = ρ (1 -δρ) with δρ > 0. A simple computation shows that p exp (ρ ) = p (ρ num ) + p (ρ num )(ρ -ρ num ) = O ε(δρ) -(γ+1) ,
which leads us to set

δρ = ε 1/(γ+1) and ρ num = ρ 1 -ε 1/(γ+1) .
For the velocity offset (VO2), we point out that ρ num is a truncation parameter of numerical nature, designed to ensure a certain stability property of the scheme, while ρ trans relies on modeling consideration and does not prevent the blow up of the velocity offset. In any case, we have 0 < ρ num < ρ trans < ρ .

b) For the law (VO3), we require that p exp (ρ num ) remains bounded when γ → +∞. Denoting again ρ num = ρ (1 -δρ), a simple computation leads to

p exp (ρ num ) = γ ρ (1 -δρ) γ-1 + γ(γ -1) ρ 2 δρ (1 -δρ) γ-2 .
So, we are looking for δρ such that δρ → 0 and both terms in this sum are O(1) when γ → +∞. A simple study of the two sequences of functions (γ exp (-γx)) γ∈N and γ 2 x exp (-γx) γ∈N shows that we should have δρ = O(γ -α ) with α ∈ (0, 1) in order to satisfy the required properties.

In [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equations[END_REF] the pressure term is split as p = p exp + p imp with p exp = αp and

p imp = (1 -α)p.
Here, due to the singularity of p at ρ , this approach does not permit to keep p exp bounded and we have to define p exp and p imp as explained above.

A time-splitting scheme

As said above, it is convenient to work on the conservative form (7) of the system (1), dealing with the unknowns U = ρ, y where y = ρ(v + p(ρ)).

With p = p exp + p imp , we arrive at

       ∂ t ρ + ∂ x (y -ρp exp (ρ) -ρp imp (ρ)) = 0, ∂ t y + ∂ x y 2 ρ -yp exp (ρ) -yp imp (ρ) = 0.
We use now a time-splitting scheme. Knowing some approximate values U n = (ρ n , y n ) at time t n , we proceed as follows.

•

Step 1: Solve with an explicit scheme the system of conservation laws

       ∂ t ρ + ∂ x (y -ρp exp (ρ)) = 0, ∂ t y + ∂ x y 2 ρ -yp exp (ρ) = 0.
As said above, this system has the same structure as the original problem [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF]. In particular the invariant domains are non-convex. It can be solved with the Glimm scheme adapted for the pressure p exp . More details will be given in Section 3.3. It defines some intermediate values (ρ n+1/2 , y n+1/2 ).

• Step 2: Solve implicitly the system

       ∂ t ρ -∂ x (ρ p imp (ρ)) = 0, ∂ t y -∂ x (y p imp (ρ)) = 0. ( 12
)
Note that the system has a simple structure and the two equations decouple. The first equation is a non linear scalar conservation law for the density ρ, the second is a linear transport equation for y where the velocity -p imp (ρ) can be considered as given. The numerical method to solve the system ( 12) is thus not that complicated. More details will be given in Section 3.4.

A few words about the Glimm scheme for (1)

In order to use a Glimm scheme, we need to know the Riemann solutions of the problem. This computation has already been done in [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] and in [7, Section 6] where all the details can be found. We refer the reader to some classical books [START_REF] Smoller | Shock waves and reaction-diffusion equations[END_REF][START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF] for general discussions about the role of Riemann problems in the theory of conservation laws and to [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF][START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF] for the specific case of the traffic flow system. We recap here only the Riemann solutions, omitting the details on the elementary waves and on the admissibility of solutions.

A brief overview on the Riemann problem for (1)

Let us just recall that the second eigenvalue λ 2 given by (6) of system (1) is always linearly degenerate, leading to contact discontinuities and that λ 1 is genuinely non-linear, leading to shocks or rarefaction waves. One of the difficulties of the computations is that vacuum regions may appear. Therefore, the Riemann solution of system (1), with an initial datum

(ρ, v)(x, 0) =    ρ L , v L , for x < 0, ρ R , v R , for x > 0,
can be computed according to the five following cases:

• if ρ L > 0, ρ R > 0 and 0 ≤ v R ≤ v L
, the solution consists of a 1-shock that connects ρ L , v L to the intermediate state (ρ * , v * ) and a contact discontinuity between (ρ * , v * ) and ρ R , v R ;

• if ρ L > 0, ρ R > 0 and v L < v R ≤ v L + p(ρ L
), the solution consists of a 1-rarefaction wave that connects ρ L , v L to (ρ * , v * ) and a contact discontinuity between (ρ * , v * ) and ρ R , v R ;

• if ρ L > 0, ρ R > 0 and v L + p ρ L < v R , a vacuum region appears;
the solution consists of a 1-rarefaction wave that connects ρ L , v L to (0, v * ), then a vacuum region between (0, v * ) and (0, v R ) and a contact discontinuity between 0, v R and ρ R , v R ;

• if ρ L > 0, ρ R = 0, the solution is only a 1-rarefaction wave connecting ρ L , v L and 0, v R ;

• if ρ L = 0, ρ R > 0, the solution is only a 2-contact discontinuity con- necting 0, v L and ρ R , v R .
The intermediate state (ρ * , v * ) of the Riemann solution is computed with the Riemann invariants, that is to say

       v * = v R , ρ * = p -1 v L -v R + p ρ L .
In the case of a shock, the speed of the shock between (ρ * , v * ) and ρ L , v L is given by

s = ρ * v * -ρ L v L ρ * -ρ L .
In the case of a rarefaction wave, the self-similar solution (ρ, v)(ξ) with ξ = x/t is given by the following formulae

       p(ρ(ξ)) + ρ(ξ)p (ρ(ξ)) = p ρ L + v L -ξ, v(ξ) = v L + p ρ L -p(ρ(ξ)), (13)
which apply for ξ ∈ λ 1 ρ L , v L , λ 1 ρ R , v R .

Remark 3

In practice, we compute the self-similar solutions of equation [START_REF] Chalons | Godunov scheme and sampling technique for computing phase transitions in traffic flow modeling[END_REF], by using the Newton algorithm. The method requires that p has the C 2 regularity. This remark explains the construction of the explicit part of the velocity offset (11) (we have observed bad behaviors of the scheme when p has jumps).

Glimm scheme

Hence, we have at hand formula to compute the solution of the Riemann problems, which are the elementary brick of the Glimm's scheme (like for Godunov's scheme). This scheme has been introduced for theoretical purposes [START_REF] Glimm | Solutions in the large for nonlinear hyperbolic systems of equations[END_REF], and its implementation for hyperbolic systems is further discussed in [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF][START_REF] Colella | Glimm's method for gas dynamics[END_REF][START_REF] Liu | The deterministic version of the Glimm scheme[END_REF]. Let U n j = ρ n j , y n j be the approximated mean value on the cell C j of U = (ρ, y) at time t n . We proceed as follows to update the numerical unknown:

• We solve the associated Riemann problem at each interface x j+1/2 , namely all the Riemann problems with

U L = U n j , U R = U n j+1 .
• Let a n be a number picked randomly in [0, 1]. We define the value U n+1 j of the numerical unknown in the cell C j at time t n+1 as to be the solution of the Riemann problem evaluated at the point x j-1 2 +a n ∆x ∈ C j . The scheme does not use any averaging or projection procedure and the obtained solution, by construction, remains in the invariant region of the PDE system. In practice, we use the Van Der Corput quasirandom sequence (a n ) n∈N (see [START_REF] Colella | Glimm's method for gas dynamics[END_REF]) defined by

a n = m k=0 i k 2 -(k+1)
where n = m k=0 i k 2 k , with i k ∈ {0, 1}, denotes the binary expansion of the integer n.

Treatment of the implicit part

Let us now discuss how we handle the system [START_REF] Chalons | Transport-equilibrium schemes for computing contact discontinuities in traffic flow modeling[END_REF] where we remind the reader that p imp contains the stiff part of the velocity offset. As said above, the system decouples and it has a very simple structure. Let us set Φ :

ρ ∈ [0, ∞) → -ρp imp (ρ) ∈ (-∞, 0]
. We solve the scalar conservation law for ρ with the classical Engquist-Osher scheme. The Engquist-Osher flux is defined by

F (ρ j+1 , ρ j ) = R [Φ ] + (ξ)1 0≤ξ≤ρ j dξ + R [Φ ] -(ξ)1 0≤ξ≤ρ j+1 dξ,
where [X] + = max(X, 0) and [X] -= min(X, 0). The implicit scheme takes the following form

ρ n+1 j = ρ n+1/2 j - ∆t ∆x F ρ n+1 j+1 , ρ n+1 j -F ρ n+1 j , ρ n+1 j-1 , (14)
where ρ n+1/2 j is the result from the first (explicit) step of the scheme.

Here, the velocity field Φ (ρ) = -d dρ ρp imp (ρ) = -p imp (ρ) -ρp imp (ρ) of this scalar equation is always non positive. Accordingly, the numerical flux reduces to a mere function of the right density

F(ρ j+1 , ρ j ) = R Φ (ξ)1 0≤ξ≤ρ j+1 dξ = Φ(ρ j+1 ).
Consequently, the non-linear equation ( 14) becomes

ρ n+1 j = ρ n+1/2 j - ∆t ∆x Φ ρ n+1 j+1 -Φ ρ n+1 j .
It forms a triangular system of non linear scalar equations. If J stands for the number of grid points, we have to solve J non linear scalar equations, which means J executions of the scalar Newton algorithm to update the density. Then a mere linear system defines the updated y. Indeed, once ρ n+1 is determined, we solve the transport equation for updating y. To this 20 end, we use the standard implicit upwind scheme

y n+1 j = y n+1/2 j - ∆t ∆x G n+1 j+1/2 -G n+1 j-1/2 with G n+1 j+1/2 = y n+1 j -p imp ρ n+1 j + + y n+1 j+1 -p imp ρ n+1 j+1 - .
The specific case p imp ≥ 0, yields

y n+1 j = y n+1/2 j - ∆t ∆x -p imp ρ n+1 j+1 y n+1 j+1 + p imp ρ n+1 j y n+1 j .
It forms a triangular linear system of equations that can be solved by backward substitution, leading to the straightforward formula

y n+1 j = y n+1/2 j + ∆t ∆x p imp ρ n+1 j+1 y n+1 j+1 1 + ∆t ∆x p imp ρ n+1 j .

Simulations

As indicated in the introduction, it is far from clear how to design a "natural" scheme for a direct simulation of the constrained model ( 4). We can only mention the recent approach for crowd dynamics proposed in [START_REF] Maury | Pressureless Euler equations with maximal density constraint : a time-splitting scheme[END_REF]; here we are rather motivated by the asymptotic issues. Our aim is two-fold. On the one hand we wish to discuss the asymptotic behavior of the different models (VO1), (VO2) and (VO3) for the velocity offset, which are all expected to capture asymptotically (for ε → 0 or γ → ∞) the features of the limit system (4). On the other hand, we shall discuss the numerical difficulties and the ability of the time-splitting strategy, which will be compared to the standard Glimm scheme with a small enough time step, in handling the asymptotic behaviour.

The simulations presented below are thought of as Riemann problems and we impose boundary conditions that maintain constant the inflow conditions.

Of course, the method can be adapted to treat further boundary conditions.

In particular imposing zero-influx produces vacuum regions, a numerical difficulty that our method is able to handle, as shown with the decongestion case below.

Case of simple transport

To begin with, we test the case of a simple transport: the computational domain is the interval [0, 1] and for the initial data we set

v 0 (x) = 1, ρ 0 (x) =    0.4, if x ∈ [0, 0.5[, 0.95, if x ∈ [0.5, 1], (15)
see Figure 1 (cyan curves). We compare the six following situations:

• system (1) with pressure (VO1), using the Glimm scheme,

• system (1) with pressure (VO1), using the scheme presented in Section 3,

• system (1) with pressure (VO2), using the Glimm scheme,

• system (1) with pressure (VO2), using the scheme presented in Section 3,

• system (1) with pressure (VO3), using the Glimm scheme,

• system (1) with pressure (VO3), using the scheme presented in Section 3.

The solution at T = 0.4 should be

v(x, T ) = 1, ρ(x, T ) =    0.4, if x ∈ [0, 0.9[, 0.95, if x ∈ [0.9, 1]. ( 16
)
The space step is equal to ∆x = 10 -3 and the time step is computed in order to satisfy the stability condition (evaluated with the full p for the Glimm scheme, and with p exp for the implicit-explicit method). The parameters are taken as follows:

• Pressure (VO1) with γ = 2, ε = 10 -3 , ρ = 1. For the explicit-implicit scheme, the numerical threshold is chosen as

ρ num = ρ 1 -1 5 ε 1/(γ+1) .
• Pressure (VO2) with γ = 2, ε = 10 -3 , ρ = 1 and ρ ε tr = ρ -ε.

For the explicit-implicit scheme, the numerical threshold is chosen as

ρ num = ρ 1 -1 5 ε 1/(γ+1) .
• Pressure (VO3) with γ = 4. For the explicit-implicit scheme, the numerical threshold is chosen as

ρ num = ρ 1 -10 -2 .
The results of the numerical simulations for the three different pressures performed with the Glimm scheme are displayed at Figure 1a-1b, whereas the same simulations using the scheme constructed in Section 3 are exhibited at Figures 1c-1d. All the results are equivalent and agree with the exact solution.

Case of decongestion

Next, we study the case of a decongestion in the traffic. The data are defined by

v 0 (x) =    1, if x ∈ [0, 0.5[, 2, if x ∈ [0.5, 1], ρ 0 (x) = 0.95. ( 17
)
The initial density is close to the threshold. Since the vehicles ahead are going faster, a decongestion occurs. The expected solution at T = 0. (d) Velocity -implicit-explicit scheme be equal to

ρ(x, T ) =            0.95, if x ∈ [0, 0.7[, 0, if x ∈ [0.7, 0.9[, 0.95, if x ∈ [0.9, 1]. (18)
Note the formation of a vacuum region, where v does not make sense. It is indeed particularly interesting and relevant to check the ability of the models and of the numerical methods to handle the formation of vacuum regions, where the density vanishes. The velocity offsets are defined as in the previous Section and we work with the same numerical parameters. The results can be found in Figure 2. The Glimm scheme and the explicitimplicit scheme give the same results for all three pressures (VO1), (VO2) and (VO3). Note that the density and velocity are the same for (VO1) and (VO2), while, for the considered parameters, (VO3) provides significantly different profiles. We observe that pressures (VO1) and (VO2) on the one hand and pressure (VO3) on the other hand give different results, especially in the vacuum region, none of them being totally in agreement with the "expected" result. This can be explained by the moderate value of the parameters ε or γ. Indeed, the constrained behavior (18) can be obtained by changing the parameters, see Figure 3, where we test (VO3) for different values of γ → +∞ and Figure 4, where we make ε → 0 vary for (VO2).

In the former case, the limit behavior is captured with γ = 100 and for the latter, we get a satisfactory result with ε = 10 -5 if γ = 2, and ε = 10 -7 if γ = 3. We notice also at Figure 4 that if we take γ = 3 for pressure (VO2) instead of γ = 2, we need to take a smaller value of ε, namely ε = 10 -7

instead of ε = 10 -5 . Note that the results for the original model (VO1) and for the modified model (VO2) are totally equivalent. The simulations are performed with the implicit-explicit scheme and the initial condition is plotted in cyan.

The two schemes behave equivalently in that case, with some advantage in terms of time step for the implicit-explicit method, see Table 1 (the smaller ε, resp. the larger γ, the more important the gain). Note that in the case when γ = 3, the results are exactly the same because the density is below the numerical threshold; consequently, the implicit part of the scheme is not used and the final result is the same as the one given by the explicit step, that is to say the final result is the same as the one given by the Glimm scheme.

As expected, making ε → 0 for (VO2) or γ → +∞ for (VO3) allows us to obtain a result compatible with (ρ 1 , v 1 ). In particular, we point out that the approach of ( 4) as an asymptotic model from (3) provides the solution where the fastest cars should brake behind the slow vehicles, independently of the density of slow vehicles ahead. This behavior corresponds to the derivation originally introduced in [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF].

A fast cluster reaching a slower one

We now consider the situation of a slow car cluster reached by a faster cluster of vehicles. The initial conditions are

ρ 0 (x) =                0.95 if x ∈ [0.2, 0.3] , 0.9 if x ∈ [0.35, 0.5] , 0 otherwise, v 0 (x) =                2 if x ∈ [0.2, 0.3] , 1 if x ∈ [0.35, 0.5] ,
0 otherwise. and the expected solution at T = 0.3 is (blue), ε = 10 -4 (red), ε = 10 -5 (green) and ε = 10 -6 (black). On top, simulations are performed with the Glimm scheme and on the bottom, with the implicit-explicit scheme. We display the densities on the left and the velocities on the right. The initial condition is plotted in cyan. Figure 13 represents the results for the sub-case (AIII) for the three velocity offsets and for different times. In order to ease the comparison between the Glimm scheme and explicit-implicit scheme presented in Section 3, the results for the two numerical approaches are displayed simultaneously.

ρ(0.3, x) =                1 if x ∈ [0.
According to Figure 13, all the different numerical approaches give similar results which coincide with the solution given by [START_REF] Grenier | Existence globale pour le système des gaz sans pression[END_REF].

Conclusion

The model ( 4) is intended to describe the formation and the dynamics of traffic jams, through a Lagrange multiplier that accounts for a density threshold.

This model can be motivated, at least formally, through asymptotic arguments from the Aw-Rascle-Zhang system with a rescaled velocity-offset. It raises the question of simulating efficiently the Aw-Rascle-Zhang system with potentially stiff velocity offsets. Depending on the values of the parameters it can be seen either as the simulation of a model for traffic flows with stiff parameters or as a way to access the limiting behavior described by (4), alternative for instance to the approach of [START_REF] Maury | Pressureless Euler equations with maximal density constraint : a time-splitting scheme[END_REF]. However the scaling induces fast propagation waves and, in turn, severe stability conditions. In this paper, we propose several approaches to obtain asymptotically (4) and we introduce an implicit-explicit method in order to cope with the large characteristic speeds of the system.

This study exhibits numerical difficulties, related to both the lack of convexity of the invariant domains of (1) and the large characteristic speeds.

We have proposed a time-splitting method, based on a decomposition of the

Figure 1 :

 1 Figure 1: Numerical results in the case of transport (15)-(16). Density (left) and velocity (right) with the Glimm scheme (top) and the explicitimplicit scheme (bottom). The results are given for the three different pressures under consideration: pressure (VO1) in blue, pressure (VO2) in red and pressure (VO3) in green. The initial conditions are plotted in cyan.

  Density -Power-type pressure (VO3).

  Velocity -Power-type pressure (VO3).

Figure 2 :

 2 Figure 2: Numerical results in the case of decongestion (17)-(18).Comparison of the two schemes. Density (left) and velocity (right) at final time T = 0.2. Top: pressure (VO1) with Glimm scheme (blue), implicitexplicit scheme (red) and pressure (VO2) with implicit-explicit scheme (green). Bottom: pressure (VO3) with Glimm scheme (blue) and with the explicit-implicit scheme (red). The initial conditions are plotted in cyan.Parameters : γ = 2, ε = 10 -3 for pressures (VO1) and (VO2) and γ = 4 for pressure (VO3).

Figure 3 :Figure 4 :

 34 Figure 3: Numerical results in the case of decongestion (17)-(18). Pressure (VO3) for different values of γ. Density (left) and velocity (right) at final time T = 0.2: γ = 4 (blue), γ = 20 (red) and γ = 100 (green). The simulations are performed with the implicit-explicit scheme and the initial condition is plotted in cyan.

Figure 5 :

 5 Figure 5: Numerical results in the case of congestion -Comparison of the two schemes. Density (left) and velocity (right) at final time T = 0.01. Glimm scheme (top) and implicit-explicit scheme (bottom), with pressure (VO1) (blue), (VO2) (red) and (VO3) (green). The initial conditions are plotted in cyan. Parameters : γ = 2, ε = 10 -3 for pressures (VO1) and (VO2) and γ = 4 for pressure (VO3).

Figure 6 :

 6 Figure 6: Numerical results in the case of congestion -Comparison of the two schemes -Different parameters. For (VO1) and (VO2), we use γ = 2 and ε = 10 -5 ; for (VO3), we take γ = 50. Pressure (VO1) in blue, pressure (VO2) in red and pressure (VO3) in green. The initial conditions are plotted cyan.

Figure 7 :

 7 Figure 7: Numerical results in the case of congestion -Pressure (VO3) for different values of γ. Density (left) and velocity (right) at final time T = 0.01. Pressure (VO3) for γ = 20 (red), γ = 50 (blue) and γ = 100 (green). The simulations are performed with the Glimm scheme (top) and the implicit-explicit scheme (bottom). The initial condition is plotted in cyan.

Figure 8 : 5 (

 85 Figure 8: Numerical results in the case of congestion -Pressure (VO2) for different values of γ and ε. Density (left) and velocity (right) at final time T = 0.01. We compare the pressure (VO2) for γ = 2, ε = 10 -5(blue), γ = 3, ε = 10 -5 (red) and γ = 3, ε = 10 -7 (green). The simulations are performed with the Glimm scheme (top) and the implicit-explicit scheme (bottom). The initial condition is plotted in cyan.

Figure 9 : 3 (

 93 Figure 9: Numerical results in the case of a shock between two blocks -Pressure (VO1) for γ = 2 and different values of ε: ε = 10 -3 (blue), ε = 10 -4 (red), ε = 10 -5 (green) and ε = 10 -6 (black). Figure 9a (on the left) represents the densities whereas figure 9b (on the right) represents the velocities. The initial condition is plotted in cyan.

Figure 10 :

 10 Figure 10: Numerical results in the case of a shock between two blocks -Pressure (VO2) for γ = 2 and different values of ε: ε = 10 -3

Figure 11 :

 11 Figure 11: Numerical results in the case of a shock between two blocks -Pressure (VO3) for different values of γ: γ = 16 (blue), γ = 32 (red), γ = 64 (green) and γ = 128 (black). On top, simulations are performed with the Glimm scheme and on the bottom, with the implicitexplicit scheme. We display the densities on the left and the velocities on the right. The initial condition is plotted in cyan.

Acknowledgements

We thank Frédéric Coquel for friendly advices and warm encouragements during the preparation of this work.

Case of congestion

A jump of velocity creating a congestion

Finally, we turn to the simulation of a congestion in the traffic. The initial conditions are given by ρ 0 (x) = 0.95,

The density is initially close to the threshold; since the cars ahead are slower, a congestion might occur and the Lagrange multiplier becomes active to prevent an excess of vehicles density.

Indeed, discontinuous solutions are characterized by the Rankine-Hugoniot conditions: with t → s(t) the speed of the discontinuity curve, we have ṡ ρ = ρv , and ṡ ρ(v + π) = ρv(v + π) .

We can check that

and a Lagrange multiplier active only in the congestion domain

is solution of (4). The presence of slow vehicles ahead of the fast ones instantaneously creates a congestion behind the velocity jump: the slow vehicles ahead make the faster ones behind brake. This is typical of the Follow-the-Leader approach, which has led to a derivation of the Aw-Rascle-

Zhang system [START_REF] Aw | Derivation of continuum traffic flow models from microscopic follow-the-leader models[END_REF][START_REF] Gazis | Nonlinear follow-theleader models of traffic flow[END_REF]. However, it is likely that solutions of the constrained model ( 4) are not uniquely defined for such data; we refer the reader to [START_REF] Berthelin | A model for the evolution of traffic jams in multilane[END_REF] for such considerations.

The parameters are defined as in Section 4.1 and we show the solutions obtained at the final time T = 0.01 in Figure 5. We observe exactly the same behaviors between the Glimm scheme and the implicit-explicit scheme with (VO1) or (VO2). We observe that with these parameters, the three models do not find the solution [START_REF] Degond | Numerical simulations of the Euler system with congestion constraint[END_REF]. The time steps for the Glimm scheme are smaller than with the explicit-implicit scheme, but, quite surprisingly, by a factor 3 or 4 only.

Regarding the velocity offsets, (VO3) overshoots the maximal value of the density, equal to 1, whereas the two other pressures (VO1) and (VO2)

underestimate it. This is not surprising since (VO3) allows values larger than the threshold, but it contrasts with the behavior of the model (VO2) which has the same feature. In Figure 6, we make the parameters vary as follows:

• Pressure (VO1) with γ = 2, ε = 10 -5 ,

• Pressure (VO2) with γ = 2, ε = 10 -5 ,

• Pressure (VO3) with γ = 50.

We observe that these parameters provide a result closer to the explicit solution (ρ 1 , v 1 ). The results for (VO3) with different values of γ are given at When using the Glimm scheme, we obtain the expected limit solution as we make the parameters vary: for (VO1) as ε goes to 0, see Figure 9, for (VO2) as ε goes to 0, see Figure 10, and for (VO3) as γ goes to infinity, see Figures 11c and11d. The explicit-implicit scheme is able to compute precisely the car density with the velocity offset (VO3), see Figure 11a.

However, difficulties arise for the computation of the velocity, especially in the back of the jam, see Figure 11b: the velocity of the last cars in the congestion (coordinate x = 0.555) is over estimated. Moreover, we observe in the case when γ = 128 a loss of the quantity of cars between initial and final times of around 9%. This behavior is even worse with (VO2), see Figures 10c and10d. Here again, the quantity of mass is not conserved: the difference between the initial density and the final one is around 23%. We point out that this simulation is quite tough; the results obtained with the Glimm scheme are neater but at the price of a significative numerical cost.

Riemann problems: comparison with [16]

In this section we present the numerical results obtained with the Glimm scheme and the implicit-explicit schemes of Section 3 for the sub-cases AI and AIII of [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF][START_REF] Degond | Modelling and simulation of vehicular traffic jam formation[END_REF]. So we consider in the following a road [0, 1] with the initial data

and

where we use the same labels as in [START_REF] Berthelin | A model for the formation and evolution of traffic jams[END_REF][START_REF] Degond | Modelling and simulation of vehicular traffic jam formation[END_REF]. A solution for (AI) is

and

x ∈ [0.5 + 0.1t, 0.5 + 0.5t] , 0.5 if x 0.5 + 0.5t,

is a solution for (AIII) (note that v does not make sense in the vacuum region x ∈ [0.5 + 0.1t, 0.5 + 0.5t]). In order to compare our scheme with the results of [START_REF] Degond | Modelling and simulation of vehicular traffic jam formation[END_REF] we take the same parameters, namely for the velocity offset (VO1) and (VO2) we use γ = 1 and ε = 10 -3 . For the velocity offset (VO3) we use γ = 64.

Figure 12 represents the results for the sub-case (AI) computed with the three velocity offsets and with the Glimm scheme or the explicit-implicit scheme, at time t = 0.2, t = 0.4 and t = 0.6. According to Figures 12a,12c and 12e, the explicit-implicit scheme for the velocity offset (VO2) overestimates the length of the congestion. For example at t = 0.4 using ( 21), the tail of the congestion should be at x = 0.166 instead of x = 0.13. Moreover,


Khalil Laghmari
email: laghmari.khalil@gmail.com

Christophe Marsala
email: christophe.marsala@lip6.fr

Mohammed Ramdani
email: ramdani@fstm.ac.ma

Graded multi-label classification: compromise between handling label relations and limiting error propagation

In graded multi-label classification (GMLC), each data can be assigned to multiple labels according to a degree of membership on an ordinal scale, and with respect to label relations. For example, in a movie catalog web page, a five stars action movie should be at least a one star suspense movie. Ignoring those relations can lead to inconsistent predictions, but if they are considered, then a prediction error for one label will be propagated to all related labels. Most of existing approaches either ignore label relations, or can learn only relations fitting a predefined imposed structure. This paper is motivated by the lack of a study analysing the compromise between handling label relations and limiting error propagation in GMLC, and by the fact that there is no known approach giving a control on that compromise to allow such a study. In this paper, a new meta-classifier with two main advantages is proposed for GMLC. Firstly, no predefined structure is imposed for learning label relations, and secondly, the meta-classifier is based on three measures giving control on the studied compromise. The studied compromise is analysed according to its impact on the classifier complexity and on hamming-loss evaluation measure. A comparison to three existing approaches shows that the proposed meta-classifier is competitive according to hammingloss evaluation measure, and it is the most stable classifier according to hamming-loss standard deviation.

I. INTRODUCTION

Graded multi-label classification (GMLC) is the task of assigning one or more labels to each data according to an ordinal scale of membership degrees M . It was recently introduced [START_REF] Cheng | Graded multilabel classification: The ordinal case[END_REF] as a generalization of the multi-label classification task (MLC) [START_REF] Tsoumakas | Multi-label classification: An overview[END_REF].

The most known source for GMLC data is catalogue web pages, where movies or animes are assigned to different categories such as action, suspense, and humour using a one-to-five star rating. Data involving the task of GMLC can be found in many other domains such as chemistry, where molecules have multiple odours with different intensities ranging from very weak to very strong.

Relations can exist between labels according to M . For example, high ranked action movies should contain at least little suspense, and molecules having a very strong intensity of jasmine odour can not have the smell of musk simultaneously.

Those relations, if well learned, are supposed to give a better understanding of hidden knowledge in data, but actually they lead to a serious dilemma:

-On the one hand, considering label relations has the advantage of making consistent predictions, but it has the disadvantage of allowing error propagation because a prediction error for one label will be propagated to all related labels. -On the other hand, learning a classifier while ignoring those relations has the advantage of making independent predictions which prevents error propagation, but it has the disadvantage of making inconsistent predictions since label relations are ignored. Label relation dilemma is an inherent challenge in all MLC tasks and not only specific to GMLC. However, this challenge has not received enough attention in the literature: first MLC approaches assume label independence, and later approaches focused only on how to handle label relations.

Learning all label relations can lead to cyclic dependencies between labels. Most of existing MLC approaches avoid this problem by setting a non cyclic dependence structure, and then try to learn only relations fitting the predefined structure.

Another limitation of existing MLC approaches is that they can handle only co-occurrence relations, while GMLC data can encapsulate also order relations based on the ordinal scale of membership degrees.

To overcome the limitations of existing MLC approaches, and to answer the challenge of label relation dilemma in GMLC, we propose a new meta-classifier named PSI-MC with two main advantages: Firstly, it allows label relation learning from GMLC data without imposing a predefined structure, and secondly, based on three measures named preselection, selection, and interest measures, it allows controlling the compromise between learning label relations and limiting error propagation.

The paper is organized as follows: in Section II we review the most used MLC approaches, in Section III we review some well-known GMLC approaches, and in Section IV we describe the proposed meta-classifier PSI-MC. Experiments on real datasets are discussed in Section V. Conclusions and ideas for future works are presented in Section VI.

II. MULTI-LABEL CLASSIFICATION

Let C = {c l } 1≤l≤k be the set of labels, and X = {x i } 1≤i≤n be the set of data. Each data x i is a vector (x ij ) 1≤j≤p where x ij is the value corresponding to the j th attribute. y i ⊆ C is the set of labels associated to x i . The MLC task is to learn a classifier H : X → P(C) mapping each data to the correct set of associated labels. The MLC task is answered either by extending a mono-label classifier, or by transforming multilabel data to mono-label data.

Many mono-label classifiers were adapted to handle multilabel data, such as the K-Nearest Neighbours algorithm [3] [4], the naive bayes algorithm [START_REF] Yan | A Double Weighted Naive Bayes for Multi-label Classification[END_REF], the Support Vector machines [START_REF] Wang | SimplifiedConstraints Rank-SVM for Multi-label Classification[END_REF] [7], the decision trees [START_REF] Madjarov | Hybrid Decision Tree Architecture Utilizing Local SVMs for Multi-Label Classification[END_REF] [START_REF] Wang | Fuzzy Rough Decision Trees for Multi-label Classification[END_REF], and the neural networks [10] [11]. In this paper we are more interested in transformation methods [START_REF] Tsoumakas | Data Mining and Knowledge Discovery Handbook[END_REF] because they can be used with any mono-label classifier.

There are three families of learning approaches based on transformation methods: Learning one multi-classes classifier, learning a binary classifier for each label, and learning a binary classifier for each two different labels.

A. Learning one multi-classes classifier

Label power set (LP) is a straightforward approach considering each label set as a single distinct label. One disadvantage of this approach is the class imbalance problem encountered when some label sets are not frequent. The idea behind pruned problem transformation approach (PPT) [START_REF] Read | A Pruned Problem Transformation Method for Multi-label classification[END_REF] is to remove infrequent label sets in order to overcome the class imbalance problem, but this leads to an information loss due to label truncation. A pruned problem transformation with no information loss (PPT-n) [START_REF] Read | A Pruned Problem Transformation Method for Multi-label classification[END_REF] is an improvement of PPT. Instead of removing infrequent label sets, PPT-n divide them to frequent subsets of labels. A common disadvantage between LP, PPT, and PPT-n approaches is that they can not predict an unseen label set in the training set. This limitation is overcome by the PPT-ext approach [START_REF] Read | A Pruned Problem Transformation Method for Multi-label classification[END_REF].

B. Learning a binary classifier for each label

In binary relevance approach (BR), a binary classifier H l : X → {0, 1} is learned for each label c l ∈ C. The labels relevant to predict are given by H(x) = {c l , H l (x) = 1} 1≤l≤k . This approach learns independent classifiers and therefore can not learn label relations. Classifier chains approach (CC) [START_REF] Read | Classifier chains for multi-label classification[END_REF] overcome this limitation using a chained structure depending on a predefined order for labels. Each classifier H l is allowed to use the prediction result of all previous classifiers {H l } 1≤l <l in order to make its own prediction. In classifier treillis approach (CT) [START_REF] Read | Scalable multi-output label prediction: From classifier chains to classifier trellises[END_REF] a directed graph is built instead of a chained structure. Label correlation is computed between each pair of labels, then classifiers are placed on nodes according to label correlation. Each classifier is allowed to use the prediction result of all classifiers corresponding to parent nodes.

C. Learning a binary classifier for each two different labels

In ranking by pairwise comparison approach (RPC) [START_REF] Hüllermeier | Label ranking by learning pairwise preferences[END_REF], a binary classifier H ll is learned for each label pair c l = c l . Each H ll is used to predict whether a data is associated to the label c l or c l . A label ranking is outputted using a majority vote aggregation but with no separation between relevant and irrelevant labels. This problem is solved by calibrated label ranking approach (CLR) [START_REF] Fürnkranz | Multilabel classification via calibrated label ranking[END_REF]. The idea is to train k more classifiers using an additional label c 0 . For each classifier H l0 , l ∈ [1, k], data not associated with the label c l is considered as associated with the additional label c 0 . Labels predicted more times then c 0 are relevant labels, and the remaining are considered as irrelevant labels.

III. GRADED MULTI-LABEL CLASSIFICATION A. GMLC and fuzzy sets

GMLC data can be viewed as fuzzy sets [START_REF] Zadeh | Fuzzy sets[END_REF]. The only difference is that in fuzzy sets a numeric scale is used for membership degrees M = [0, 1], while in GMLC an ordinal scale is used M = {m 1 < . . . < m s }. The interest of GMLC is related to data acquisition because it is easier for annotators to give their opinions on an ordinal scale. Each data x i is associated with a fuzzy set y i where λ i : C → M is the mapping function of x i associating each label c l ∈ C to its membership grade.

A fuzzy set can be described by the vertical representation using the membership function, or by the horizontal representation using α-cuts. In a similar way, the GMLC task can be solved using either a vertical or an horizontal decomposition, or both [START_REF] Cheng | Graded multilabel classification: The ordinal case[END_REF].

B. Decomposing the GMLC problem

Using the vertical decomposition, k classifiers are trained, one for each label c l ∈ C, so that each classifier H l : X → M predicts the membership grade for the label c l .

Using the horizontal decomposition, s -1 classifiers are trained, one for each grade m g ∈ {m 2 , . . . , m s }, so that each classifier H g : X → P(C) predicts the set of associated labels according to a membership grade at least equals to m g . There is no need to train a classifier for m 1 because the corresponding label set is

H 1 (x) = C.
If a label c l is predicted by a classifier H g , it is expected from all classifiers {H g } 1≤g ≤g to predict c l as well. In regard to this hierarchical property, the trained classifiers are not completely independent, hence the horizontal decomposition in theory, can learn better label relations than the vertical decomposition. However, it is not guaranteed that the hierarchical property is satisfied by trained classifiers. An aggregation function is generally used to handle this case.

Combining both decompositions can be done starting first by the vertical decomposition, then applying an horizontal decomposition, or the opposite. Indeed, the task of each vertical classifier is an ordinal classification [START_REF] Eibe | A simple approach to ordinal classification[END_REF] [START_REF] Cardoso | Learning to classify ordinal data: the data replication method[END_REF] and can be solved horizontally using s -1 binary classifiers, and the task of each horizontal classifier can be solved vertically using the BR approach.

C. Solving the GMLC problem using a pairwise approach

The task of GMLC can be solved also using a pairwise approach like CLR. The three methods named Horizontal CLR, Full CLR, and Joined CLR are all based on the idea of using multiple calibration labels [START_REF] Brinker | Graded multilabel classification by pairwise comparisons[END_REF]. The key idea is to use s -1 virtual labels:

V = {v g } 1≤g≤s-1 with fixed membership grades: ∀i ∈ [1, n] : m g < λ i (v g ) < m g+1 .
v g is used to denote (for simplicity) both the virtual label and its corresponding membership grade:

λ i (v g ) = v g .
For the Horizontal CLR approach, an horizontal decomposition is first performed, then each MLC task g from the obtained s -1 tasks is solved using CLR, with v g as the cutting point between relevant and irrelevant labels.

For the Full CLR approach, a classifier is trained for each label pair in C ∪V . All virtual labels are considered as cutting points. The membership grade m g is predicted for a label c l ∈ C if m g is the highest membership grade for which c l is predicted more times than v g .

One drawback of Full CLR approach, is that for a classifier H ll , a data x i is considered positive if λ i (c l ) > λ i (c l ), and negative otherwise, regardless of the difference between membership grades.

Joined CLR answers this problem by combining both Horizontal CLR and Full CLR approaches. Indeed, each MLC task for the Horizontal CLR approach is solved using all virtual labels as cutting points, instead of using only one cutting point.

The three CLR based approaches discussed in this section are also used in experiment section (V), where our proposed meta-classifier is compared to them.

IV. A NEW META-CLASSIFIER FOR GRADED MULTI-LABEL CLASSIFICATION

A. Key ideas

The first key idea of the proposed meta-classifier is to learn an initial set of k multi-class classifiers H 0 = {H l } 1≤l≤k , one for each label, where the training set for a classifier H l includes membership grades of labels {c l } l =l as descriptive numeric attributes, which means that the prediction by a classifier H l can depend on the prediction result of other classifiers {H l } l =l . This allows learning label relations considering membership grades without fixing a predefined structure, but it does not prevent learning cyclic dependencies.

The second key idea is to avoid cyclic dependencies by replacing involved classifiers. The compromise between handling label relations and limiting error propagation is defined by the way we select which classifier to replace first, and by the way we select alternative classifiers that the new classifier can depend on. The final set of classifiers without cyclic dependencies is the one used to make predictions

H = {H L } 1≤L≤k .
Initially, the final classifier set is empty: H ← ∅, and the set of classifiers not yet added to H is H 0 -H = H. Our proposed meta-classifier PSI-MC removes iteratively classifiers from H, and adds them directly or after being replaced to H, until |H| = 0 and consequently |H| = k.

B. Measures

A pre-selection measure P : H → {0, 1} is used to fill the set of candidate classifiers for replacement: {H l ∈ H, P(H l ) = 1} 1≤l≤k . Note that our objective is not only to solve cyclic dependencies but also to answer the challenge of label relation dilemma. Hence, in order to reduce the effect of error propagation, even if a classifier is not involved in a cyclic dependency, it can be a candidate for replacement if it depends on too many other classifiers.

Classifiers not candidates for replacement are added directly to the final classifier set:

H ← H ∪ {H l ∈ H, P(H l ) = 0} 1≤l≤k .
A selection measure S : P(H) → H is used to select one classifier to be replaced from the set of candidate classifiers. For example, it can be either the classifier that depends on the lowest number of classifiers to reduce information loss, or the classifier that depends on the highest number of classifiers to reduce error propagation. It could also be the classifier whose the highest number of classifiers depend on to resolve more cyclic dependencies.

A measure of chaining interest I : H → {0, 1} is used to decide whether the new classifier to build can depend on a final classifier H L : I(H L ) = 1 or not: I(H L ) = 0.

pre-selection, selection, and chaining interest measures are called the PSI-measures. In the following, we investigate their impact on synthetic data before analysing experiment results on real datasets. In the following, we use decision trees as base classifiers because they are easily interpreted, however the proposed meta-classifier can be used with any other base classifier.

C. Analysing the impact of PSI-measures using synthetic data

In TABLE I, X = {x i } 1≤i≤10 is the training set, C = {c l } 1≤l≤5 is the label set, M = {m g } 1≤g≤4 = {0, 1, 2,
Fig. 1 shows the obtained decision trees for H 0 using the weka implementation [START_REF] Hall | The weka data mining software: An update[END_REF] of the C4.5 algorithm [START_REF] Quinlan | C4.5: Programs for Machine Learning[END_REF], and Fig. 2 shows the corresponding dependency graph. In this example, we choose to make all non independent classifiers candidates for replacement:

P(H l ) = 0 if D ← (H l ) = ∅ , 1 otherwise.
H 1 is independent, hence it is removed from H and added directly to H (Fig. 3). The remaining classifiers are all dependent. One of them should be selected according to a selection measure. In this example we choose to select the classifier allowing us to solve the highest number of cyclic dependencies:

S(H) = argmax H l ∈H (|D → (H l )|). We have |D → (H 3 )| = 1 < |D → (H 2 )| = |D → (H 4 )| = |D → (H 5 )| = 2.
The selection measure outputs the first classifier with the highest value which is H 2 in this case.

H 2 is to be replaced by another classifier H 2 , and we have the choice to chain it with H 1 or to make it independent. In this example the chaining interest measure is given by I(H L ) = 0 , ∀H L ∈ H.

Fig. 4 shows the updated dependency graphs for H and H after replacing H 2 by H 2 . Note that after removing H 2 from H, H 4 becomes independent and then it is added directly to H (Fig. 5). Also note that since H 4 was depending on H 2 in H, it is now depending on H 2 in H. This shows that the selection measure choosing which classifier to replace first, has also an impact on the learned label relations.

→ (H 3 )| = |D → (H 5 )| = 1, then S(H) = H 3 .
H 3 is replaced by an independent classifier H 3 according to the chaining interest measure I (Fig. 6).

After removing H 3 from H, H 5 becomes independent and it is added directly to H (Fig. 7). Since H 5 was depending on both H 2 and H 3 in H, it is now depending on H 2 and H 3 in H.

Fig. 8 shows the dependency graph for H 0 , compared to the one for H using a chaining interest measure always equals to 0 (the same used in the example), and to the one for H using a chaining interest measure always equals to 1. Note 1 is collected from the answers of 1930 graduate students about the importance of 48 properties of their future jobs. Descriptive attributes are age and sex of students, and labels are the 48 properties. To overcome the 1 http://www.ke.tu-darmstadt.de/resources/GMLC problem of insufficient attributes, 50 datasets are randomly generated considering only k labels as target labels and the remaining labels as descriptive attributes. In order to compare our proposed meta-classifier results to reported results of Horizontal CLR, Full CLR, and Joined CLR approaches [START_REF] Brinker | Graded multilabel classification by pairwise comparisons[END_REF], we used the same 50 datasets for k = 5, the same 50 datasets for k = 10, and the same base classifier (weka implementation of C4.5) as in [START_REF] Brinker | Graded multilabel classification by pairwise comparisons[END_REF]. Results averaged over 10 folds crossvalidation are shown in TABLE III, where hamming-loss is extended to the GMLC case as described in [START_REF] Brinker | Graded multilabel classification by pairwise comparisons[END_REF]. III), and almost as good as Joined CLR, while it is the most stable one according to hamming-loss standard deviation.

We generated 15 descriptive attributes for odorous molecule data [START_REF] Arctander | Perfume and Flavor Chemicals: (aroma Chemicals), ser. Perfume and Flavor Chemicals: Aroma Chemicals[END_REF] based on the name and formula of molecules. TABLE IV as TABLE III shows that there is no significant difference between using the measure I = 0 and the measure I = 1 according to hamming-loss. Indeed, by allowing more label relations to be learned (I = 1) more true positive labels are predicted, hence the sensitivity (true positive rate) is increased and hamming-loss is decreased, but due to error propagation, more false positive labels are predicted and consequently hamming-loss is increased. This explains why using I = 1 is not always better.

According to the low sensitivity obtained in TABLE IV, the 15 generated attributes are not relevant to discern labels (molecule odours). The prediction error rate in overlapping regions may be high. This problem can be answered by outputting for each data a set of possible labels for each membership grade [START_REF] Laghmari | A distributed graph based approach for rough classifications considering dominance relations between overlapping classes[END_REF], because experts prefer to have one true prediction in a small set of possible predictions instead of one prediction that could be totally wrong.

TABLE IV THE IMPACT OF CHAINING INTEREST MEASURE IN MOLECULE DATA

Label relation dilemma has also an impact on classifier complexity. TABLE V shows the average number of nodes, leafs, and dependent nodes over all decision trees in each meta-classifier, including the initial meta-classifier with cyclic dependencies H 0 . When descriptive attributes are good enough to discern labels (BelaE-5 and BelaE-10), the classifier complexity (number of nodes and leafs) is reduced by learning more label relations (number of dependent nodes), but the opposite happens when descriptive attributes can not discern labels well enough (molecule data), because more dependent labels will be needed and consequently the complexity is increased. In this paper, we proposed a new meta-classifier for GMLC named PSI-MC. It is based on three measures named PSImeasures for pre-selection, selection and chaining interest measure. PSI-MC has two main advantages: it allows learning label relations without fixing a predefined relation structure, and it allows controlling the compromise between handling label relations and limiting error propagation. Experiment results on real datasets shows that PSI-MC is competitive with other existing approaches according to hamming-loss, and it is the most stable one according to hamming-loss standard deviation.

In this paper, the impact of PSI-measures is not fully analysed, there is many interesting measures to be studied, and there is even more if we use base classifiers outputting attribute weights such as neural networks. For future work, we plan to study further our proposed classifier, by analysing different PSI-measures, with different base classifiers and for different datasets.

Fig. 1 .

 1 Fig. 1. Decision trees for H 0

Fig. 2 .

 2 Fig. 2. Dependency graph for H 0 Let D → : H → P(H) be the function giving for each classifier H l the set of classifiers depending on it, and D ← : H → P(H) be the function giving for each classifier H l the set of classifiers that it depends on.In this example, we choose to make all non independent classifiers candidates for replacement:P(H l ) = 0 if D ← (H l ) = ∅ , 1 otherwise.H 1 is independent, hence it is removed from H and added directly to H (Fig.3).

Fig. 3 .

 3 Fig. 3. Moving H 1 from H to H

Fig. 4 .

 4 Fig. 4. Replacing H 2 by H 2

Fig. 5 .

 5 Fig. 5. Moving H 4 from H to H We have D ← (H 3 ) = D ← (H 5 ) = 1, hence P(H 3 ) = P(H 5 ) = 1. H 3 and H 5 are both candidates for replacement, and since we have |D→ (H 3 )| = |D → (H 5 )| = 1, then S(H) = H 3 .H 3 is replaced by an independent classifier H 3 according to the chaining interest measure I (Fig.6).After removing H 3 from H, H 5 becomes independent and it is added directly to H (Fig.7). Since H 5 was depending on both H 2 and H 3 in H, it is now depending on H 2 and H 3 in H.Fig.8shows the dependency graph for H 0 , compared to the one for H using a chaining interest measure always equals to 0 (the same used in the example), and to the one for H using a chaining interest measure always equals to 1. Note

Fig. 6 . 3 Fig. 7 .

 637 Fig. 6. Replacing H 3 by H 3

Fig. 8 .

 8 Fig. 8. Dependency graphsIn summary, by marking classifiers to be replaced, the preselection measure P can reduce learned label relations to minimize the risk of error propagation. The learned relation structure can change according to the order of replaced classifiers, which is determined by the selection measure S. The chaining interest measure I controls allowed label relations: predictions are supposed to be more consistent by allowing replaced classifiers to learn many relations, and by doing the opposite the risk of error propagation is reduced.

  ± 0.95 PSI-MC (I = 1) 0.1894 ± 0.0721 PSI-MC (I = 0) 0.1884 ± 0.0709

  3} is the set of membership grades, and {a 1 , a 2 } is the set of descriptive attributes.

		a 1	a 2	c 1	c 2	c 3	c 4	c 5
	x 1	20	20	0	0	3	0	0
	x 2	30	40	1	0	3	0	0
	x 3	20	30	0	0	3	0	0
	x 4	20	10	0	0	0	0	3
	x 5	50	40	2	3	0	1	2
	x 6	50	20	2	3	0	1	2
	x 7	10	10	0	1	2	2	3
	x 8	10	30	0	3	1	2	2
	x 9	10	10	0	1	2	2	3
	x 10	10	50	0	3	1	2	2
				TABLE I
				GMLC DATA

TABLE II

 II

	describes datasets used in our experiments. The
	original BelaE data

TABLE III EVALUATION

 III OF GRADED MULTI-LABEL CLASSIFIERS USING BELAE

DATASET

Our proposed meta-classifier is at least better than Full CLR according to hamming-loss (TABLE

TABLE V THE

 V IMPACT OF CHAINING INTEREST MEASURE ON CLASSIFIER COMPLEXITY VI. CONCLUSIONS AND FUTURE WORK

978-1-5090-5781-8/16/$31.00 2016 IEEE