Adrien Zerbinati
email: adrien.zerbinati@inria.fr

Jean-Antoine Désidéri
email: jean-antoine.desideri@inria.fr

Régis Duvigneau
email: regis.duvigneau@inria.fr

Comparison between two multi objective optimization algorithms : PAES and MGDA. Testing MGDA on Kriging metamodels

In multi-objective optimization, the knowledge of the Pareto set provides valuable information on the reachable optimal performance. A number of evolutionary strategies (PAES [4], NSGA-II [3], etc), have been proposed in the literature and proved to be successful to identify the Pareto set. However, these derivativefree algorithms are very demanding in computational time. Today, in many areas of computational sciences, codes are developed that include the calculation of the gradient, cautiously validated and calibrated. Thus, an alternate method applicable when the gradients are known is introduced presently. Using a clever combination of the gradients, a descent direction common to all criteria is identified. As a natural outcome, the Multiple Gradient Descent Algorithm (MGDA) is defined as a generalization of the steepest-descent method and compared with PAES by numerical experiments. Using MGDA on a multi objective optimization problem requires the evaluation of a large number of points with regard to criteria, and their gradients. In the particular case of CFD problems, each point evaluation is very costly. Thus here we also propose to construct metamodels and to calculate approximate gradients by local finite differences.

Introduction

The numerical treatment of a multi-objective minimization is usually aimed to identify the Pareto set or a convenient subset of it. In the literature, several authors have proposed to achieve this goal by various algorithms, each one adapting a particular Evolution Strategy (ES). Such approaches are compared in the book of Deb [START_REF] Deb | Meyarivan : A Fast and Elitist Multi-Objective Genetic Algorithm : NSGA-II[END_REF]. Using a sufficiently diverse initial sample, these methods produce a discrete set of 2 by 2 non-dominated points. However, the most commonly used methods are very demanding in terms of computational time, as ES do in general.

In the particular case in which the gradients of the objective functions are at reach, at the current design point, faster algorithms can be developed. In the convex hull of the gradients of the objective functions, a direction exists along which all criteria diminish [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF]. The MGDA results in utilizing this direction as search direction and optimizing the stepsize appropriately. In this way, the classical steepest-descent method is generalized to multi-objective optimization. Applying MGDA thus corresponds to a phase of cooperative optimization.

In section 2, theoretical aspects leading to MGDA are briefly recalled. A complete presentation is available in [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF]. In section 3, results of a numerical experimentation on a classical test case are presented and commented.

Theoretical aspects

Cooperative-optimization phase : Multiple-Gradient Descent Algorithm (MGDA)

Here, to be complete, we review briefly the notions developed in [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF]. The general context is the simultaneous minimization of n (n ∈ N) smooth criteria (or disciplines) J i (Y ) (Y : design vector, Y ∈ R N ). Starting from an initial design point that is not Pareto optimal, a cooperative optimization phase is defined that is beneficial to all criteria.

Pareto concepts

Following [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF], we introduce the notion of Pareto stationarity: a design point Y 0 is said to be Pareto stationary if there exists a convex combination of the gradients of the smooth criteria J i that is equal to 0 at this point. Thus :

Definition 1. The smooth criteria J i (Y ) (1 ≤ n ≤ N)
are said to be Pareto stationary at the design point Y 0 if:

• ∀i = 1, .., n, u 0 i = ∇J i Y 0 ; • ∃ (α i ) i=1,..,n , α i ≥ 0, n ∑ i=0 α i = 1, n ∑ i=0 α i u 0 i = 0.

Existence and uniqueness of the minimal-norm element

Consider a family of vectors, denoted (u i ) i∈I , 1 ≤ i ≤ n. The following lemma holds :

Lemma 1 (Existence and uniqueness of the minimal-norm element). Assume :

• {u i }(1 ≤ i ≤ n) a family of n vectors in R N ;
• U be the set of strict convex combinations of these vectors :

U = w ∈ R n /w = n ∑ i=0 α i u 0 i ; α i > 0, ∀i ; n ∑ i=0 α i = 1 .
and U its closure, or convex hull.

Then,

∃!ω ∈ U , ∀ ū ∈ U : ( ū, ω) ≥ (ω, ω) = ω 2 .
(The element ω exists since U is closed, and it is unique since U is convex; as a result, ∀ ū ∈ U , and ∀ε ∈ [0, 1], ω + ε(uω) ∈ U , and ω + ε(uω) ≥ ω , and this yields the conclusion [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF]).

In the case of two criteria, three configurations of the two gradients can be considered, as illustrated below 1: This result applies in particular to u i for all i. But, (u i , ω) is the Frechet-derivative of J i in the direction ω. Hence, if ω = 0, the Frechet-derivatives of all the criteria are bounded from below by the strictly positive number ω 2 . The direction -ω is therefore a descent direction common to all criteria. These considerations yield the following:

Theorem 1. Let J i (Y ) (1 ≤ i ≤ n ≤ N, N ∈ N)
be n smooth functions of the vector Y ∈ R N . Assume Y 0 is an admissible design-point. We denote u 0 i = ∇J i (Y 0 ) and :

U = w ∈ R N , w = n ∑ i=1 α i u 0 i ; ∀i, α i > 0; n ∑ i=1 α i = 1 (1)
Let ω be the minimal-norm element of the convex hull U , closure of U . Then :

1. Either ω = 0, and the criteria J i (Y ) (1 ≤ i ≤ n) are Pareto-stationary ;

2. Or ω = 0 and -ω is a descent direction common to all the criteria; additionally, if ω ∈ U , the inner product ( ū, ω) is equal to ω 2 for all ū ∈ U .

Based on these results, when the gradients of all the criteria can be computed, the following algorithm (MGDA) proceeds by successive steps that are beneficial to all criteria. In the practical implementation, one specifies a tolerance ε T OL on ω below which the linesearch is not performed.

Convergence of the MGDA

Provided that the criteria are formulated to be smooth, positive and infinite at infinity, the sequence of iterates produced by the MGDA has been proved to admit a subsequence converging to a Pareto-optimal point [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF]. One main purpose of this report is to illustrate this convergence by numerical experiments using testcases of variable complexity.

Practical determination of the vector ω

In the general case (n > 2), ω can be calculated by numerical minimization of the quadratic form that expresses ω 2 in terms of the coefficients {α i } of the con- vex combination, subject to the inequality constraints α i ≥ 0 (∀i), and the linear equality constraint ∑ i α i = 1. Many routines are effective to perform this optimiza- tion, for instance certain evolution strategies. However, the problem may become ill-conditioned for large dimensions.

However, in the particular case of two objectives, ω can be expressed explicitly. Recall Figure 1, for which u = u 1 = ∇J 1 and v = u 2 = ∇J 2 . In this figure, the gradient vectors, elements of R N are represented as vectors of R 2 with same origin O. This results in no loss of generality since only the norms of the two vectors, and the angle between them do matter. Eliminating the trivial case in which u = v (for which ω = u = v), the convex hull is then represented by the segment uv connecting the extremities of these representative vectors. Let ω ⊥ be the vector whose origin is O, and extremity is the orthogonal projection of O onto the line that supports the segment uv (convex-hull). If the vector ω ⊥ is in the convex hull, that is, if its representative points on the segment uv, it is ω; otherwise, ω is the vector of smallest norm between u and v. Thus let:

ω = (1 -α)u + αv (2)
and compute α ⊥ for which the above convex combination is orthogonal to uv, that is :

α ⊥ = (u, u -v) (u -v, u -v) If α ⊥ ∈ [0, 1], α = α ⊥ ; otherwise, α = 0 or 1, that is, ω = u or v, depending on whether α ⊥ < 0 or > 1.

Line-search

This part deals with the determination of the step length (line-search). In multi criterion optimization, it is not easy to compute a satisfactory step with respect to all criteria producing a significant evolution. An adaptative method to compute a satisfactory step for each multiobjective problem would be convenient.

At the current design point, the Frechet-derivatives of all the criteria are strictly negative (and equal if ω ∈ U ). For each criterion, a surrogate quadratic model is constructed after computing three function values, and a related optimum stepsize ρ i is calculated corresponding to the location of the ith surrogate model's minimum.

ρ = min i,1≤i≤n ρ i
The vector ω is such that, ∀i, ρ i ≥ 0 and ρ ≥ 0. Whenever ρ = 0, MGDA is interrupted.

Numerical experimentation

In this section, we conduct numerical experiments to demonstrate the convergence of MGDA to Pareto optimal solutions, and to compare this algorithm with PAES [START_REF] Knowles | Approximating the nondominated front using the Pareto Archived Evolution Strategy[END_REF].

Fonseca test case

This testcase corresponds to the two-objective unconstrained minimization of the functions

               f 1 (x) = 1 -exp - 3 ∑ i=1 x i - 1 √ 3 2 f 2 (x) = 1 -exp - 3 ∑ i=1 x i + 1 √ 3 2
The design variable is x = (x 1 , x 2 , x 3 ) ∈ R 3 . This testcase is known to yield a continuous but concave Pareto set in function space. Here, the Pareto set is not known analytically, but has been well identified by Deb using the well-known genetic algorithm NSGA-II [START_REF] Deb | Meyarivan : A Fast and Elitist Multi-Objective Genetic Algorithm : NSGA-II[END_REF]. To obtain an accurate discrete representation of the Pareto set by MGDA, we have applied the method starting from a set of some 50 initial design points located on a sphere in the design-space (Figure 3). In all cases, MGDA converges and provides an accurately defined point on the Pareto set (Figure 4).

In the next experiment, we have first applied PAES twice, each time starting from a different design point and generating 50 others. Then the remaining dominated design points have been discarded. Thus less than one hundred design points have been archived. This set is compared on Figure 5 with the result of applying MGDA starting from 12 well-distributed initial design points, so that the number of function evaluations is the same in the two cases. MGDA again produces design points closer to the Pareto set (improved accuracy), but here in a fewer number.

However, at identical computational cost, generally, PAES introduces more diversity in the final result. Thus it appears interesting to combine the accuracy of MGDA with the robustness of PAES in a hybrid method. To check this, we have used the two methods sequentially: PAES first to generate 15 design points, retaining 8 nondominated design points, then used as intial points for MGDA. In each case about 3 to 4 iterations are sufficient to converge and produce the accurate result indicated on Figure 6.

Applying MGDA on a Kriging metamodel

In this section, we conduct numerical experiments to demonstrate the convergence of MGDA to Pareto optimal solutions in conjunction with Kriging metamodels. A first Kriging metamodel is constructed with an initial database. From each initial point, MGDA yields to a better point used subsequently to update the metamodel.

- 6 (Fonseca test case) First step with a large PAES followed by MGDA iterates on each non dominated point found. Design space on the left, functional space on the right.

Kur test case

This testcase corresponds to the two-objective unconstrained minimization of the functions

g 1 (x) = - 2 ∑ i=1 -10 exp -0.2 x 2 i + x 2 i+1 , g 2 (x) = 3 ∑ i=1 |x i | 0.8 + 0.5 sin(x 3 i )
The design variable is x = (x 1 , x 2 , x 3 ) ∈ R 3 . This testcase is known to yield a nonconvex discontinuous Pareto set in function space. Two generations of non dominated points applying PAES from different initial configurations gives a good discrete approximation of the Pareto front obtained by Deb [START_REF] Deb | Meyarivan : A Fast and Elitist Multi-Objective Genetic Algorithm : NSGA-II[END_REF]. Figure 7 shows that the Pareto set is here discontinuous, especially in the design space, where three distinct groups of points are evident. In the next experiment, we have first applied PAES once from one initial design point to generate 100 new points. PAES sorts out 11 non dominated points from these 100. For each point obtained, MGDA produces a new one closer to the Pareto front, as illustrated by Figure 8.

Because of the sine in the second function, this test case is a multi modal problem. Thus optimization algorithms based on Gradient descent methods have experience difficulties. To asses the MGDA, a clever strategy must be adopted to generate a sufficiently diverse set of initial points. Presently we use an initial small and diverse set of design points forming a sample of a latin hypercube. This set gives a Kriging metamodel on which MGDA drives each initial point to a better one in terms of function values. If the MGDA points are sufficently widespread, a new metamodel is constructed with the initial set augmented. Whenever a new point is found too close to another one in the database, it is not considered to update the metamodel. In a few iterations of this method, best points obtained are close to the Pareto front. The following experiment (Figure 9) is based on a set of 10 initial design points from [-5, 5] 3 evaluated with respect to g 1 and g 2 , after 10 iterations of the process described above. The computational cost corresponds to 43 point evaluations, including the database but not the metamodel construction.

CFD test case

The last experiment is an optimum-shape design in compressible aerodynamics. The transonic flow (M ∞ = 0.83, α = AoA = 2 o ) about a generic aircraft wing is simulated by the solution of the 3D Euler equations by an upwind finite volume method over an unstructured mesh of some 200,000 points generated by the software GMSH [START_REF] Geuzaine | a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The cross sections of the wing are made homothetic with a linear variation in the spanwise direction. Thus only the shape of these sections, an airfoil, is optimized. This airfoil is represented by 7 B-spline functions for the upper surface, and 7 other ones for the lower surface. The leading and trailing edges are fixed, and this permits us to introduce a total of 10 geometrical design variables. Intially, these variables are set to define a cross section close to the classical NACA0012 airfoil.

MGDA is used here to solve the two-criterion optimization problem consisting of maximizing the lift coefficient and minimizing the drag coefficient simultaneously, starting from the specified initial geometry.

An initial set of 40 design points forming a sample of a latin hypercube in R 10 haa been considered. This first set of datapoints is employed for two purposes. First, it is used to construct initial Kriging metamodels of both functions (lift and drag). Second, it is used throughout the following cycle to provide starting points to initiate the MGDA iteration in different conditions. This iteration is conducted until convergence using at every iteration, gradients that are calculated by local finite-differences of the metamodels. Each converged point belongs to the Pareto set associated with the two-criterion problem related to the metamodels. It is then reealuated by a flow computation, and added to the database unless it is found too close to an existing point. At completion of this database enrichment process, the metamodels are updated, and this completes the cycle. In practice, in what follows, only two cycles were performed.

Figure 10 represents the convergence of MGDA from a particular initial database point. The figure indicates the converged point and the point obtained by the same design reevaluated by an Euler flow computation (actual lift and drag).

Figure 11 represents the initial database of 40 points and the ultimate database. With only 95 calls to the CFD solver, a significant improvement of both criteria is achieved and visibly, an approximate Pareto front begins to form.

Figure 12 represents pressure fields on the wing and the symmetry plane, corresponding to 3 particular non dominated points of the ultimate database (points 1 and

Conclusion

In this article, we have tested by numerical experiment a recently proposed gradientbased algorithm for multiobjective optimization, MGDA [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF].

First, the convergence to Pareto-optimal solutions has been demonstrated in an analytical testcase corresponding to a continuous, but concave Pareto front (Fonseca test case). Additional information on this comparison can be found in [START_REF] Zerbinati | Comparison between MGDA and PAES for multi objective optimization[END_REF].

Second, MGDA has been compared with the well-known PAES algorithm. Both the Fonseca and Kur testcases have been considered in this comparison. We found that the two algorithms have particular merits of their own. PAES is very effective to converge to a very diverse dataset, whereas MGDA achieves this only if the initial set of design points is itself diverse. However, the iterative convergence of MGDA which makes use of (approximate) gradients is much faster. Thus both algorithms are complementary.

Thirdly, a hybrid method has been proposed and tested over the above mathematical testcases, demonstrating promising potentials.

Lastly, in the context of a two objective aerodynamic wing shape optimization in which the 3D Euler equations have been solved, MGDA has been used to define a strategy to progressively enrich the database associated with metamodels of drag and lift. With less than 100 calls to the flow solver, both lift and drag have been improved significantly from an initial design of a wing whose cross section was close to the classical NACA0012 airfoil.

Fig. 1

 1 Fig. 1 Various possible configurations of the two gradient-vectors u = u 1 and v = u 2 and the minimal-norm element ω.

Fig. 2

 2 Fig. 2 Variation of normalized the cost functions with the stepsize ρ in -ω direction.

Fig. 3

 3 Fig.[START_REF] Deb | Meyarivan : A Fast and Elitist Multi-Objective Genetic Algorithm : NSGA-II[END_REF] Convergence of MGDA to the Pareto front, for several initial design points., in design space (x, y, z) (left) and in function space ( f 1 , f 2 ) (right)

Fig. 4 (Fig. 5 (

 45 Fig. 4 (Fonseca test case) Convergence of MGDA to the Pareto front, for several initial design points.

Fig. 7 (

 7 Fig. 7 (Kur test case) Discrete Pareto front produced by 2 generations of PAES optimization.

Fig. 8 (

 8 Fig. 8 (Kur test case) Applying MGDA to each non dominated point from one PAES generation of 100 points (sort 11 non dominated).

Fig. 9 (

 9 Fig. 9 (Kur test case) Evolution of points gived by MGDA on an evoluting Kriging metamodel. 10 initial design points lead to 43 points

Fig. 10 (

 10 Fig. 10 (Eulerian flow test case) Example of the convergence of MGDA from an initial database point on the metamodel and the corresponding simulation result point.

Fig. 11 (

 11 Fig. 11 (Eulerian flow test case) Evolution of datapoints by MGDA applied to Kriging metamodels of lift and drag; the dataset is made of 40 design points initially, and 95 ultimately. All points are Eulerian simulation result.

Fig. 12

 12 Fig. 12 Pressure field associated with design points 1, 2 and 3 of Figure 11


John-David Yoder

Jeffrey West

Eric Baumgartner

Mathias Perrollaz

Michael Seelinger

Matthew Robinson

Experiments comparing precision of stereo-vision approaches for control of an industrial manipulator Authors

Keywords:

Motivation

Since the 1970s, researchers have been proposing ways to use vision in order to solve this problem. While the purpose of this paper is not to provide a complete review of visionbased robotic manipulation, an excellent overview of early work is found in [START_REF] Corke | Visual Control of Robot Manipulators, --A review[END_REF]. There has been success in the application of such technologies, especially in 2-D and 2.5-D problems [START_REF] Morel | System helps prepare sintered metal parts[END_REF]. Despite the fact that the theory for the solution of this problem has been established, there are limited implementations of full 3-D applications. Basically, the reasons for this come down to "the devil is in the details" and, in real 3-D applications, the accuracy, robustness, and cost-effectiveness of vision-based systems have been insufficient to justify widespread use.

A variety of efforts have been put forth to solve the problem of vision-based robotic manipulation. Visual Servoing offers promise by "closing the loop" in the image plane [START_REF] Hutchinson | A tutorial on visual servo control[END_REF]. However, in many assembly tasks, the robot end effector obscures the target when it gets close to task completion. Several methods have been proposed to get past this limitation, such as outfitting both the robot and the workpiece with fiducials that remain visible [START_REF] Hammer | An autonomous mobile manipulator for assembly task[END_REF], or "eye-in-hand" in cases where it is possible to place the camera(s) in the robot end effector. Another limitation of visual servoing is that in general it provides the least accuracy in the direction of the focal axis. Camera Space Manipulation [START_REF] Skaar | Three-dimensional camera space manipulation[END_REF] [START_REF] Gonzalez-Galvan | Application of a Precision-Enhancing Measure in 3D Rigid-Body Positioning Using Camera-Space Manipulation[END_REF] offers an alternative methodology and has been shown to allow completion of high-precision 3-D tasks. Unfortunately, this method does not allow standard stereo cameras to be used, as it requires widely separated and highly vergent cameras and utilizes the simplified orthographic model (though later work created accuracy similar to the pinhole model) [START_REF] Gonzalez-Galvan | Efficient camera-space manipulation using moments[END_REF].

2

The Mars Exploration Rovers used stereo vision with calibration for placing instruments on rock and soil targets [START_REF] Baumgartner | The Mars Exploration Rover Instrument Positioning System[END_REF], but this showed some limitations in accuracy [START_REF] Robinson | Hybrid image plane/stereo (HIPS) manipulation for robotic space applications[END_REF]. Recently, in space robotics applications, two additional approaches have been offered, HIPS [START_REF] Robinson | Hybrid image plane/stereo (HIPS) manipulation for robotic space applications[END_REF] and AGATE [START_REF] Yoder | Long-range autonomous instrument placement[END_REF] [START_REF] Seelinger | Autonomous Go-And-Touch Exploration (AGATE)[END_REF]. Both of these approaches have shown promise on space-related platforms to produce high-precision, vision-based manipulation using stereo cameras, with an application of instrument placement. These papers also reported large numbers of experiments. In particular, [START_REF] Robinson | Hybrid image plane/stereo (HIPS) manipulation for robotic space applications[END_REF] showed the improved accuracy using HIPS compared to precalibrated stereo on NASA testbeds. However, these testbeds typically have inaccurate kinematics (backlash, inconsistent zero offsets, lack of rigidity, etc.).

Technical Approach

In response to a request from a local manufacturer, several of the authors were involved in a project to use computer vision to automate the teaching of a robot task for parts with small lot sizes. The details of that task are found in [START_REF] Perrollaz | Teachless Teach-Repeat : Toward Vision-based Programming of Industrial Robots[END_REF]. While the corporation was pleased with the results of the project, the accuracy obtained was not sufficient for all of their parts. As such, the authors have begun a detailed investigation of the accuracy obtained, and how that compares to other approaches. It should be noted that the industry-focus of this project has influenced the technical approach.

The goal of this paper is to answer the question of whether techniques such as AGATE can create significant improvement in accuracy when applied to an industrial robot by using local data to modify the parameters of the camera-robot model. Also, while AGATE was developed for mobile manipulation (as it can control the mobile base and the manipulator relative to a visual target), this paper will use AGATE techniques for a fixed-base manipulator. The paper further includes the first direct comparison of AGATE with HIPS. Future work will compare these methods directly to visual servoing and traditional camera calibration techniques.

HIPS and AGATE both make use of the CAHVOR camera model [START_REF] Gennery | Least-Squares Camera Calibration Including Lens Distortion and Automatic Editing of Calibration Points[END_REF]. In both cases, the cameras can be calibrated using a standard checkerboard, but are typically calibrated by moving the manipulator through a series of poses in camera-space. In summary, a least-squares minimization is completed to compute the CAHVOR parameters based on data sets made up of 3-D robot positions and the 2-D image-plane appearance of the end effector in each camera. This minimization is completed after the manipulator moves through a pre-programmed set of positions designed to cover a significant portion of the robot's workspace and view of the camera(s). The CAHVOR parameters are then updated with additional measurements when they are available. Such measurements can be obtained while the end effector is approaching the target point. Addition of the localized samples results in a mapping between image space and physical space that is not globally accurate, but is very accurate relative to the manipulator in the region of the target point. This general approach has been shown to work well even in the presence of large kinematic errors [START_REF] Robinson | Hybrid image plane/stereo (HIPS) manipulation for robotic space applications[END_REF] [START_REF] Seelinger | Autonomous Go-And-Touch Exploration (AGATE)[END_REF]. Both HIPS and AGATE have been shown to provide the following advantages compared to traditional calibration methods: Robustness to poor kinematics, the ability to deal with changes to the internal and external parameters of the camera, and the ability to achieve high precision manipulation relative to visually specified points. In addition, AGATE has been shown to be able to control mobile manipulators, and use more than two cameras [START_REF] Seelinger | Autonomous Go-And-Touch Exploration (AGATE)[END_REF].

Robinson [START_REF] Robinson | Hybrid image plane/stereo (HIPS) manipulation for robotic space applications[END_REF] showed the advantage of HIPS (which finds the CAHVOR model parameters based on manipulator observations of fiducials on the manipulator) over precalibrated stereo, and the advantage of adding the local informationboth in simulation and in a large number of experiments with a mobile manipulator. However, the manipulator in question was not as accurate as a typical industrial manipulatorit was less rigid and had significant kinematic errors (it was a prototype on a mobile manipulator testbed). The primary goal of this study is to see if a similar improvement takes place using an industrial manipulator, or whether the improved rigidity and kinematic accuracy negates any improvement due to on-the-fly CAHVOR model updates. For consistency with [START_REF] Robinson | Hybrid image plane/stereo (HIPS) manipulation for robotic space applications[END_REF], we will call the version of the code that does not update the model during approach "Static AGATE" and the version that does update "Dynamic AGATE".

Experiments

Experiments were conducted with the 6-axis ABB IRB-140 robot shown in Figure 1. The vision sensor is a pair of PointGrey Flea2 cameras, delivering 640 x 480 RGB images at 15FPS. As shown in Figure 2, Cameras are placed on a common portable support with a baseline separation of 15cm, and are placed approximately 2.5m from the workspace. While the cameras are effectively parallel, no particular efforts were taken to ensure a proper rectified configuration. Central to HIPS or AGATE is the ability to locate the end effector of the robot as it moves through the workspace. These image-plane locations, along with the corresponding location of the robot, are used to update the CAHVOR camera parameters.

To simplify this task, the robot's end-effector is equipped with a LED blinking system as shown in Figure 3. The blinking frequency is roughly a third of the frequency of the cameras to avoid possible aliasing. We chose a rather simple design of the LED lamp, which does not require synchronization between the emitter (LED) and the receiver (camera). This is possible because the arm can be static during the acquisition of a set of images. Once the images are acquired, two images, respectively representing the mean value and the standard deviation over time of the intensity of each pixel, are computed.

The cyan component of both images is computed using the green and blue planes (because our LEDs produce cyan rather than blue light). Then connected components of high standard deviation are extracted. These regions are detected as LEDs if their mean value over the cyan image is also sufficiently high.

All control is done by the ABB's controller. A directly-connected serial connection is used to transfer commands and coordinates to the ABB. More details about the system, and the implementation of a precise manipulation task, can be found in [START_REF] Perrollaz | Teachless Teach-Repeat : Toward Vision-based Programming of Industrial Robots[END_REF].

A typical experiment consists of placing a fiducial in the workspace. For Dynamic AGATE, the robot is commanded to move half the distance to the target, and an additional measurement is made. If the results of this measurement are within normal bounds, this point is used to update the CAHVOR parameters. The process is repeated several times, until the robot is less than 4mm from the final goal, at which point it is moved to the final goal position. For Static AGATE, once the fiducial is found, the robot is moved directly to the target.

Results

To provide a measure of accuracy, the robot end effector was equipped with a standard whiteboard marker. The marker was placed in a holding device that provides approximately 1cm of compliance along the axis of the marker (shown in Figure 3). Circular fiducials were placed in the workspace, and the robot was instructed to approach the center of the fiducial and stop 5mm above the fiducial (above corresponds to the Z axis) of the base frame of the robot. The jog mode of the robot is then used to move straight down to mark the fiducial. The vertical travel during the jog is used to measure the error in the Z direction. The distance from the mark to the center of the fiducial circle can be used to measure the error in the X-Y plane in the robot base frame, since the fiducials are placed on this plane. The X-Y-Z coordinate frame of the robot is shown in Figure 4. A sample of the marked fiducial is given in Figure 5 To provide a sense of scale, the outer diameter of the black circle is 76mm.

These tests were repeated 7 times for both the Static AGATE and Dynamic AGATE with the fiducials placed at varying locations in the workspace. One immediate finding was that Static AGATE did not consistently touch the fiducial (the marker often stopped above the fiducial). In these cases, the distance from the marker tip to the fiducial was measured and recorded as an error in the z direction.

Thirty tests were conducted using Static AGATE and Dynamic AGATE. The tests consisted of 20 different physical locations throughout a region of the workspace (ten of the positions are reached twice). Table I shows a summary of these measurements. In summary, the mean error with Dynamic AGATE was 11.0mm, while the mean error with Static AGATE was 14.6mm. A t-test showed this to be significant at the 95% confidence level (t=-2.169, significance =0.038). Dynamic AGATE showed a standard deviation of 7.8mm, reduced from 9.9mm with Static AGATE.

Looking at the data in Table I, it is clear that most of the difference between the two is due to the XY-plane error. Comparing only this error, Dynamic AGATE resulted in a mean error (distance from the center of the cue) of 9.6mm compared to 13.5mm for Static AGATE. A t-test showed this to be significant at the 95% confidence level (t=-2.135, significance =0.041).

Using the data sets obtained during the calibration process and during the tests, a direct comparison was completed between AGATE and HIPS. Both methods were used to project the observed locations of the LEDs from camera-space into physical space. The methods use different approaches for this process. While full details are outside the scope of this paper (see [START_REF] Robinson | Hybrid image plane/stereo (HIPS) manipulation for robotic space applications[END_REF] and [START_REF] Seelinger | Autonomous Go-And-Touch Exploration (AGATE)[END_REF]) HIPS find the 3-D point at the center of the shortest line between the projecting rays. AGATE finds the 3-D point that minimizes the imageplane error between the measured and projected points. Despite these differing approaches, the results were very similar. The average difference (as measured by the norm of the vector between the 3-D projections) was 0.021mm, and the standard deviation was 0.025mm. Given the scale of the other errors involved, the conclusion is that AGATE and HIPS would produce the same accuracy of control.

repeat mode, this error would not cause any positioning errors. Because of the plan to use this in industry, the standard ABB functions were used to "teach" the ABB the location of the points of interest (the tip of the marker, and the point between the two LEDs). This involves taking the point being taught to the same point in physical space from four different directions. The ABB then completes a least-squares minimization to locate the point relative to the robot coordinate system. For the data set shown in Table I, the mean error for the marker was 0.73mm, while the mean error for the LED center position was 0.41mm.

Error #7, the correspondence error, is very important, and is related to the detection errors. Essentially this error is due to the fact that a given item (LED, fiducial, etc.) in X-Y-Z space will produce 4 values in camera space (assuming two cameras). If both cameras detect EXACTLY the same physical point, then there is no inconsistency between the camera-space measurements and the physical space measurements. However, this is impossible to do (as stated earlier, finding the center of the fiducial, for example, has a standard deviation of 0.2 pixels). An error of 1 pixel in the horizontal location of an item in one of the cameras, for instance, will create an error of over 20mm when it is projected back into physical space using the AGATE method. This was found to be identical using HIPS.

When computing the parameters of the CAHVOR camera model, many of these errors are included. Recall that this computation is based on the image-plane location of the LEDs and the corresponding physical location of the robot. Thus, errors #1-3, and #5-7 are included (all errors other than the marker kinematics and the fiducial errors). In our experiments over multiple calibrations, these lead to a mean reconstruction error in the CAHVOR model of approximately 5mm. The reconstruction error is found by using the CAHVOR parameters to project the LED locations back into 3-D space and comparing those locations to the actual, recorded 3-D locations of the robot.

Main Experimental Insights

The main experimental insight is that Dynamic AGATE outperforms Static AGATE, even on this system with high-quality cameras and a very repeatable industrial manipulator. Furthermore, AGATE and HIPS produce effectively the same results.

It should be noted that the Dynamic AGATE discussed to this point in the paper typically resulted in six or seven model updates. Experiments were conducted to see if additional updates would further improve the accuracy. In order to do this, the approach was modified so that the robot moved one third of the way to the target each time rather than half way. This resulted in, on average, five additional updates per experiment. Another 30 experiments were conducted in this manner. However, no statistically significant change was observed. Thus it was concluded that this number of updates slowed the system without a significant advantage in accuracy.

The accuracy is not as high as that shown in previous work with AGATE and HIPS.

Examining the sources of error described above, there are a variety of reasons for this increased error. Error #1 is much smaller than other systems tested, and errors #2 and #3 are comparable. Error #6, however, is larger. In previous work with HIPS and AGATE, highly accurate kinematic models of the end effectors (specifically the relationship between the end effector and the fiducials) were used. It is more realistic for an industrial tool which will have to be changed regularly to have an operator use the tool-teaching algorithms built-in to the ABB. However, this introduces more error. Taking the mean of errors #1, #3, and #6 would predict a mean error of approximately 2mm.

Errors #4, #5, and #7, however, should be expected to be larger in the experiments described in this paper than in previous work using HIPS or AGATE. This is because the distance from the cameras to the workspace is much greater (2.5m in this case compared to approximately 0.5m). This leads to increased error in several ways. First, it increases the error in locating the target and the LEDs in the image plane. Second, it increases the size of a voxel. Thirdly, and most importantly, it increases the important of correspondence error. The importance of this effect was shown in [START_REF] Baumgartner | Sensitivity Analysis for a Remote Vision-Guided Robot Arm under Imprecise Supervisory Control[END_REF]. But in that set of experiments, the error shows was only 1.5mm/pixel. In the experiments shown in this paper, the error was over 20mm/pixel. Given the magnitude of these errors, showing an overall error approximately 11mm is consistent with expectations.

It is also worth noting that a major advantage of this methodology is that it allows the placement of the stereo rig on a simple tripod. This tripod can be moved relative to the robot without impacting the performance of the system, since the system can re-estimate the CAHVOR camera model at any time. In industrial settings, the cameras would certainly be fixedbut this would provide the ability to adjust if the cameras are bumped or moved during operation or maintenance activities. The method furthermore does not require that the images are rectified, and the alignment of one camera relative to the other is arbitrary, as long as the workspace can be seen by both cameras.

Figure 1 Figure 1 :Figure 4 :

 114 Figure 1: Experimental Setup Figure 2: Camera System

Figure 5 :

 5 Figure 5: Marked fiducial

Acknowledgment

The authors would like to thank American Trim Corporation for their support of this project. Furthermore, thanks are due to Amber Cool and Patrick Whitten, undergraduate students without whose assistance we could not have gathered sufficient data for this paper. Further thanks to Dr. Sami Khorbotly, who developed the hardware for the LED system.

Understanding the Source of Positioning Errors

To draw any conclusions from the data presented, there must be an understanding of the sources of error that are included in the errors reported in Table I. The basic sources of error will be enumerated here: 1) Errors due to the limitations on the mechanical system. That is, a motor can only be commanded to within one encoder countthis limitation on the precision of control will manifest itself as an error in the final measurement. 2) Error due to inaccuracies in the models and the limitations of the algorithms.

3) Error in making the measurement. 4) Error in the detection of the target. 5) Error in the detection of the LEDs. 6) Error due to the kinematic modeling of the LEDs and the marker. 7) Error in ensuring that the same physical point is selected in both cameras (correspondence error).

While these are all sources of error, error #1 is very small for an industrial robot. ABB reports repeatability, not accuracy, of their robots. The repeatability of the ABB 140 used in these tests is ± 0.03mm [14], which is insignificant compared to all other errors.

Error #2 is related to the approach, and will be discussed later.

Error #3 has been minimized by physically marking the fiducial for later measurement. However, since the marker makes a point that has a diameter of over 1mm, it is likely that an error on the order of 1mm is reasonable for this item. A detailed examination of this could be undertaken by having multiple people measure the location of each mark, but this was not done for the experiments described here.

Using a fiducial limits the magnitude of error #4, the detection of the target. Multiple measurements were taken (though not a statistically significant number), and it was found that the standard deviation on the location of the center of the fiducial was approximately 0.2 pixels. While the volume of a voxel varies throughout the workspace, it is approximately 1.67x5.00x1.67mm (in X, Y, Z, respectively). Thus, simply due to the resolution, one should expect errors with a standard deviation of approximately 1mm. To mitigate the effects of this error, the fiducial location was found in ten consecutive images, and the average center location was used.

Regarding error #5, the LEDs are detected using image-differencing in order to reduce the likelihood of false positives. A center-of-gravity approach is used to find the center of each LED, and the two LED locations are combined into an "average" of their locations in order to reduce errors and eliminate the ambiguity between the two. Multiple tests were run and the standard deviation for the detection of the LEDs was found to be negligible.

Error #6 would always be presentthe kinematics of any tool that is added on to the robot must be modeled and would always include error. Note that when used in teach-

Continuing Work

While the work to date has shown that Dynamic AGATE shows improvement over Static AGATE, and that AGATE and HIPS produce effectively the same result, much work remains to be done. In particular:  Build a more accurate tool holder to improve the accuracy of measurements. It is unclear how much of the current error is due to the methodology, and how much is due to inaccuracy in the tool holder. As can be seen in Figure 3, the current setup is simply clamped on, in a manner that is neither accurate nor repeatable.  Comparison with traditional CAHVOR stereo (using a checkerboard for calibration)  The use of more than two cameras (which is directly supported in AGATE)  Comparison with an "off the shelf" disparity-based stereo system.  Comparison with visual servoing. When complete, it is expected that this will provide a rich dataset of direct comparison of four different, previously-published approaches to the problem of visually-guided manipulation. The one approach not specifically examined is Camera-Space Manipulation. While this approach has shown very good results, it requires widely separated and highly vergent cameras, and therefore would require a completely different experimental setup than the other methods. This would make a direct comparison difficult.


Paola Goatin

TRAFFIC FLOW MODELS WITH PHASE TRANSITIONS ON ROAD NETWORKS

Keywords: 2000 Mathematics Subject Classification. 90B20, 35L65. Key words and phrases. Hyperbolic Conservation Laws, Riemann Problem, Phase Transitions, Continuum Traffic Models

.

1. Introduction. The interest in traffic dynamics has considerably increased in the last decades. The modelling of pedestrian and vehicular traffic can be developed under different approaches. We can distinguish between microscopic (particle-based), mesoscopic (gas-kinetic) and macroscopic (fluid-dynamic) models. We refer the reader to the review paper [START_REF] Helbing | Traffic and related self-driven many-particle systems[END_REF] for an overview on the possible approaches, and the analysis and interpretation of various interesting phenomena occurring in traffic. A survey of the available mathematical models is given in [START_REF] Bellomo | On the mathematical theory of vehicular traffic flow I -fluid dynamic and kinetic modelling[END_REF][START_REF] Garavello | Traffic flow on networks[END_REF][START_REF] Klar | Mathematical models for vehicular traffic[END_REF].

The fluid-dynamic approach considers the evolution of macroscopic variables, such as the density of vehicles and their average velocity. Historically, one of the first continuous models introduced to describe traffic flow is the well known Lighthill-Whitham [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF] and Richards [START_REF] Richards | Shock waves on the highway[END_REF] (LWR) model, which reads

∂ t ρ + ∂ x [ρv(ρ)] = 0, (1)
where ρ ∈ [0, R] is the mean traffic density, and v(ρ), the mean traffic velocity, is a given non-negative non-increasing function. The maximal density R > 0 corresponds to a traffic jam. This scalar model expresses conservation of the number of cars, and relies on the assumption that the car speed depends only on the density (more complex closure relations between speed and density, involving the density gradient, can be assumed, see [START_REF] Bellomo | On the mathematical theory of vehicular traffic flow I -fluid dynamic and kinetic modelling[END_REF] and references therein). This phenomenological relation is valid in steady state conditions, and is not realistic in more complicated situations. In particular, as shown in Figure 1, the corresponding fundamental diagram in the (ρ, ρv)-plane does not qualitatively match experimental data at high densities.

Later on, several second order models, i.e. models with two equations, were considered, see [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF][START_REF] Greenberg | Congestion on multilane highways[END_REF][START_REF] Payne | Models of freeway traffic and control[END_REF][START_REF] Siebel | On the fundamental diagram of traffic flow[END_REF][START_REF] Zhang | A non-equilibrium traffic model devoid of gas-like behavior[END_REF]. A third order model was presented in [START_REF] Helbing | Improved fluid-dynamic model for vehicular traffic[END_REF].

The diagram showed in Figure 1 suggests that a good traffic flow model should exhibit two qualitative different behaviors:

• for low densities, the flow is free and essentially analogous to that of the LWR model; [START_REF] Kerner | Phase transitions in traffic flow[END_REF]; here q denotes the flux ρv.

• at high densities the flow is congested and covers a 2-dimensional domain in the fundamental diagram; a "second order" model seems more appropriate to describe this dynamic.

Traffic flow models with phase transitions have been considered in the literature since the 60-ties (see Helbing [29, Section II] for a description of the features recovered by a detailed analysis of the fundamental diagram). In particular, we refer the reader to the scalar model by Drake, Schofer and May [START_REF] Drake | A statistical analysis of speed density hypotheses[END_REF]. Another model has been introduced more recently by the author [START_REF] Goatin | The Aw-Rascle traffic flow model with phase transition[END_REF].

In the present paper we concentrate on the second order model with phase transitions introduced by Colombo [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF]. This model has been conceived in order to reproduce the empirical flow-density relations showed in Figure 1 (right). From the analytical point of view, the model is well posed in the space of functions with bounded total variation. More precisely, when considering the Cauchy or the Initial-Boundary value problem, it is possible to construct a Lipschitz continuous semigroup of solutions, which is defined on a domain of functions with bounded total variation, see [START_REF] Colombo | Global well posedness of a traffic flow model with phase transitions[END_REF].

Here we are interested in extending the theory to road networks. The results available for networks concern the LWR model or the Aw-Rascle model, see [START_REF] Bastin | A second order model of road junctions in fluid models of traffic networks[END_REF][START_REF] Chitour | Traffic circles and timing of traffic lights for cars flow[END_REF][START_REF] Coclite | Traffic flow on a road network[END_REF][START_REF] Garavello | Conservation laws with discontinuous flux[END_REF][START_REF] Garavello | Source-destination flow on a road network[END_REF][START_REF] Garavello | Traffic flow on a road network using the Aw-Rascle model[END_REF][START_REF] Helbing | Self-organized network flows[END_REF][START_REF] Herty | Optimization criteria for modelling intersections of vehicular traffic flow[END_REF][START_REF] Holden | A mathematical model of traffic flow on a network of unidirectional roads[END_REF]. No result was available for models with phase transitions up to now.

In [START_REF] Colombo | Road network with phase transitions[END_REF] authors prove the existence of weak solutions on the whole network for initial data of bounded variation under the (technical) assumption that traffic keeps away from the zero velocity. The construction is based on the wave-front tracking method: We first consider Riemann problems at nodes, which are Cauchy problems with constant initial data on each road converging to a given junction. Notice that the only conservation of cars is not sufficient to determine a unique solution. Thus one has to prescribe solutions for every initial data and we call the relative mapping a Riemann solver at nodes. In the case studied here, we consider two Riemann solvers that are defined by generalizing to the phase transition model two Riemann solvers previously presented for the LWR model: The Riemann solver R 1 J was proposed for vehicular traffic in [START_REF] Coclite | Traffic flow on a road network[END_REF], while the Riemann solver R 2 J was introduced for telecommunication networks in [START_REF] D'apice | Packet flow on telecommunication networks[END_REF]. The first prescribes a fixed distribution of traffic in outgoing roads, and then the maximization of the flux through the junction. The second maximizes the flux through the junction and then prescribes a distribution of traffic. In the case the Riemann solver R 1 J is considered, we are able to construct piecewise constant approximations via wave-front tracking algorithm (see [START_REF] Bressan | Hyperbolic Systems of Conservation Laws[END_REF] for the general theory and [START_REF] Garavello | Traffic flow on networks[END_REF] in the case of networks), using classical self-similar entropic solutions for Riemann problems inside roads and an assigned Riemann solver at junctions. To pass to the limit we rely on an estimate on the total variation of the flux. The paper is organized as follows. Section 2 is devoted to the description of the model, and the classical Riemann solver. Section 3 collects the results on the well posedness of the Cauchy and the Initial-Boundary Value problem. Road networks are introduced in Section 4, where we describe the sets of attainable values at junctions, and the Riemann solvers at junctions are described in Section 5.

2. Description of the model. We consider the model introduced in [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF]. It consists of a scalar LWR model coupled with the 2 × 2 system presented in [START_REF] Colombo | A 2 × 2 hyperbolic traffic flow model[END_REF]. The former applies to the states of free flow, while the latter to the congested states. A phase transition is a discontinuity separating a state of free traffic from one in the congested phase. More precisely, the model in [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF] reads

Free flow: (ρ, q) ∈ Ω f Congested flow: (ρ, q) ∈ Ω c ∂ t ρ + ∂ x [ρ • v] = 0 ∂ t ρ + ∂ x [ρ • v] = 0 q = ρV ∂ t q + ∂ x [(q -Q) • v] = 0 v = v f (ρ) = 1 -ρ R V v = v c (ρ, q) = 1 -ρ R q ρ . (2)
Here, R is the maximal traffic density, V is the maximal traffic speed and Q is a parameter of the road under consideration related to the phenomenon of wide jams, see [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF][START_REF] Kerner | The physics of traffic[END_REF]. The weighted linear momentum q is originally motivated by gas dynamics. It approximates the real flux ρv for ρ small compared to R.

It is assumed that if the initial data are entirely in the free (resp. congested) phase, then the solution will remain in the free (resp. congested) phase for all times. Thus we take Ω f and Ω c to be invariant sets for the corresponding equations. The resulting domain is given by Ω = Ω f ∪ Ω c , where

Ω f = {(ρ, q) ∈ [0, R] × [0, +∞[ : v f (ρ) ≥ V f , q = ρ • V } , Ω c = (ρ, q) ∈ [0, R] × [0, +∞[ : v c (ρ, q) ≤ V c , q-Q ρ ∈ Q --Q R , Q + -Q R .
Here, V f and V c are the threshold speeds, i.e. above V f the flow is free, and below V c the flow is congested. The parameters

Q -∈]0, Q[ and Q + ∈]Q, +∞[ depend on
the environmental conditions and determine the width of the congested region.

Figure 2 shows that the shape of the invariant domain is in good agreement with experimental data. Notice that the sets are represented in the (ρ, ρv)-plane.

Following [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF][START_REF] Colombo | Global well posedness of a traffic flow model with phase transitions[END_REF], throughout the present note we assume that the various parameters are strictly positive and satisfy

V > V f > V c , Q + ≥ Q ≥ Q -, Q + -Q RV < 1 , V f = V -Q + /R 1 -(Q + -Q)/(RV ) . ( 3
) ρv ρ 0 Ω f Ω c R Figure 2
. Left: invariant domain for [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF]. Right: experimental data, taken from [START_REF] Kerner | Phase transitions in traffic flow[END_REF]. The dotted straight lines exiting the origin are respectively ρv = ρV f and ρv = ρV c . The continuous curves that border Ω c are ρv

= (1 -ρ/R)(Q + ρ(Q ± -Q)/R). ρv ρ 0 Ω f Ω c R R - c u - c ũ w 1 w 2 0 V c V f W + 2 W - 2 W W o

Figure 3. Notation used in the paper

We recall the basic informations on the 2 × 2 system on the right hand side of (2):

r 1 (ρ, q) = ρ q -Q , r 2 (ρ, q) = R -ρ R ρ q , λ 1 (ρ, q) = 2 R - 1 ρ • (Q -q) - Q R , λ 2 (ρ, q) = v c (ρ, q) , ∇λ 1 • r 1 = 2 Q -q R , ∇λ 2 • r 2 = 0 , L 1 (ρ; ρ o , q o ) = Q + q o -Q ρ o ρ , L 2 (ρ; ρ o , q o ) = ρ ρ o R -ρ o R -ρ q o , w 1 = v c (ρ, q) , w 2 = q -Q ρ , (4)
where r i is the i-th right eigenvector, λ i the corresponding eigenvalue and L i is the i-Lax curve. In the Riemann coordinates (w

1 , w 2 ), Ω c = [0, V c ] × [W - 2 , W + 2 ]
. For (ρ, q) ∈ Ω f , we extend the corresponding Riemann coordinates (w 1 , w 2 ) as follows. Let ũ = (ρ, ρV ) be the point in Ω f defined by ρ = Q/(V -W - 2 ). Define

w 1 = V f and w 2 = V -Q/ρ if ρ ≥ ρ , v f (ρ) -v f (ρ) + V -Q/ ρ if ρ < ρ , (5)
so that, in the Riemann coordinates,

Ω f = {V f } × [W o , W + 2 ]
, see Figure 3. The 2 × 2 system describing the congested flow turns out to be hyperbolic, the second characteristic field is linearly degenerate but the first has an inflection point along the curve q = Q. Analogies between the solutions to (2) and real traffic features are given in [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF].

For notational convenience, we introduce the following short form

∂ t u + ∂ x f (u) = 0 (6)
for the model of phase transitions under consideration, with

u = (ρ, q) and f (u) = (ρv f (ρ), qv f (ρ)) , if (ρ, q) ∈ Ω f , u = (ρ, q) and f (u) = (ρv c (ρ, q), (q -Q)v c (ρ, q)) , if (ρ, q) ∈ Ω c .
2.1. The Riemann problem. We recall in this section the description of the classical Riemann solver for (2), i.e. the self-similar solution of the Cauchy problem

   ∂ t u + ∂ x f (u) = 0 , u 0 (x) = u l , if x < 0 , u r , if x > 0 . (7)
If the initial data u l , u r are in the same phase, standard Lax solutions to the corresponding Riemann problem can be considered. Otherwise, following [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF], admissible solutions are defined as follows.

Definition 2.1. If u l ∈ Ω f and u r ∈ Ω c , then an admissible solution to ( 7) is a self-similar function u : R × [0, +∞[ → Ω f ∪ Ω c such that, for some Λ ∈ R, we have:

1. u(] -∞, Λt[) ⊆ Ω f and u(]Λt, +∞[) ⊆ Ω c ; 2. the functions u -(x, t) = u(x, t) if x < Λt , u(Λt-, t) if x > Λt , (8)
u + (x, t) = u(Λt+, t) if x < Λt , u(x, t) if x > Λt , (9) (10)
are Lax solutions to corresponding Riemann problems for (2) left, right, respectively; 3. the Rankine-Hugoniot condition

ρ(Λt+, t) v c (u(Λt+, t)) -ρ(Λt-, t) v f (ρ(Λt+, t)) = Λ (ρ(Λt+, t) -ρ(Λt-, t))
holds for all t > 0.

If u l ∈ Ω c and u r ∈ Ω f , the conditions are obtained by interchanging the roles of Ω f , Ω c and v f , v c .

Notice that condition 3 above ensures that the total number of car is conserved across phase transitions. Definition 2.1 does not assure uniqueness. We are then led to introduce the notion of consistency [START_REF] Colombo | Hyperbolic phase transitions in traffic flow[END_REF]. Definition 2.2. Let R : (u l , u r ) → R(u l , u r ) denote a Riemann solver, i.e. x → R(u l , u r )(x) is the solution of (7) computed at time t = 1. R is consistent if the following two conditions hold for all u l , u m , u r ∈ Ω f ∪ Ω c , and x ∈ R:

(C1) R(u l , u m )(x) = u m R(u m , u r )(x) = u m ⇒ R(u l , u r ) = R(u l , u m ) , if x < x , R(u m , u r ) , if x ≥ x , (C2) R(u l , u r )(x) = u m ⇒        R(u l , u m ) = R(u l , u r ) , if x ≤ x , u m , if x > x , R(u m , u r ) = u m , if x < x , R(u l , u r ) , if x ≥ x .
Essentially, (C1) states that whenever two solutions to two Riemann problems can be placed side by side, then their juxtaposition is again a solution to a Riemann problem. (C2) is the vice-versa.

We are now ready to construct the Riemann solver. We consider several different cases:

(A) The data in [START_REF] Chitour | Traffic circles and timing of traffic lights for cars flow[END_REF] are in the same phase, i.e. they are either both in Ω f or both in Ω c . Then the solution is the standard Lax solution to (2), left, resp. right, and no phase boundary is present. (B) u l ∈ Ω c and u r ∈ Ω f . We consider the points u c ∈ Ω c and u m ∈ Ω f implicitly defined by

1 - ρ c R Q + w 2 (u l )ρ c = ρ c V c , 1 - ρ m R Q + w 2 (u l )ρ m = ρ m V 1 - ρ m R .
If w 2 (u l ) > 0, the solution is made of a 1-rarefaction from u l to u c , a phase transition from u c to u m and a Lax wave from u m to u r . If w 2 (u l ) ≤ 0, we have a shock-like phase transition from u l to u m and a Lax wave from u m to u r .

(C) u l ∈ Ω f and u r ∈ Ω c with w 2 (u l ) ∈ [W - 2 , W + 2 ]
. Consider the points u c and u m ∈ Ω c implicitly defined by

1 - ρ c R Q + w 2 (u l )ρ c = ρ c V c , 1 - ρ m R Q + w 2 (u l )ρ m = ρ m w 1 (u r ) .
If w 2 (u l ) > 0, the solution is made of a shock-like phase transition from u l to u m and a 2-contact discontinuity from u m to u r . If w 2 (u l ) ≤ 0, the solution displays a phase transition from u l to u c , a 2-rarefaction from u c to u m and a 2-contact discontinuity from u m to u r . (D) u l ∈ Ω f with w 2 (u l ) < W - 2 and u r ∈ Ω c . Let u m ∈ Ω c be the point on the lower boundary of Ω c implicitly defined by

1 - ρ m R Q + W - 2 ρ m = ρ m w 1 (u r ) ,
and consider the speed of the phase boundary joining

u l ∈ Ω f to u m ∈ Ω c Λ(u l , u m ) = ρ l v f (ρ l ) -ρ m w 1 (u r ) ρ l -ρ m . Let U c = (R c , Q c ) ∈ Ω c be the point whose Riemann coordinates are (V c , W - 2 ). If λ 1 (U c ) ≥ Λ(u l , U c
), the solution is a phase transition from u l to U c , a 1-rarefaction from U c to u m and a 2-contact discontinuity from u m to u r . Otherwise:

-If λ 1 (u m ) ≤ Λ(u l , u m ), the solution is a phase transition from u l to u m followed by a 2-contact discontinuity from u m to u r . -If λ 1 (u m ) > Λ(u l , u m ), let u c = (ρ c , q c ) ∈ Ω c be implicitly defined by

λ 1 (u c ) = Λ(u l , u c ) ,
i.e. ρ c is the bigger root of the equation

(Q -Q -)ρ 2 -2ρ l (Q -Q -)ρ + R 2 (ρ l v f (ρ l ) -Q) + ρ l R(2Q -Q -) = 0 and q c = Q -ρ c (Q -Q -)/R.
Then the solution shows a phase transition from u l to u c , an attached 1-rarefaction from u c to u m and a 2-contact discontinuity from u m to u r .

3. Well posedness. In the literature, several results deal with the solution to Riemann problems in presence of phase transitions, see for instance [START_REF] Colombo | Characterization of Riemann solvers for the two phase p-system[END_REF][START_REF] Lefloch | Nonclassical Riemann solvers and kinetic relations. I. A nonconvex hyperbolic model of phase transitions[END_REF][START_REF] Lefloch | Non-classical Riemann solvers and kinetic relations. II. An hyperbolic-elliptic model of phase-transition dynamics[END_REF]. Other works prove the global in time well posedness of the Cauchy problem, but with initial data that are perturbations of a given phase boundary, see for instance [START_REF] Colombo | Continuous dependence in conservation laws with phase transitions[END_REF][START_REF] Colombo | Sonic hyperbolic phase transitions and Chapman-Jouguet detonations[END_REF]. On the contrary, the results presented in this section do not require a priori bounds on the number of phase boundaries that are present in the data and in the solution. From the analytical point of view, this is a first example of a system of conservation laws developing phase transitions whose well posedness is proved globally, i.e. for all initial data attaining values in a given set and with bounded total variation.

From the traffic point of view, well posedness allows to consider various control and optimization problems, see [START_REF] Colombo | Minimising stop and go waves to optimise traffic flow[END_REF].

3.1. The Cauchy Problem. In this case, ( 6) is supplemented with a given value of the solution at time t = 0. More precisely, we assume that an initial datum u 0 ∈ Ω is given and we set u(., t = 0) = u 0 . (11) We introduce the notations:

X = L 1 (R; Ω) , TV(u) = TV(ρ) + TV(q) . (12
S t1 u 1 -S t2 u 2 L 1 ≤ L • ( u 1 -u 2 L 1 + |t 1 -t 2 |) ;
(RS4) if u ∈ D is piecewise constant, then for t small, S t u coincides with the gluing of solutions to Riemann problems.

By "solutions to Riemann problems" we refer here to those defined in Section 2. Properties (RS1)-(RS4) provide the natural extension of [6, Definition 9.1] to the present case.

We are now ready to state the main result of this section, namely the existence of an M -RS generated by the Cauchy problem for (2). The proof can be found in [15, § 4.2]. Observe that the description of several realistic situations requires suitable source terms in the right hand sides of model [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF]. The techniques in [START_REF] Bagnerini | On the role of source terms in continuum traffic flow models[END_REF][START_REF] Colombo | On a class of hyperbolic balance laws[END_REF] can then be applied.

The Initial-Boundary Value Problem.

From the point of view of traffic flow, it is natural to consider Initial-Boundary Value problems (IBVP). We start considering the case of a road starting at x = 0 where the inflow f (t) is regulated. This leads to study the following Riemann problem with boundary

∂ t u + ∂ x f (u) = 0 t ≥ 0 , x ≥ 0 u(0, x) = ū x ≥ 0 , (ρv)(t, 0) = f t ≥ 0 . (13)
We denote the maximum possible flow along the considered road by F = R f V f . When considering model ( 2), we assume that, besides (3), also

1 - Q + RV • Q + Q -1 < 1 (14)
holds. Condition [START_REF] Colombo | Road network with phase transitions[END_REF] ensures that sup Ω f ∪Ωc λ 1 < 0, hence all waves of the first family are exiting the domain x ≥ 0, t ≥ 0 and the definition of solution introduced in [START_REF] Goodman | Initial Boundary Value Problems for Hyperbolic Systems of Conservation Laws[END_REF], see also [1, Definition NC], applies: the boundary data f is attained in the sense that lim The thresholds f min and f max are showed in Figure 4 and are given explicitly in [START_REF] Colombo | Global well posedness of a traffic flow model with phase transitions[END_REF], where the proof of Proposition 1 is detailed. Note that, as remarked in [START_REF] Colombo | Global well posedness of a traffic flow model with phase transitions[END_REF], the incoming flow f can be slightly greater than the flow ρ v(ρ) present on the road. This is not allowed by the LWR model. Once the Riemann solver is available, well posedness for the Initial-Boundary Value Problem can be proved as in [START_REF] Colombo | Global well posedness of a traffic flow model with phase transitions[END_REF], for all initial and boundary data with bounded total variation. Remark that due to the presence of phase boundaries, the number of waves entering the domain (t, x) ∈ R + × R + can not be a priori established.

Q - Q + Q 0 Ω f Ω c R v = V c (ρ + f , q + f ) ρ ρv ρ 0 Ω f Ω c R (ρ, q) f min f max
We consider now the somewhat symmetric case of a road whose outflow at x = 0 is regulated. At the level of Riemann problem, this can be modeled by

∂ t u + ∂ x f (u) = 0 t ≥ 0 , x ≤ 0 u(0, x) = ū x ≤ 0 , (ρv)(t, 0) ≤ f t ≥ 0 . (15)
Proposition 2. With reference to [START_REF] Colombo | Global well posedness of a traffic flow model with phase transitions[END_REF], conditions (3) and [START_REF] Colombo | Road network with phase transitions[END_REF] imply that for all ū ∈ Ω f ∪ Ω c , and for all possible flows f ∈ [0, F ] the Riemann problem ( 15) admits a solution in the sense of [START_REF] Amadori | Continuous dependence for 2 × 2 conservation laws with boundary[END_REF]Definition NC]. More precisely, there exists a unique state ũ ∈ Ω f ∪ Ω c such that the flow at ũ is less than or equal to f and the standard solution to the Riemann problem for model [START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] with data u L = ū and u R = ũ consists only of waves having negative speed. If the flow at ū is less than or equal to f , then ũ = ū, otherwise:

1. If ū ∈ Ω f , then ũ is in Ω c .
The solution consists of a phase transition possibly followed by a 1-wave in the congested phase. 2. If ū ∈ Ω c , then the solution consists of a simple 1-wave, ũ being the intersection of the 1-Lax curve through ū and the line ρv = f .

A proof is given in [START_REF] Goatin | Modeling a bottleneck by the Aw-Rascle model with phase transitions[END_REF].

4. Road network with phase transitions. A road network is a couple (I, J ), where I is a finite collection of unidirectional roads and J is the set of junctions. Each road is modelled by real intervals I i , i = 1, . . . , N , while each junction J consists of two sets Inc(J) ⊂ {1, . . . , N } and Out(J) ⊂ {1, . . . , N } corresponding to incoming and outgoing roads of J. Given a junction J, a Riemann problem at J is a Cauchy problem with initial data constant on each incoming and outgoing road.

As for classical Riemann problems on a real line, we look for self-similar, centered solutions, which are the building blocks to construct solutions to Cauchy problems. We are now ready to introduce the key concept of Riemann solver at J. Definition 4.1. Consider a junction J and assume for simplicity Inc(J) = {1, . . . , n}, Out(J) = {n + 1, . . . , n + m}. A Riemann solver R J is a function

R J : (Ω f ∪ Ω c ) n+m -→ (Ω f ∪ Ω c ) n+m (u 1,0 , . . . , u n+m,0 ) -→ (û 1 , . . . , ûn+m ) satisfying the following 1. n i=1 f 1 (û i ) = n+m j=n+1 f 1 (û j ), where f 1 is the first component of f ; 2. for every i ∈ {1, . . . , n}, the classical Riemann problem      u t + f (u) x = 0, x ∈ R, t > 0, u(0, x) = u i,0 , if x < 0, ûi , if x > 0,
is solved by waves with negative speed; 3. for every j ∈ {n + 1, . . . , n + m}, the classical Riemann problem

     u t + f (u) x = 0, x ∈ R, t > 0, u(0, x) = ûj , if x < 0, u j,0 , if x > 0,
is solved by waves with positive speed.

To effectively describe a solution to Riemann problems at J, a Riemann solver needs to satisfy the following consistency condition: Definition 4.2. We say that a Riemann solver R J satisfies the consistency condition if R J (R J (u 1,0 , . . . , u n+m,0 )) = R J (u 1,0 , . . . , u n+m,0 )

for every (u 1,0 , . . . ,

u n+m,0 ) ∈ (Ω f ∪ Ω c ) n+m .
In what follows we will assume [START_REF] Colombo | Road network with phase transitions[END_REF], in order to have in the congested phase 1-waves always moving with negative speed. 4.1. Incoming roads: attainable values at the junction. To respect condition 2 of Definition 4.1 only waves with negative speed can be produced on incoming roads. Thus we determine all states which can be connected to an initial state (to the right) by waves with negative speed. In particular, we determine the maximum flux γ max i that can be reached from an initial datum u i,0 = (ρ i,0 , q i,0 ) by means of waves with negative speed only.

We start describing the sets of fluxes corresponding to states that can be connected to u i,0 on the right using non positive waves only. We use the notations introduced in Section 2.1, Cases (B)-(D), where we set u i,0 = u l . Moreover, we introduce the velocities V 1 and V 2 defined as follows:

•

V 1 := v f (ρ 1 )
, where ρ 1 ∈ Ω f is the smaller root of the equation

ρ 1 v f (ρ 1 ) = R c V c ; • V 2 := v f (ρ 2 ), where ρ 2 ∈ Ω f is the smaller root of the equation 1 - ρ 2 R Q + Q --Q R ρ 2 = ρ 2 V 1 - ρ 2 R .
We refer the reader to Figure 5 for help in understanding notations. The sets of reachable fluxes are then given by ρv ρ 0

Ω f Ω c R R c ρ i,0 R c V c v = V 1 ρ 1 ρv ρ 0 Ω f Ω c R R c ρ i,0 R c V c v = V 1 v = V 2 ρ 2 ρv ρ 0 Ω f Ω c R ρ i,0 ρ c V c ρ c v = V 2 ρv ρ 0 Ω f Ω c R ρ i,0 ρ c V c ρ c ρ m
O i =        [0, ρ i,0 v f (ρ i,0 )] if u i,0 ∈ Ω f , v f (ρ i,0 ) ≥ V 1 , [0, R c V c ] ∪ {ρ i,0 v f (ρ i,0 )} if u i,0 ∈ Ω f , V 2 ≤ v f (ρ i,0 ) ≤ V 1 (Case (D), Sec. 2.1) , [0, ρ c V c ] ∪ {ρ i,0 v f (ρ i,0 )} if u i,0 ∈ Ω f , v f (ρ i,0 ) ≤ V 2 (Case (C), Sec. 2.1) , [0, ρ c V c ] ∪ {ρ m v f (ρ m )} if u i,0 ∈ Ω c (Case (B), Sec. 2.1) , (16)
for i = 1, . . . , n. We observe that the sets O i are non convex. In order to have continuous dependence of solutions, we have to get convexity removing the metastable states from the attainable sets. This choice is consistent with the idea that such states should appear in a transient situation, which should not happen at a junction. Hence we define the corresponding maximum fluxes as follows:

γ max i =        ρ i,0 v f (ρ i,0 ) if u i,0 ∈ Ω f , v f (ρ i,0 ) ≥ V 1 , R c V c if u i,0 ∈ Ω f , V 2 ≤ v f (ρ i,0 ) ≤ V 1 (Case (D), Sec. 2.1) , ρ c V c if u i,0 ∈ Ω f , v f (ρ i,0 ) ≤ V 2 (Case (C), Sec. 2.1) , ρ c V c if u i,0 ∈ Ω c (Case (B), Sec. 2.1) . ( 17
)
Proposition 3. Given an initial datum u i,0 on an incoming road and γ ∈ [0, γ max i ], there exists a unique ûi such that the Riemann problem (u i,0 ûi ) is solved by waves with negative speed and f 1 (û i ) = γ. 4.2. Outgoing roads: maximal flux at the junction. To respect condition 3 of Definition 4.1 only waves with positive speed can be produced on outgoing roads. Thus we determine all states, and the corresponding set of fluxes, which can be connected to an initial state u j,0 (to the left) using waves with positive speed. We introduce the fluxes F and f max defined as follows (see Figure 6):

ρv ρ 0 Ω f Ω c R R f ρ i,0 F ρv ρ 0 Ω f Ω c R ρ i,0 f max Figure 6. Notations used in the definition of O j , j = 3, 4. • F = R f v f (R f ) = max ρ∈Ω f ρv f (ρ) > max (ρ,q)∈Ωc ρv c (ρ, q)
is the maximal flux supported by the road; • for u j,0 ∈ Ω c , f max = f max (u j,0 ) = ρ max v c (ρ max , q max ), where ρ max is the bigger root of the equation

1 - ρ max R Q + Q + -Q R ρ max = ρ max v c (ρ j,0 , q j,0 ) , and
q max = Q + ρ max (Q + -Q)/R.
The sets of reachable fluxes are given by

O j = [0, F ] if u j,0 ∈ Ω f , [0, f max ] if u j,0 ∈ Ω c , (18)
for j = n + 1, . . . , n + m. Since the sets O j are convex, the corresponding maximum fluxes are defined accordingly:

γ max j = F if u j,0 ∈ Ω f , f max if u j,0 ∈ Ω c . (19)
Proposition 4. Given an initial datum u j,0 on an outgoing road and γ ∈ [0, γ max j ], there exists a unique ûj ∈ O j such that the Riemann problem (û j , u j,0 ) is solved by waves with positive speed and f 1 (û j ) = γ. 5. Riemann solvers at junctions. In this section we describe two Riemann solvers introduced in [START_REF] Coclite | Traffic flow on a road network[END_REF][START_REF] D'apice | Packet flow on telecommunication networks[END_REF] for a LWR model, which can be used also for the phase transition model (2).

Riemann Solver R 1

J . We define a Riemann solver similar to that introduced in [START_REF] Coclite | Traffic flow on a road network[END_REF] for vehicular traffic.

First, we need to define a suitable set of matrices. Consider the set

A :=    A = {a ji } i=1,...,n, j=n+1,...,n+m : 0 < a ji < 1 ∀i, j, n+m j=n+1 a ji = 1 ∀i    . (20)
Let {e 1 , . . . , e n } be the canonical basis of R n . For every i = 1, . . . , n, we denote H i = {e i } ⊥ . If A ∈ A, then we write, for every j = n + 1, . . . , n + m, a j = (a j1 , . . . , a jn ) ∈ R n and H j = {a j } ⊥ . Let K be the set of indices k

= (k 1 , ..., k ℓ ), 1 ≤ ℓ ≤ n -1, such that 0 ≤ k 1 < k 2 < • • • < k ℓ ≤ n + m and for every k ∈ K define H k = ℓ h=1 H k h .
Writing 1 = (1, . . . , 1) ∈ R n and following [START_REF] Coclite | Traffic flow on a road network[END_REF] we define the set

N := A ∈ A : 1 / ∈ H ⊥ k for every k ∈ K . ( 21
)
Notice that, if n > m, then N = ∅. The matrices of N will give rise to a unique solution to Riemann problems at J. 1. Fix a matrix A ∈ N and consider the closed, convex and not empty set

Λ =    (γ 1 , • • • , γ n ) ∈ n i=1 [0, γ max i ] : A • (γ 1 , • • • , γ n ) T ∈ n+m j=n+1 [0, γ max j ]    . (22)
2. Find the point (γ 1 , . . . , γn ) ∈ Λ which maximizes the function

E(γ 1 , . . . , γ n ) = γ 1 + • • • + γ n , (23)
and define (γ n+1 , . . . , γn+m ) T := A • (γ 1 , . . . , γn ) T . Since A ∈ N, the point (γ 1 , . . . , γn ) is uniquely defined. 3. For every i ∈ {1, . . . , n}, set ûi either to u i,0 if f 1 (u i,0 ) = γi , or to the solution to f 1 (u) = γi given by Proposition 3. For every j ∈ {n + 1, . . . , n + m}, set ûj either to u j,0 if f 1 (u j,0 ) = γj , or to the solution to f 1 (u) = γj given by Proposition 4. Finally, set R 1 J (u 1,0 , . . . , u n+m,0 ) = (û 1 , . . . , ûn+m ) .

It is easy to verify that R 1 J satisfies the consistency condition (CC).

5.2. Riemann Solver R 2 J . In this subsection we define a Riemann solver similar to that introduced in [START_REF] D'apice | Packet flow on telecommunication networks[END_REF]. 3. Denote with (γ 1 , . . . , γn ) the orthogonal projection on the convex set Inc of the point (Γθ 1 , . . . , Γθ n ) and with (γ n+1 , . . . , γn+m ) the orthogonal projection on the convex set Out of the point (Γθ n+1 , . . . , Γθ n+m ). 4. For every i ∈ {1, . . . , n}, set ûi either to u i,0 if f 1 (u i,0 ) = γi , or to the solution to f 1 (u) = γi given by Proposition 3. For every j ∈ {n + 1, . . . , n + m}, set ûj either to u j,0 if f 1 (u j,0 ) = γj , or to the solution to f 1 (u) = γj given by Proposition 4. Finally, set R 2 J (u 1,0 , . . . , u n+m,0 ) = (û 1 , . . . , ûn+m ) .

Also in this case it is easy to verify that R 2 J satisfies the consistency condition (CC).

Figure 1 .

 1 Figure 1. Left: standard flow for the LWR model. Right: experimental data, taken from[START_REF] Kerner | Phase transitions in traffic flow[END_REF]; here q denotes the flux ρv.

) Definition 3 . 1 .

 31 Fix M > 0 and X as above. A map S : R + × D → D is an M -Riemann Semigroup (M -RS) if the following holds: (RS1) D ⊇ {u ∈ X : TV(u) ≤ M }; (RS2) S 0 = Id and S t1 • S t2 = S t1+t2 ; (RS3) there exists an L = L(M ) such that for t 1 , t 2 in R + and u 1 , u 2 in D,

Theorem 3 . 2 .

 32 For any positive M , the system (2) generates an M -RS S : R + × D → D. Moreover (CP1) for all u 0 ∈ D, the orbit t → S t u 0 is a weak entropic solution to[START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF] with initial datum u 0 ; (CP2) any two M -RS coincide up to the domain; (CP3) the solutions yielded by S can be characterized as viscosity solutions, in the sense of[START_REF] Bressan | Hyperbolic Systems of Conservation Laws[END_REF] Theorem 9.2].(CP4) D ⊆ u ∈ X : TV(u) ≤ M for a positive M depending only on M .

Proposition 1 . 1 .

 11 x→0+ ρ(t, x) • v(t, x) = f for a.e. t ≥ 0 . With reference to problem[START_REF] Colombo | On a class of hyperbolic balance laws[END_REF], if (3), (14) hold for[START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF], then for all ū ∈ Ω f ∪ Ω c , there exists a threshold f max = f max (ū) such that for all f ∈ [0, f max ] the Riemann problem for (13) admits a solution in the sense of[START_REF] Amadori | Continuous dependence for 2 × 2 conservation laws with boundary[END_REF] Definition NC]. More precisely, there exists a unique state ũ ∈ Ω f ∪ Ω c such that the flow at ũ is f and the solution to the standard Riemann problem for (2) with data u L = ũ and u R = ū consists only of waves having positive speed. If ū ∈ Ω f , then f max = F and ũ is in Ω f . The solution consists of a 2-wave in the free phase. 2. If ū ∈ Ω c , then there exists a f min = f min (ū) such that:(a) If f min ≤ f ≤ f max , ũ is the unique intersection between the curve ρv(u) = f and the 2-wave through ū. The solution consists of a simple 2-wave. (b) If f < f min , then ũ is the unique state in Ω f such that ρv f (ρ) = f . The solution consists of a phase boundary and a 2-wave.

  q

Figure 4 .

 4 Figure 4. Notation used in Proposition 1

Figure 5 .

 5 Figure 5. Notations used in the definition of O i , i = 1, 2.

i=1 θ i = n+m j=n+1 θ j = 1 . 2 .

 12 First let us defineΘ = θ = (θ 1 , . . . , θ n+m ) ∈ R n+m : θ 1 > 0, • • • , θ n+m > 0,n possible through-flow at the crossing is Γ = min {Γ inc , Γ out } . Introduce the closed, convex and not empty sets Inc = (γ 1 , . . . , γ n ) ∈ n+1 , . . . , γ n+m ) ∈

E-mail address: goatin@univ-tln.fr