Enrico Lovisari

C Canudas De Wit
email: carlos.canudas-de-wit@gipsa-lab.grenoble-inp.fr

A Y Kibangou

Density/Flow reconstruction via heterogeneous sources

Keywords: Road Transportation systems, Dynamical flow network, Density reconstruction, Floating Car Data, Optimal Sensor Placement

published or not. The documents may come

Introduction

The last decades have witnessed a considerable increase of traffic volumes, especially due to urbanisation in big metropolis, which was not matched by a comparable extension of road infrastructures. As a consequence, crucial freeways, highways and arterial roads have been steered to a state of near saturation, and experience on daily basis periods of congested traffic [START_REF] Papageorgiou | ITS and traffic management[END_REF]. In turn, congestion causes increased travel times and stop-and-go phenomena, leading to decreased safety, economical losses, and environmental and psychological hazards in terms of pollution and road rage [START_REF] Bilbao-Ubillos | The costs of urban congestion: estimation of welfare losses arising from congestion on cross-town link roads[END_REF]. Increasing road capacity by extending road infrastructures, such as construction of new arterial roads, has been the standard way to cope with congestion problems, but it is infeasible when existing roads lie on built-in areas. Intelligent Transportation Systems (ITSs), on the contrary, are expected to provide robust techniques for real-time monitoring, prediction and actuation of traffic networks, and to better integrate with road and rail public transportation, by leveraging recent technological and theoretical advancements in distributed computation and communication.

(Please place Figure 1 about here.)

The present paper investigates the two interconnected problems of sensor network design and estimation in traffic networks. As illustrated in Fig. 1, we propose three procedures to address the problems of Optimal Sensor Placement, and of Traffic Reconstruction, via the two modules of Fundamental Diagram calibration and Density and Flow Estimation.

Standard devices to obtain information on the state of the network are fixed sensors such as induction loops and magnetometers. Placed over a section of road, they provide rich information on the vehicles that cross such a section over a pre-fixed period of time: 1) their number, or flow, 2) their average speed, and 3) their average density, or more precisely their occupancy (see Section 2). Current technology allows for very precise measurements, with relative errors of measured quantities against ground truth often being below 1 ∼ 2%. However, deployment and maintenance of a sensing network requires considerable investment and manpower, and consequently sensing networks are usually designed to be as sparse as possible. In this paper we formulate and provide a solution to the Optimal Sensor Placement problem, that is, positioning sensors on the cells of the network given partial information on the system and in such a way to trade off between performance and cost. To this aim, we assume that traffic managers know the splitting ratios of the network, namely, the turning percentages at each road intersection, an information that can be obtained via standard surveys using optical count or radar sensors.

Optimal Sensor Placement is a ubiquitous problem that has received a high degree of attention in several communities due to its importance for network design. In Transportation Systems, it is of interest both in the dual-problem of best placement of hubs for cost-efficient transportation of goods and people [START_REF] Shahabi | Robust hub network design problem[END_REF] and Origin-Destination coverage [START_REF] Ehlert | The optimisation of traffic count locations in road networks[END_REF][START_REF] Hu | A generalized sensor location model for the estimation of network origindestination matrices[END_REF][START_REF] Antoniou | Towards a generic benchmarking platform for origindestination flows estimation/updating algorithms: Design, demonstration and validation[END_REF]. In these works, and dif-ferently from the present paper, the problem is cast as a mixed integer problem which corresponds to determining the minimal set of locations from which the flows on the whole network can be determined, and sensor measurements are assumed to be perfect. Finally, recent contributions focused on Origin-Destination flows and travel time estimation from an information theoretical point of view and under several possible cost functions [START_REF] Zhou | An information-theoretic sensor location model for traffic origin-destination demand estimation applications[END_REF][START_REF] Xing | Designing heterogeneous sensor networks for estimating and predicting path travel time dynamics: An information-theoretic modeling approach[END_REF], and on probabilistic approaches to the sensor placement problem that take into account sensors failures and other random events [START_REF] Fei | Vehicular network sensor placement optimization under uncertainty[END_REF]; [START_REF] Danczyk | A probabilistic optimization model for allocating freeway sensors[END_REF].

Once the sensing network has been designed an implemented, flow and density measurements are available. Sensors are not, however, the only source of information that we exploit. In fact, the recent spread of wireless devices allows sensing and communication capabilities unforeseeable up to few years ago. Limiting the attention to traffic applications, vehicles equipped with positioning devices (such as GPS) and able to communicate with an ITS monitoring system can act as a probes in the traffic and provide Floating Car Data (FCD), namely, information on the vehicles' positions and speeds. The collected data can be used to estimate the speed in the network, thus offering a second source of information. Due to privacy reasons, single vehicles traces are usually not directly used, but rather aggregated as average speed of vehicles in segments of road. Advanced methodologies ensure fine spatial partitions of the network, with segments as short as 250 meters (INRIX, 2014). Compared to fixed sensors, a service based on this technology can only make use of information coming from her customers, which are a fraction of the total vehicles on the road (the penetration rate of the system). This implies that speed measurements are less precise and flow measurement are unavailable. On the other side, since it exploits existing communication systems it is relatively inexpensive and, more important, already covers all major traffic networks.

Fixed sensor measurements and Floating Car Data provide rich information that we employ to address the problem of estimating road usage in terms of density and flow of vehicles in a traffic network. The latter are commonly considered a good representation of the state of the system, providing more information than average speed alone. In particular, they are of crucial importance for 1) forecasting travel time and traffic evolution, along with historical data; 2) informing in real-time drivers about the state of the network through navigation systems; 3) providing public authorities with statistical data to monitor the state of the network and predict dangerous scenarios; 4) computing and actuating control actions through traffic lights, ramp metering and speed limits, or, in the future, lane change and semi-autonomous routing and navigation [START_REF] Papageorgiou | Review of road traffic control strategies[END_REF][START_REF] Papageorgiou | Alinea: A local feedback control law for on-ramp metering[END_REF][START_REF] Pisarski | Optimal balancing of road traffic density distributions for the cell transmission model[END_REF][START_REF] Como | Throughput optimal distributed routing in dynamical flow networks[END_REF].

Traffic models for analysis and control synthesis date back to the first half of the 20th century. The most celebrated macroscopic model is the PDE based Lighthill-Whitham and Richards (LWR) model [START_REF] Lighthill | On kinematic waves. ii. a theory of traffic flow on long crowded roads[END_REF], which, based on fluid kinematics, is able to reproduce crucial phenomena such as traffic shock waves. Discretization of the LWR-PDE is not straightforward but stable numerical schemes have been proposed, the most well known being the Cell Transmission Model (CTM) [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF][START_REF] Daganzo | The cell transmission model, part II: network traffic[END_REF]. Huge efforts have been put in the last 15 years to calibrate the CTM [START_REF] Muñoz | A piecewise-linearized cell transmission model and parameter calibration methodology[END_REF] and to unveil its system-theoretical properties [START_REF] Morbidi | A new robust approach for highway traffic density estimation[END_REF]. Fusion of flow, density and speed measurements has also been addressed, even though mostly considering single vehicles traces. Approaches range from signal processing techniques such as the generalized Treiber-Helbing filter (van Lint and Hoogendoorn, 2010), nonlinear versions of the Kalman filter in the context of Lagrangian sensing [START_REF] Work | A traffic model for velocity data assimilation[END_REF], stochastic versions of the three-detector model [START_REF] Deng | Traffic state estimation and uncertainty quantification based on heterogeneous data sources: A three detector approach[END_REF], and Extended Kalman Filter [START_REF] Nantes | Realtime traffic state estimation in urban corridors from heterogeneous data[END_REF]. Alternative approaches do not rely on discretization of the LWR-PDE model and allow to cast problems of estimation and control as convex problems [START_REF] Li | Optimal control of scalar conservation laws using linear/quadratic programming: Application to transportation networks[END_REF]. Finally, an extremely recent contribution [START_REF] Wright | Fusing loop and gps probe measurements to estimate freeway density[END_REF] addresses the density estimation problem from a point of view very similar to ours and making use of particle filtering strategies.

In the first part of this paper we propose an algorithm that aims to reconstruct the traffic density and flow by fusing fixed sensors measurements and Floating Car Data. We employ a macroscopic model, partitioning the network in cells and assigning to each cell a density of vehicles. The latter evolves dynamically according to a first order mass-conservation law, similarly to the Cell Transmission Model. We inherit from the CTM the assumption that the inflow in a cell is a fixed linear combination of the outflows of the preceding cells. Differently from CTM, however, inflows and outflows in all the cells are estimated on the basis of the available flow measurements only. In addition, using the concept of Fundamental Diagram and the speed measurements, we compute an instantaneous (namely, only based on the latest available measurements) pseudo-measurement of the density. These quantities are then the inputs for the density observer. As an intermediate step, we propose a gradient descent method to calibrate the Fundamental Diagram which can be run offline and only requires fixed sensor measurements.

In the second part of the paper we address the problem of Optimal Sensor Placement, namely, the problem of finding the best location where to physically place sensors. This is based on trading off between two contrasting objectives: the first, to maximize the performance of state reconstruction; the second, to minimize the total economic cost of the network. The performance of the state reconstruction is usually related to the ability to properly estimate the density of vehicles in the roads. Unfortunately, nonlinearity and complexity of traffic systems make it hard to evaluate the performance of nonlinear observers. In order to simplify the setting, we consider the related problem of reconstruction of flows in a static setting. In particular, we consider as performance metric the error covariance of the Best Linear Unbiased Estimator of the cumulative flows in the network over a long period of time under a simple linear measurement model. The resulting cost has an analytical formulation and can be easily studied. The rationale of such a choice is the following: In the proposed den-sity estimation algorithm availability of speed measurements via FCD makes the density estimation directly dependent on the quality of flows estimation. As such, sensor placement for optimal flow reconstruction is also expected to yield good performance for density estimation. The Optimal Sensor Placement problem can be then seen as trading off between the performance of such a flow estimator, and a cost that depends on the dimension of the sensing network. Since this is a combinatorial problem, we relax it using a method that we call Virtual Variance algorithm, based on the idea to associate to each sensor a virtual variance which is large when the sensor is not relevant for good reconstruction of the flow vector. The only input that the algorithm needs is the matrix of splitting ratios, that prescribes how vehicles split at each junction, and the nominal noise variance of each sensor. Furthermore, we discuss in detail two extensions of the proposed algorithm dealing with important scenarios. In the first, Optimal Sensor Placement with geographical constraints, we address the scenario in which sensors cannot be placed in a subset of cells of the network. In the second, Optimal Sensor Placement with Number of Sensors constraints, we deal with the case in which the maximum number of sensors is pre-specified, for example due to budget limitations.

To summarize, the contributions of this paper are the following: 1) we formulate the problem and we design an easily implementable approach to data fusion of fixed sensors measurements and Floating Car Data; 2) we propose a gradient descent calibration algorithm of the underlying macroscopic model; 3) we formulate the problem of Optimal Sensor Placement in terms of positions of sensors in a network when sensors are noisy and we provide an approximate solution using the concept of Virtual Variances; 4) we show the prowess of the devised Optimal Placement procedure on a regular grid, for which we offer a comparison between the solution found with our approach and the true optimal placement, found by exhaustive search; 5) we illustrate the performance of the Optimal Placement procedure and of the Reconstruction algorithm through extensive numerical experiments in the real-world scenario of Grenoble Traffic Lab (GTL) (Canudas de Wit et al., 2015), a sensing network deployed along the freeway "Rocade Sud" in Grenoble, France, with FCD provided by INRIX, one of the most well known traffic solutions companies.

The remainder of the paper is organized as follows: after setting up the notation, Section 2 describes the model for a Road Transportation Network. Section 3 formulates the problem of flow and density reconstruction and describes the proposed nonlinear observer, while the problem of optimal sensor placement is formulated and a solution based on the heuristic Virtual Variance algorithm is presented in Section 4. Finally, Section 5 illustrates the solutions on a regular grid and on the real world scenario of the freeway Rocade Sud in Grenoble, and Section 6 draws the conclusions and provides several future research directions.

Notation

The symbol R A (R A + ) for a finite set A denotes the set of real vectors (real positive vectors) indexed by elements of A. For x ∈ R n , x ≥ 0 is meant component-wise. The identity matrix of suitable dimensions is denoted I, and |A| is the cardinality of the set A. The symbol E denotes the expected value.

A graph G is a pair (V, E) where V is called the set of nodes and E the set of edges. For an edge e ∈ E, t(e) ∈ V and h(e) ∈ V are the tail note and the head node of e, so that e = (t(e), h(e)). Denote by E + e := {j ∈ E : h(e) = t(j)} and E - e := {j ∈ E : h(j) = t(e)} the set of edges that follow or precede e, respectively. A path of length n ≥ 2 is a sequence of edges e 1 , . . . , e n such that e i+1 ∈ E + ei for all i = 1, . . . , n -1. A path of length 1 is a path made of a single link. The matrix L ∈ R V×V is a sublaplacian of G if L ej > 0 only if (e, j) ∈ E, e = j, and j L ej ≤ 0.

Table 1 presents the notation that will be used throughout the paper.

(Please place Table 1 about here.)

Road Transportation System Model

We adopt a macroscopic approach by partitioning the lanes of the roads in a traffic network in cells. Cells that lie on the same section of a road and on different lanes are said to be parallel one each other (see Figure 2) We interpret each cell as an edge e ∈ E in a graph G = (V, E), which models the whole network. In such a graph, nodes v ∈ V represent intersections between cells, for example, junctions at which roads intersect. Long roads might be partitioned into shorter segments, and are thus represented as a series of consecutive cells. Among the set of cells E, we denote by R i and R o the set of onramps and offramps, respectively. In this paper, we shall call onramp (offramp) any gate, be it a real ramp, a connector, a secondary road, etc, which lets vehicles enter into (exit from) the network.

(Please place Figure 2 about here.)

We make the following connectivity assumption, which formalizes the mild requirement that any cell can be reached from an onramp and that vehicles from any cell can exit from the network.

Assumption 1. For any cell e in E, there is at least one onramp j ∈ R i and one offramp k ∈ R o such that e is an edge of a path from j to k.

Time is discrete and slotted in intervals of duration T > 0. On each cell e ∈ E, denote by ρ e (t) the density of vehicles, in number of vehicles per km2 , at time t, and let ρ(t) = [ρ e (t)] e∈E , which we call the state of the network. The density of vehicles in a cell changes dynamically in time according to the following mass-conservation first-order model

ρ e (t + 1) = ρ e (t) + 1 e (f in e (t) -f out e (t)), ∀e ∈ E (1)
where e is the length of cell e, and f in e (t) and f out e (t) are the inflow and the outflow at cell e during the t-th time slot. As such, f in e (t) and f out e (t) represent a number of vehicles, not, as more common, a number of vehicles per unit of time, e.g., per hour.

To relate inflows and outflows we resort to the standard concept of splitting ratios. Indeed, denote by R ek ≥ 0 the fraction of vehicles that turn into cell k when they exit from cell e, which is the splitting ratio of e towards k. Clearly, R ek = 0 if e and k are not consecutive, and k R ek = 1, if e ∈ R o . Such a cell change might correspond to continuing on the same lane, lane change in the same road, or taking a different road, as exemplified in Figure 2. From this moment on, we make the following assumption.

Assumption 2. The set of splitting ratios {R ej } (e,j)∈E×E is known.

Splitting ratios establish the relation f in e (t) = j∈E R je f out j (t), for all t ≥ 0 and for any cell e ∈ R i , while f in e (t) = λ e (t) for e ∈ R i , where λ e (t) is an exogenous external inflow. We set for sake of convenience λ e (t) = 0 for e ∈ R i . By stacking inflows and outflows into vectors f in (t) and f out (t), respectively, we can rewrite the previous relation in matrix form as

f in (t) = R T f out (t) + λ(t) (2)
where the matrix R = [R ej ] e∈E,j∈E\R i is the matrix of splitting ratios. In the present paper it is assumed that the matrix of splitting ratios is fixed, predetermined, and known. Its calibration, closely related to the estimation of Origin-Destination pairs, can be performed on single-lane freeways with onramps and offramps by taking ratios of flows on main line and ramps [START_REF] Muñoz | A piecewise-linearized cell transmission model and parameter calibration methodology[END_REF]. We plan to extend this setting to networks by casting the problem as an optimization problem in future research.

On modelling of cell flows

Macroscopic models such as the CTM postulate that the flow f out e (t) that exits from cell e at time t is a deterministic function of the densities in the cells around e. In the simplest case, where e and j are two consecutive cells and E + e = {j} and E - j = {e}, then in the CTM

f out e (t) = min{d e (ρ e (t)), s j (ρ j (t))}
where d e (ρ e (t)) and s j (ρ j (t)) are the demand of cell e and the supply of cell j, and represent the maximum outflow from e and the maximum inflow into cell j, respectively. The resulting system is a Godunov scheme for discretization of the LWR-PDE model, and can be extended to the network case in various ways [START_REF] Daganzo | The cell transmission model, part II: network traffic[END_REF][START_REF] Lovisari | Data fusion algorithms for density reconstruction in transportation networks[END_REF][START_REF] Coogan | Dynamical properties of a compartmental model for traffic networks[END_REF]. While these models reproduce important phenomena that must be taken into account when modelling and controlling traffic networks, such as the movement of shockwaves, each of them can only partially represent traffic dynamics in networks.

For this reason, we will avoid to explicitly model the relation between flows and density, and we will limit to the standard

f out e = ρ e v e , ∀e ∈ E (3)
namely that the volume that exits from a cell in a period is proportional to the density of the vehicles in the cell, and to their speed. Again, we leave unmodelled the relation between these two quantities, because, as we will make clear in the following, we assume to have a direct measurement of the average speed in each cell (in the form of Floating Car Data). In conclusion, we consider from now on the dynamics of the real system to be dictated by Eqs. ( 1)-( 3), where v e , e ∈ E, is an unmodelled quantity which depends on the local state of the network.

While in the model of our network we do not use an explicit relation between flows and densities, we shall use it for data fusion and estimation purposes. In particular, we write ϕ e = ϕ e (ρ e ), ∀e ∈ E (4)

where ϕ e is the flow of vehicles at the sensor locations, which we shall always assume to be at the end of the cell. The graph of the function ϕ e (•), which is the Fundamental Diagram on cell e, is a concave function with ϕ(0) = ϕ(ρ jam ) = 0, where ρ jam is the jam density, at which vehicles are too close one each other to move. The value of ρ jam varies from 150 to 300 vehicles per km, and is approximately related to average vehicle length by the formula ρ jam = 1/ ave , where ave is measured in km.

In this paper, and consistently with (1), we make the steady assumption that ϕ e measures the flow of vehicles through a certain section of road per sampling time T , as opposed to standard choices such as flow per minute or per hour. This is done without loss of generality, as the two functions only only by a multiplicative constant, and solely for sake of notation and to simplify the discussion in the numerical section, where all relevant variables will be normalized to the Grenoble Traffic Lab sampling time T = 15 seconds.

Available measurements

In this paper we consider two heterogeneous sources of information: flow and density measurements from sensors, and Floating Car Data.

Flow and density measurements

Standard measurement devices for traffic are loop detectors or magnetometers, radar traffic detectors, or video detection systems. They are positioned at fixed and predefined positions in the network, monitor a section of road, and are able to detect and assign a timestamp to the event "a vehicle crossed the section". Information is then aggregated in time slots. For sake of simplicity and without loss of generality, we assume that such time slots correspond to the time discretization of the system (1). As such, measurements of the following two quantities are available at every time slot of duration T :

• Flow of vehicles ϕ e (t), by counting the number of vehicles crossing the section during the t-th time slot;

• Density of vehicles ρ e (t) over the section. The quantity that is actually measured by the aforementioned devices is the occupancy of vehicles o e (t), defined as the percentage of time any vehicle was standing over the section in the t-th time slot. Since occupancy is approximately related to the density of vehicles by ρ e ≈ oe 100 ave , where ave is the average length of a vehicle in km, we assume from now on that density can be directly measured.

Measurements are unavoidably noisy, with sources of noise ranging from temporary inability to detect changes of the magnetic field, too fast or too slow vehicles, blurred videos, etc. We adopt a simple additive noise model

ϕ m e (t) = ϕ e (t) + ω ϕ e (t), e ∈ E m ρ m e (t) = ρ e (t) + ω ρ e (t), e ∈ E m (5)
where ϕ m e (t) and ρ m e (t) are flow and density measurements at time t, and ω ϕ e (t) and ω ρ e (t) are measurement errors whose stochastic properties depend on the performance of the sensor as well as on road and weather conditions, and E m ⊆ E is the set of cells equipped with sensors. Due to installation and maintenance costs, usually

|E m | << |E|.
Remark 1. It should be mentioned that measurement systems based on magnetometers can be used to measure the average speed of vehicles crossing the section they monitor, in addition to flow and density, simply by deploying them in pairs monitoring sections that are at a fixed and known distance. The two consecutive instants at which the same vehicle crosses the two sections provides then a measurement of its speed. This type of installation is however not standard and more expensive due to additional hardware and software. For this reason, and to show that our approach does not need this additional information, we shall assume that no speed measurement from static sensors is available.

Speed measurements

As already mentioned, recent technological advancements provide public and private companies with speed measurements in the form of Floating Car Data (FCD). In this paper we assume that the traffic network is partitioned in segments, let S be the set of all segments, and that aggregate speed measurements are available on each segment s ∈ S. In particular, such an aggregate information is the average speed of the vehicles on each segment, where averaging is required for privacy reasons. While some classes of vehicles, such as taxi and buses, can indeed be traced, we do not employ this additional information, and leave the possibility to use it as future research. Notice that segment partitioning is, in general, different from cell partitioning, since the former depends on communication and sensing capabilities, while the latter on modeling choices. In this paper we assume that FCD segments are coarser than cells, so that each segment usually consists of several cells, as in Figure 3.

(Please place Figure 3 about here.) Floating Car Data are less expensive and require less maintenance effort with respect to fixed sensors since they exploit existing communication architectures. For the same reason, they are already accessible almost everywhere once a data collecting mechanism is deployed. Despite such advantages, FCD also have drawbacks. Aside from the already mentioned possibly low penetration rate, information provided via FCD is usually averaged over a relatively long period of time. As an example, within the Grenoble Traffic Lab fixed sensors yield flow and density (and speed) measurements every T = 15 seconds. The FCD provided by INRIX are instead aggregated per minute, with standard practice ranging between 5 and 10 minutes. A comparison between sensors speed measurements from the GTL and FCD is provided in Figure 4, in which we plot FCD measurement and average of GTL measurements for slow and fast lanes at location Taillat on the Rocade Sud in Grenoble from 07:00:00 to 11:00:00 on April 24th, 2014 (see Paragraph 5.1.2 for details), both sampled at a rate of one measurement every 15 seconds (left panel) and averaged with a 5 minutes rate (right panel). The left plots clearly illustrates the difference between the high measurement rate of sensors and the effect of temporal averaging of FCD, the latter yielding a smoother signal. Such a phenomenon is much less apparent in the right plot, where the sampling time is the same. As a further comparison, we plot both speed measurements for all the cells on the main line at 09:00:00 on Figure 5, showing again that on average the measurements are in agreement. In fact, the average difference between sensor and FCD speed measurement is m e = -5.8, with standard deviation s e = 13.5, whereas the average speed on the Rocade is 82 km/h (according to FCD data).

(Please place Figure 4 about here.) (Please place Figure 5 about here.)

We provide now a more formal description of the FCD speed measurements. We partition time into intervals of duration N time slots (N T seconds), that is, into FCD acquisition intervals [0, N -1], [N, 2N -1], . . . . At the beginning of each interval, namely at time t = N, 2N, 3N, . . . , new aggregate data are available in the form of the average speed during the whole last interval and over each segment (which, recall, can consist of several cells). For a segment s ∈ S, we denote such aggregate datum by v m s (k), where k refers to the k-th FCD acquisition interval. These speed measurements are kept for the whole subsequent FCD acquisition interval (namely, for N time slots), when they are again updated. During the first interval, namely, for t = [0, N -1], no data are available, and hence the speed measurement can be set to an arbitrary value. In this paper we set it to the freeflow speed, which is a reasonable choice if the first FCD acquisition interval falls during night time, when traffic is very low. In formulae, the speed measurement on cell e at time t is

v m e (t) = v ff e , t ∈ [0, N -1] v m s(e) (k), t ∈ [kN, (k + 1)N -1] (6)
where

• v ff e > 0 is the freeflow speed on cell e, namely, the speed of vehicles in low density regime.

• v m

s(e) (k), the aggregate datum, is given by

v m s(e) (k) = 1 N |s(e)| j∈s(e),τ ∈It v j (τ ) + ω m s(e) (k)
where s(e) denotes the segment of which e is one of the cells (see Figure 3), ω m s(e) (k) is a measurement error whose stochastic properties depend on the performance of the sensor as well as on road and weather conditions, and 6 about here.)

I t = {τ : t N -1 ≤ τ N < t N }. (Please place Figure

Flow and Density Reconstruction

The first problem we address in this paper is Flow and Density Reconstruction on the basis of heterogeneous sources of information. Our aim is to build an observer for the densities of vehicles in all cells of the network given static sensor measurements and Floating Car Data.

We start by observing that Eqs. ( 1)-( 2) cannot be directly used to observe the system except for the ideal scenario in which ideal measurements of the outflows f out e (t), for all e ∈ E and for all times t ≥ 0, and of the initial conditions of the system, are available. Such a naive observer would however be very sensitive to noise, as notice that while errors in the initial conditions correspond to offsets during the evolution of the system, noises in the flow measurements are integrated by the system's dynamics, thus possibly producing unbounded and/or unrealistic results. Since real systems are never error free, Eqs. ( 1)-( 2) cannot be directly used to observe the system.

We address this difficulty by considering the following standard Luenbergerlike observer

ρ(t + 1) = ρ(t) + L -1 ( f in (t) -f out (t)) + Γ(ρ(t) -ρ(t)) (7)
where

• ρ(t) is the stacked vector of {ρ e (t)} e∈E , the density estimates on cells e ∈ E at time t;

• f in (t) and f out (t) are the stacked vectors of { f in e (t)} e∈E and { f out e (t)} e∈E , the inflow and outflow estimates on cells e ∈ E at time t. Inflow and outflow estimates are to-be-designed functions of the flow measurements ϕ m (t) = [ϕ m e (t)] e∈E m , as detailed in the next section;

• ρ(t) is the stacked vector of {ρ e (t)} e∈E , which are density pseudo-measures on cells e ∈ E at time t. Density pseudo-measures are functions of both flow measurements and speed measurements ϕ m (t), v m (t) = [v m e (t)] e∈E , as detailed in the next section;

• Γ is a tunable gain trading off between flow and density pseudo-measure;

• L is the diagonal matrix of cell lengths, with L ee = e for all e ∈ E.

We are interested in the following:

Problem 1 (Flow and Density Reconstruction using Heterogeneous Sources): Design the maps

f in = { f in e } e∈E : R E m + → R E f out = { f out e } e∈E : R E m + → R E ρ = {ρ e } e∈E : R E m + × R E + → R E
to minimize the absolute errors with respect to real flows and densities

a ρ (t, e) = |ρ e (t) -ρ e (t)| a ϕ (t, e) = f out e (t) -f out e (t) (8)

A nonlinear observer for traffic networks

In this subsection we describe the proposed solution to Problem 1. It consists in an offline calibration procedure and an online filtering step.

Offline calibration

This paragraph is devoted to providing a solution for calibrating the Fundamental Diagram. We employ a gradient descent strategy, similarly to [START_REF] Qu | On the fundamental diagram for freeway traffic: A novel calibration approach for single-regime models[END_REF], for which we do not require a CTM formulation as in [START_REF] Munoz | Methodological calibration of the cell transmission model[END_REF].

Recall that the Fundamental Diagram is the graph of the function ϕ e (•) that is the flow of vehicles at the point where sensors are placed on cell e and whose argument is the density of vehicles on cell e. As such, it can be only estimated on cells e ∈ E m , where measurements of flow and density are available. For a cell j ∈ E \ E m , we assume that the Fundamental Diagram can be estimated by extending by linear interpolation the parameters on the cells in E m that are close to j, with coefficients that depend on the mutual distance between the cells and on the type of road j belongs to.

The profile of Fundamental Diagram that we consider in the present paper is the following e ), in which the mutual influence of vehicles is small, from the high-density congested region (ρ c e , ρ jam e ], in which speed decreases with density due to interaction of close vehicles;

ϕ e (ρ) = v ff e ρ, ρ ≤ ρ c
• ρ jam , the jam density, is assumed to be known a priori; where ω e is the wave speed at section e. Clearly, our model recovers the latter with a e = 0, b e = -ω e , and c e = ω e ρ jam e . The choice of a quadratic Fundamental Diagram in congestion regime has been driven by the empirical observation, based on data on our experimental setting, that the triangular diagram tends to overestimate the flow in congestion, as it will be shown in Section 5. For this reason, we call it linear-quadratic Fundemental Diagram.

• v ff e >
An alternative appealing solution which fits our data is the inverted-λ fundamental diagram [START_REF] Hall | Empirical analysis of freeway flow-density relationships[END_REF]. However, the number of parameters to be estimated is higher in the latter case, and the resulting model is more complex as it involves hysteresis.

Since it is well known that deterministic Fundamental Diagrams are in any case only a rough approximation of the relation between flow and density, we chose the quadratic profile because it is simple to calibrate and to use.

We describe now the procedure for the calibration of the Fundamental Diagram. Let e ∈ E m and let {(ρ m e,k , ϕ m e,k )} k∈K , K = {1, . . . , K}, the set of K density and flow measurements used as learning set and obtained via the fixed sensor on cell e. For sake of notation, and since all variables refer to cell e only, let us write from now on ρ m k and ϕ m k instead of ρ m e,k and ϕ m e,k , and the same for the parameters of the Fundamental Diagram.

The proposed calibration procedure requires two steps

• Estimation of ρ c and C = v ff ρ c : the first step consists in estimating the critical density and the capacity C of the cell. We consider the standard least square estimation, which results into the following non-linear and non-convex minimization problem: given the set of measurements

{(ρ m k , ϕ m k )} k∈K , K = {1, . . . , K}, solve min (ρ c ,C) V (ρ c ,C) = 1 2 K k=1 (ϕ k -ϕ (ρ c ,C) (ρ k )) 2 s.t. 0 < ρ c < ρ jam C > 0 ϕ (ρ c ,C) (x) = C ρ c x, x ≤ ρ c C(ρ jam -x) ρ jam -ρ c , x > ρ c (9)
We aim to solve (9) by the following gradient descent with diminishing stepsize algorithm -Basic step: initialize ρ c 0 , C 0 . A reasonable choice is ρ c 0 = 20, which corresponds to the vehicles influencing one each other when the average distance among them is less than 50 meters, and C 0 = v limit e ρ c 0 , where v limit e is the speed limit on cell e normalized by the sampling time T ;

n-th step: let (ρ c n , C n ) descend along the gradient of the cost, namely

ρ c n+1 = ρ c n - δ n ∇ ρ c V (ρ c ,C) C n+1 = C n - δ n ∇ C V (ρ c ,C) with ∇ ρ c V (ρ c ,C) = k∈I FF (ρ c n ) ϕ k -ϕ (ρ c n-1 ,Cn-1) (ρ k ) C n-1 (ρ c n-1 ) 2 ρ k - k∈I C (ρ c ) ϕ k -ϕ (ρ c n-1 ,Cn-1) (ρ k ) C n-1 ρ jam -ρ k (ρ jam -ρ c n-1 ) 2 ∇ C V (ρ c ,C) = - k∈I FF (ρ c n ) ϕ k -ϕ (ρ c n-1 ,Cn-1) (ρ k ) ρ k ρ c n-1 - k∈I C (ρ c ) ϕ k -ϕ (ρ c n-1 ,Cn-1) (ρ k ) ρ jam -ρ k ρ jam -ρ c n-1 I FF (ρ c ) = {k ∈ K : 0 < ρ k ≤ ρ c } I C (ρ c ) = {k ∈ K : ρ k > ρ c } where the gradients ∇ ρ c V (ρ c ,C) and ∇ C V (ρ c ,C) are computed at the point (ρ c , C) = (ρ c n-1 , C n-1
), and δ > 0 is a fixed initial step size. Notice that if ρ c n = 0 then I FF (ρ c ) = ∅, and conversely when ρ c n = ρ jam then I C (ρ c ) = ∅, and thus the previous summations are always well defined. Nonetheless, for numerical reasons, additional care should be taken in order to avoid ρ c < 0 or ρ c > ρ jam , for example projecting at each step ρ c n+1 into [0, ρ jam ] after the gradient update.

-

Stopping criterion: stop if ρ c n C n - ρ c n-1
C n-1 < ε for some small threshold ε > 0.

• Calibration of the quadratic function in the congested region: the problem of calibrating the quadratic function for the congested region is cast into the quadratic problem: given the set of measurements

{(ρ m k , ϕ m k )} k∈K , K = {1, . . . , K}, solve min (a,b,c) k∈I C (ρ c ) (ϕ k -(aρ 2 k + bρ k + c)) 2 s.t. a (ρ c ) 2 + bρ c + c e = C a ρ jam 2 + bρ jam + c = 0 a ≥ 0 (10)
The problem ( 10) is computationally very simple and can be solved using off-the-shelf tools.

Notice that as side products of the previous procedure we can compute the freeflow speed as v ff = C/ρ c and, in case a bilinear Fundamental Diagram is also needed, the wave speed as ω = -C/(ρ jamρ c ), for each cell.

Finally, we remark that, in principle, one could add in (9) an upper bound on the capacity, for example C max = v limit ρ c . Such an additional constraint can be included in the cost, and hence in the gradient iterates, via a standard barrier function. However, such a choice might be in practice too restrictive (for example, drivers tend to exceed the speed limit). Furthermore, since the algorithm aims to fit the (density, flow) pairs, it is to be expected that the resulting capacity does not exceed the maximum measured flow, which can be verified a posteriori. A capacity upper bound is therefore a modeling choice, which for the aforementioned reasons has not been implemented in the present paper.

Online density reconstruction algorithm

We assume from now on that Fundamental Diagrams have either been calibrated or extended on the whole network, and that the matrix of splitting ratios has been pre-specified or estimated on the basis of field surveys.

We propose the following online algorithm for Density Reconstruction

• at the beginning of the t-th time slot, a centralized computation unit receives measurements {ϕ m e (t)} e∈E m ;

flow estimation: estimates the vector of outflows f out (t) by solving the following minimization problem

min f out ||(I -R T ) f out || 2 + γ e∈E m ( f out e -ϕ m e (t)) 2 s.t. f out ≥ 0 (11)
Problem ( 11) aims to a) match outflows and measurements where available, by penalizing the squared difference between the two, and b) to balance outflows according to the splitting ratios. The latter term provides the estimate of flows on cells in which no measurement is available, and is performed "as if" the network were at steady state, which is a simplifying assumption due to absence of a dynamical model for flows. The tunable parameter γ selects whether more weight is given to matching estimated outflows and measurements (high γ), or to estimate the flows as if the network were at steady state (low γ). Once Problem ( 11) is solved, the vector of estimate of the inflows is easily computed according to Eq. 2, by setting as a rough estimate of the density. This estimate is only based on the actual measurements of flow and speed, and is generally very noisy, especially when the cell is in congestion. Therefore, the algorithm does not directly uses it;

f in e (t) = j∈E - e R je
density estimation: for each cell e ∈ E, lets the density estimate evolve according to the observer equation (7).

Observer dynamic analysis

This section is dedicated to the analysis of the estimation error dynamics for the observer (7). Let e(t) = ρ(t) -ρ(t) ∈ R E be the estimation error at time t. By (1) and ( 7), the error evolution is given by

e(t + 1) = ρ(t) + L -1 (f in (t) -f out (t)) -ρ(t) -L -1 ( f in (t) -f out (t)) -Γ(ρ(t) -ρ(t)) = (I -Γ)e(t) -L -1 (R T -I)∆f out (t)) -Γ∆ρ(t) (12)
where we defined

∆f out (t) := f out (t) -f out (t), ∆ρ(t) := ρ(t) -ρ(t)
the flow and density errors. The quantities ∆f out (t) and ∆ρ(t) can be interpreted as noises, but we observe that they 1) depend on the measured quantities, and 2) do not follow a Gaussian distribution. For these reasons, a straightforward application of Kalman filter is not theoretically justified, and would not be optimal in the minimum variance sense. One could resort to employing extensions of Kalman filter based on system linearization, such as extended Kalman filter, but such strategies would be more complex than the proposed linear estimator, and their optimality would not be proved anyway. Finally, particle filters might provide better performance than both our algorithm and the extended Kalman filter, as they run in parallel a number of scenarios and iteratively select the most informative trajectories, but would be much more complex and computationally requiring than the simple linear observer proposed.

We observe now that the error equation ( 12) is essentially a low pass filter of the total noise signal n(t) = L -1 (R T -I)∆f out (t)) + Γ∆ρ(t) with stability and bandwidth dependent on the eigenvalues of Γ. Let in fact Λ = I -Γ. Stability, a strictly necessary condition, requires that |λ| < 1 for any eigenvalue λ of Λ. Within this range, eigenvalues of Λ close to 0 yield large bandwidth, hence a faster filtering but a more noisy error, and viceversa eigenvalues close to 1 imply small bandwidth, thus less influence of the total noise n(t) but slower decaying rate of e(t) to zero (in absence of n(t)).

Further analysis can be carried on under different assumptions on the noise terms ∆f out (t)) and ∆ρ(t). We shall consider for sake of simplicity the scenario in which Γ = gI, so that stability requires g ∈ (0, 2). We first present a deterministic analysis based on Lyapunov theory and the concept of ultimate boundedness [START_REF] Khalil | Nonlinear Systems[END_REF]. Assume that

||L -1 (R T -I)∆f out (t))|| ∞ ≤ M f and ||∆ρ(t)|| ∞ ≤ M ρ ,
namely, that the noise components of the error equation are bounded. Notice that such M f and M ρ are surely finite by definition of ∆f out (t) and ∆ρ. Let us consider the discrete time Lyapunov equation

(1 -g) 2 p -p = -q =⇒ p = q g(2 -g)
where q > 0, and the candidate Lyapunov function V (e) = p||e|| 2 . Tedious but straightforward computations3 yield

V (e(t + 1)) -V (e(t)) ≤ q -||e(t)|| 2 2 + 2(M f + gM ρ )|1 -g| √ N g(2 -g) ||e(t)|| 2 + N (M f + gM ρ ) 2 g(2 -g)
We can now establish uniform boundedness. Since q > 0, the right hand side is a concave parabola in ||e(t)|| 2 with roots (±1+|1-g|)(M f +gM ρ ) g(2-g)

. As a consequence, after possibly an initial transient, the error trajectory will be bounded within the sphere of e(t) such that ||e(t)|| 2 ≤ r(g) := (1+|1-g|)(M f +gM ρ ) g(2-g)

, since outside such a sphere V (t), and hence by construction ||e(t)|| 2 , cannot but decrease. A reasonable tuning choice is to minimize the radius of such a sphere, namely to minimize r(g) as a function of g. Straightforward computations yield that the minimum radius is r min = M f + M ρ , obtained for g min = 1. By this choice we make the filter as fast as possible and the ultimate boundedness radius as small as possible.

A second, probabilistic, analysis is based on a statistical description of ∆f out (t)) and ∆ρ(t). In particular, assume that ∆ρ(t) and ∆f out (t) are two jointly independent and identically distributed stochastic processes with ∆f =

E(L -1 (R T -I)∆f out (t))), σ 2 f = E||L -1 (R T -I)∆f out (t)) -∆f (t)|| 2 , ∆ρ = E∆ρ(t), σ 2 ρ = E||∆ρ(t)) -∆ρ(t)|| 2
. These values can be obtained numerically as empirical average and empirical second order moment of ∆f out (t)) and ∆ρ(t). We remark once again that these assumptions do not hold true in the framework under analysis, where the two processes are not independent and where their statistical properties depend in a complex way on the state of the system -this is one of the reasons for which directly employing a Kalman filter strategy is not justified. Nonetheless, they are instrumental for the tuning strategy presented in this section. Let m(t) and σ 2 e (t) be the expected value and the second moment of e(t), which evolve according to m(t + 1) = (1 -g)m(t) + ∆f (t) + g∆ρ(t) and σ 2 e (t + 1) = (1 -g) 2 σ 2 e (t) + σ 2 f + g 2 σ 2 ρ . Standard analysis establishes that asymptotically m(t) → m(g) = ∆ρ(t) + ∆f (t) g and σ 2 e (t) → σ2 e (g) =

σ 2 f +g 2 σ 2 ρ g(2-g) .
Tuning the gain g thus involves trading off between minimizing the expected value of the error me (g), for which the best choice is g * m = 2 (the maximum possible), and minimizing the error variance σe (g), for which the best choice can be found to be g

* σ = -1+ √ 1+4σ 2 ρ /σ 2 f 2σ 2 ρ /σ 2 f .
The two analysis of the error dynamics proposed in this section are based on simplifying assumptions, but provide qualitative insights on how to tune the parameter g. Additional information can be derived from the actual implementation of the algorithms. In our numerical experiments, for example, we chose a small g since we observed that the density estimation algorithm with g ≥ 1 was very noisy and would sometimes provide negative, and hence physically inconsistent, density estimates.

Optimal Sensor Placement

The second problem we tackle in this paper is Optimal Sensor Placement, namely, the problem of deciding the position of sensors yielding a good trade off between performance and cost. Since assessing in a theoretical way the performance of algorithms for density reconstruction is difficult due to the nonlinear nature of the system, we simplify the setting and limit our attention to estimation of cumulative flows, namely, of the total outflows from the cells. The resulting, static, problem is then considered as a proxy for the more complicated problem built on the dynamic density model.

We start by deriving some properties of cumulative outflows, and we proceed describing a simple linear model for cumulative flows estimation. This will help us formalizing the Optimal Placement problem.

Linear flow constraints

Let f e := t1-1 k=t0 f out e (k) be the cumulative outflow from cell e, namely, the total flow through the cell over the period of time [t 0 , t 1 ], and let f = f 1 . . . f |E| T ∈ R E + be the vector of cumulative outflows. By integrating the system's dynamics we have

e (ρ e (t 1 ) -ρ e (t 0 )) = j∈E R je f j -f e , e ∈ E \ R i
Assume now that [t 0 , t 1 ] is a period of time whose duration is high enough, and that at both times t 0 and t 1 the number of vehicles in the network is low, for example, assume that t 0 and t 1 correspond to two consecutive midnights. Then the magnitude of the vector of differences of vehicles { e (ρ e (t 1 ) -ρ e (t 0 )} e∈E\R i is small if compared to the cumulative flows in the network, and the following relation holds approximately Lf ≈ 0 ,

where L ∈ R E\R i ×E , is the matrix obtained by removing from L = R T -I the rows corresponding to onramps. This imposes a linear constraint on the cumulative flows that we shall exploit in the next subsection.

Linear measurement model and the Optimal Sensor Placement problem

In this subsection we study the performance of a linear estimator of the cumulative outflows and we show how to formulate the problem of Optimal Sensor Placement.

Let E m ⊆ E be a set of cells in which sensors are placed. We assume the following simple linear measurement model

y = H E m f + ω f (14)
where

• y s is the measurement of the s-th sensor, namely f e + ω f s if the s-th sensor is located on cell e;

• H E m ∈ {0, +1} p×n , p = |E m |, [H E m ] se = 1 if the s-th sensor is located on cell e, and [H E m ] se = 0 otherwise, so that H E m 1 = 1 and 1 T H E m 1 = p;
• ω f is a random noise vector with zero mean and covariance matrix Σ nom , related to the measurement noise ω ϕ described in the previous sections.

For sake of simplicity, we shall often assume that the components of the noise, one for each sensor, are independent with same variance σ 2 nom , so that Σ nom = σ 2 nom I. Let now n = |E| and r = rank{ L}, and consider a matrix V ∈ R n×r whose columns are an orthonormal basis of the right kernel of LT , i.e., LT V = 0 and V T V = I. From (13) we get (approximately) f = V z for some z ∈ R r , so that the measurement model can be rewritten as

y = H E m V z + ω f . ( 15
)
Given y, consider a linear estimator of z, ẑ = K z y + q z , where K z ∈ R r×p and q z ∈ R r . The Best (minimum variance) Linear Unbiased Estimator of z corresponds to the solution to

min Kz,qz E[(z -ẑ)(z -ẑ) T ] s.t. E[z -ẑ] = 0 ẑ = K z y + q z (16)
We provide an equivalent formulation and a solution of ( 16) in Lemma 1 in Appendix. An immediate consequence of Lemma 1 is that the Best Linear Unbiased Estimator (BLUE) of f is

f = K f y = V (V T H T E m Σ -1 nom H E m V ) -1 V T H T E m Σ -1 nom y
and its error covariance is

V p (E m ) = E[(f -f )(f -f ) T ] = V (V T H T E m Σ -1 nom H E m V ) -1 V T .
The quantity V p (E m ) depends on a) the (right kernel of the) matrix of splitting ratios via the matrix V and the nominal variance of the noise ω f , two parameters that are assigned, and b) on the locations of the sensors, which is the set E m , via the matrix H E m . For this reason, we will take the magnitude of V p (E m ), measured via its trace, as our metric to measure the performance of a sensor network placed on the cells E m .

Clearly, with no additional constraints the optimal placement is to equip every cell with sensors. This is straightforward as equipping all cells means setting H E = I, and from

H T E H E = I ≥ H T E m H E m immediately descends V p (E) ≤ V p (E m ), for any E m ⊆ E.
Each device has however a non-negligible purchase and maintenance cost, which has to be considered when designing the sensor network. For sake of simplicity, in this paper we make the simplifying assumption that the cost of a network over its lifetime is proportional to its number of sensors via a coefficient c > 0, so that the cost of deploying sensors on E m is c|E m |.

We thus consider the following Problem 2 (Optimal Sensor Placement): Let G = (V, E) be a traffic network with splitting ratios R and cumulative flows noise variance σ 2 nom . Find E m which solves min

E m trace {V p (E m )} + c|E m | (17)
The optimal Êm , solution to (17), trades off between the network performance, which is measured by the trace of the estimator error covariance, and the total cost of the network. Clearly, the two have a contrasting effect on the number of deployed sensors. It is however inherently combinatorial, the optimal position of the sensor being in general hard to find and requiring an exhaustive search among all the possibilities, which is intractable even for relatively small network dimensions.

We approach the problem by proposing an heuristic that relaxes it into a convex problem. Such a strategy is described after a brief discussion on the minimum required number of sensors.

Remark 4. The choice of measuring the performance via the trace of the estimator error covariance is one of several possible choices to scalarize the problem of minimizing the minimum covariance itself. Another interesting choice is the logarithm of the determinant of the estimator error covariance, proposed and studied for example in [START_REF] Zhou | An information-theoretic sensor location model for traffic origin-destination demand estimation applications[END_REF][START_REF] Xing | Designing heterogeneous sensor networks for estimating and predicting path travel time dynamics: An information-theoretic modeling approach[END_REF], which also yields a convex optimization problem following the ideas proposed in the following sections. A further possibility is the || • || 2 norm of the covariance. A comparison between these possible choices, in light of their geometrical interpretations, is due in future research.

Minimum number of sensors

Before proposing our method for solving (17), we prove that there exists a lower bound on the number of sensors |E m | in order the trace {V p (E m )} to be finite. As it will be proven, below such number, which corresponds to the number of onramps of the system, the problem of reconstruction of flows admits infinite solutions.

To this aim, relabel the cells in such a way that onramps are the first 1, . . . , |R i | cells so that we can partition L as L = L on L nn where L nn models the mutual influence of flows on non-onramp cells, and L on models the influence of onramps on non-onramp cells.

It can be shown, and we provide a formal proof in Appendix, that the minimum number of sensors for the flow estimation problem to be solved is equal to the number of onramps |R i |.

We offer two interpretations on this fact:

• Ideal measurement scenario: assume p = |E m | ideal measurements at sensors s(1), . . . , s(p) are available, y s(i) = f s(i) , i = 1, . . . , p. Then a solution to the system of equations

L f = 0 fs(i) = f s(i) , i = 1, . . . ,
= |E m | < |R i |, then rank H T E m H E m < |R i | = rank {V }, which implies that rank V T H T E m H E m V < |R i |. However, V T H T E m H E m V ∈ R |R i |×|R i |
, so the matrix is singular, and therefore the trace of the error covariance is unbounded.

Conversely, if p ≥ |R i | then rank V T H T E m H E m V = |R i
| and the trace of the error covariance is bounded.

Relaxation via Virtual Variances

The solution that we propose is based on the observation that cells that are not endowed with sensors can be interpreted as cells in which sensors have infinite noise variance. It turns out that an equivalent formulation of ( 17) is

min E m trace V (V T ΣV ) -1 V T + cσ 2 nom 1 T Σ -1 1 s.t. Σ ee = +∞, e ∈ E m σ 2 nom , e ∈ E m (18)
In this formulation

• the term H E m is the identity matrix, namely, all cells are endowed with a sensor -except, some of them have infinite noise variance and thus provide no information;

• the second term in the cost corresponds to c|E m | as

cσ 2 nom 1 T Σ -1 1 = cσ 2 nom e∈E m 1 σ 2 nom = c e∈E m 1 = c|E m |
In other terms, (18) corresponds to assigning a virtual variance Σ ee = σ 2 e to each sensor, and decide for which it should be σ 2 e = σ 2 nom , the cells in E m that actually provide information, and for which it should be σ 2 e = +∞. We call Σ the (diagonal) matrix of virtual variances, and we let the corresponding trace of error covariance be denoted, with an abuse of notation,

V p (Σ) = trace V (V T ΣV ) -1 V T .
Our approach is then based on the intuitive idea that increasing the variance on the sensors that are not crucial for the solution of ( 17) should not have a strong effect on the performance term V p (Σ). More formally, we consider the following relaxed version of the previous problem

min Σ∈Dn trace V (V T Σ -1 V ) -1 V T + cσ 2 nom 1 T Σ -1 1 s.t. Σ ee ≥ σ 2 nom , ∀e ∈ E (19)
where D n is the set of diagonal matrices of dimension n.

We now slightly rewrite the cost. First of all, by the well known property of trace trace {AB} = trace {BA}, and by

V T V = I, it follows that trace V (V T Σ -1 V ) -1 V T = trace (V T Σ -1 V ) -1
. Second, we perform the change of variables Ω = Σ -1 . In this way, we obtain the following problem

min Ω∈Dn trace (V T ΩV ) -1 + η1 T Ω1 s.t. 0 ≤ Ω ≤ Σ -1 nom ( 20
)
where η is a tunable parameter. The choice η = cσ 2 nom is a natural one due to the previous discussion, but since η influences the relative weight of performance (penalized for low η) and cost (penalized for high η) we leave it as an additional degree of freedom. Notice, finally, that the term 1 T Ω1 corresponds to the 1 norm of the inverse of the variances, a term which is commonly used term to sparsify solutions of optimization problems.

The Virtual Variance algorithm proceeds as follows:

1. solve ( 20) and denote by Ω its solution; 2. compute Σ -1 = Ω 3. discard all cells whose virtual variance is above a fixed discard threshold T d

As explained above, if the found solution provides high virtual variances at locations where sensors are redundant then this is effectively a way to select the most important cells where to place sensors.

However, this strategy does not, in general, provide good solutions to the problem. Indeed, numerical simulations have shown that in the considered scenario the solution of ( 20) can be often interpreted as endowing all cells of the network with sensors with average virtual variance, rather than keeping it low in some of them and high in others.

In order to enhance diversity between sensors, we enrich the cost of ( 20) with a term that aims to penalize homogeneity. This is reminiscent of dissensus (as opposed to consensus) strategies in multi agent networks, in which each agent possesses a value and the goal of the network is to differentiate as much as possible such values.

In this paper, we make the following simple choice. Let W ∈ R n×n-1 be an orthonormal base of the subspace orthogonal to 1, that is, W * 1 = 0 and W * W = I. We add to the cost in (20) a term that is proportional to e -1 T W * Ω1 . Since the columns of W span the orthogonal to 1, W * Ω1 is high when the element on the diagonal of Ω, which are gathered in the vector Ω1, are different one with respect to the other, and is low otherwise.

We propose the following optimization problem

min Ω∈Dn trace (V T ΩV ) -1 + η1 T Ω1 + κe -1 T W * Ω1 s.t. 0 ≤ Ω ≤ σ -2 nom I (21)
While the optimization variable in the ( 21) is Ω, the real optimization variables are clearly its diagonal entries. Observe that a) Ω is an affine function of its diagonal entries; b) V T ΩV and 1 T W * Ω1 are affine functions of Ω; c) trace X -1 is convex in the matrix X [START_REF] Boyd | Convex Optimization[END_REF]) and e -x is convex in x; d) composition of convex and affine functions, and nonnegative sum of functions, preserves convexity [START_REF] Boyd | Convex Optimization[END_REF]. The cost function of ( 21) is thus a convex function of the diagonal entries of Ω under the natural requirement that the tunable parameters η, the total variance weight, and κ, the discrepancy weight, are chosen to be nonnegative. Since the inequality constraints are clearly convex in the diagonal entries of Ω, we conclude that ( 21) is a convex optimization problem having as optimization variables the diagonal entries of Ω.

Tuning of the parameters η and κ appearing in the previous optimization problem requires some careful selection so that the resulting virtual variances are strongly separated (in order to easily discard cells not providing sufficient information) while showing sufficiently good performance. We shall provide in the following paragraphs examples of application of the previous optimization problem for different sensor placement scenarios. We anticipate here that, by trial and error, the first parameter we tuned was κ, set to the value 20. This ensures that, as required, the virtual variances solution of ( 21) are distributed in a strongly bimodal way, with low and high values being different by several orders of magnitude. As a consequence, it is usually easy to distinguish among the two groups and discard cells whose contribution to the performance metric would be negligible. The parameter η was then chosen to show good performance/cost trade off, or according to additional constraints such as the ones described in the next paragraphs. In general, high η penalizes the number of sensors, thus yielding to solutions with higher virtual variances at the expense of poor performance.

We conclude this section observing that an important issue is sensitivity of the Virtual Variance algorithm with respect to the only datum it requires, the matrix of splitting ratios. We present in Appendix the results of a numerical study suggesting that such a sensitivity might be low in some situations.

Optimal Sensor Placement with geographical or budget constraints

In this paragraph we discuss two variations of the previous procedure, which address the additional problems of geographical constraints and of strict budget limitation.

Optimal Sensor Placement with geographical constraints

The first scenario we address is concerned with the scenario in which a) some cells cannot be equipped with sensors, for example for physical reasons, and/or b) there are subsets of cells for which either all cells are endowed with sensors, or none are. An example of the latter constraints is a multi-lane road which is modelled using parallel cells (on different lanes) and on which the traffic manager can deploy induction loop. The former are buried underground and are usually required to cover the whole carriageway. As such, if i and j are, for example, two parallel cells on the two lanes on a certain section or road, then either i and j are both equipped with a sensor (i.e., the induction loop), or not.

To include this type of constraints in the proposed procedure, let E am ⊆ E, |E am | = m, be the subset of available cells, and let H Eam ∈ {0, +1} m×n be built as in Subsection 4.2. Further, let S ⊆ P E be the set of subsets of E which must be simultaneously equipped, or not, with sensors, where P E is the powerset, or set of subsets, of E.

Then the following problem

min Ω∈D k trace (V T H T Eam ΩH Eam V ) -1 + η1 T Ω1 + κe -1 T W * Ω1 s.t. 0 ≤ Ω ≤ σ -2 nom I Ω ii = Ω jj , ∀i, j ∈ σ, ∀σ ∈ S (22)
is ( 21) once we constrain sensors to be placed on cells in E am only and simultaneous sensor placement to happen according to the constraints specified by the set S. Notice that the diagonal entries of the solution Ω are the inverse of the virtual variances on the cells in E am only. The matrix W is defined as previously, but with suitable dimension (m × m -1). As in the general problem, cells are chosen only if the corresponding virtual variance is below a certain threshold, and clearly cells that are not in E am cannot be chosen. Notice that the constraint Ω ii = Ω jj in ( 22), or the less requiring |Ω ii -Ω jj | ≤ ε, for some tunable parameter ε, are convex constraints, so (22) remains convex.

Remark 5. By the discussion in Subsection 4.3, the minimum number of sen- 22) is not well posed and the solution will only have very high virtual variances. Clearly, such a solution is not acceptable and should be discarded.

sors is r = |R i |. As such, if |E em | < |R i | the problem (

Optimal Sensor Placement with budget constraints

In this second scenario we discuss budget constraints in the form of constraints on the maximum number of chosen sensors, a very common requirement in real-case applications.

We propose a solution based on the following iterative approach:

• Initialization: set η(0) and κ to some prespecified nonnegative values; t max to the maximum number of iterations; n max to the maximum number of sensors;

• t-th step -The problem ( 21) is solved with η = η(t); -Let n(t) be the number of sensors in the solution provided by the Virtual Variance algorithm. Then * If n(t) ≤ n max , or if t ≥ t max , the procedure stops; * Otherwise, set η(t + 1) = g(η(t)), where g is a monotonically increasing increasing function, and the procedure iterates.

The rationale behind this procedure is that, as previously discussed, η penalizes a low total sum of the virtual variances. Therefore, by iteratively increasing η the solution to (21) will tend to exhibit higher and higher virtual variances, thus reducing the number of sensors. Remark 6. Once again, and related to the discussion in Remark 5, the specified maximum number of sensors cannot be be less than r = |R i | by the results presented in Section 4.3. If this is not the case, numerical experiments show that the algorithm either simply iterates until the number of iterations reaches t max , or the found solution exhibits an extremely high trace (V T H T Eam ΩH Eam V ) -1 -due to the fact that the internal matrix is (numerically) almost singular. As in Remark 5, such a solution is not acceptable and should be discarded.

Numerical experiments

Numerical Experiments for the Optimal Sensor Placement

In this subsection we present the results of two numerical experiments. In the first, we solve the problem of Optimal Sensor Placement in a small (but not trivial) regular grid, in the second, we apply the procedure to the real-world case of the Rocade Sud.

(Please place Figure 7 about here.)

Regular grid

Consider the regular grid composed of 25 cells shown in Figure 7. We want to solve the problem ( 17) by ( 20) with parameters σ 2 nom = 1, nominal sensor variance, and c = 1, cost of a single sensor.

The network is small enough to run an exhaustive search to solve the problem (17). In particular, for each h = 4, 5, 6, . . . , 21, we try all possible combinations of h = |E m | sensors, thus finding the one that minimizes V p (E m ).

We also run our Virtual Variance algorithm with total variance weight η = 2 and discrepancy weight κ = 20. Further, we set the threshold on the virtual variance for discarding a sensor to T d = 100.

We illustrate the results in Figures 7 and8. In Figures 7 the 8 cells chosen by the Virtual Variance algorithm are represented as a green dot marks, and the best possible placement with 8 sensors (found by exhaustive search) as red dot marks. It can be noticed that both procedures place most of sensors at the boundaries of the network.

Figure 8 shows instead the total cost of the best placement obtained through exhaustive search for h = 4, . . . , 21, and the total cost found by the Virtual Variance algorithm with 8 sensors. By total cost we mean the sum of estimation performance and network cost V p (E m ) + c|E m |, which, notice, is not the metric that is used in the Virtual Variance algorithm. Nonetheless, it is appreciable that the Virtual Variance algorithm not only provides a solution whose number of sensors is close to the global optimum (which is with 6 sensors), but also that, using 8 sensors, the Virtual Variance algorithm places them almost in the optimal way.

(Please place Figure 8 about here.)

Rocade Sud

Our second experimental setting is the Grenoble Traffic Lab (GTL), a network of sensors deployed for monitoring and research purposes along the Rocade Sud, a peri-urban 10.5 km long freeway connecting the two highways A41 (north-west) to A480 (south) in the town of Grenoble in the south of France, see Figure 9. The network is composed of 135 magnetometers buried in the ground on both lanes along the main line every 250 meters (on average), on each onramp and offramp, and on three connectors from urban roads to three onramps, for a total of 68 sensing locations. For our purposes, each sensing location will correspond to one pair of sensors. For a detailed report on the GTL, we refer to (Canudas de Wit et al., 2015). Furthermore, and for sake of simplicity, while in the real network sensors are deployed in pairs, we shall assume from now on that each sensing location has only one sensor (as we shall always discard sensor speed measurements).

Figure 9 shows the position of each of the 22 sections of the main line in which there are sensing locations on both slow and fast lanes (and usually a ramp). In the same figure we also show a stylized representation of the freeway, including ramps and queues, the positions of the 68 sensing locations, and the distance between consecutive measurement sections along the main line.

We partition the Rocade in cells in such a way each cell includes one sensor, so that in Figure 9 each numbered circle also corresponds to one cell.

Here, we do not consider onramps and offramp, limiting our attention to the main line of the Rocade Sud. The reason for this choice is that the Rocade has 10 onramps along its main line, which, summed to the two cells in the very first section of road, imply a minimum number of sensor of 12 by the discussion in Subsection 4.3. While this number is not high on its own, numerical experiments not reported in this paper have shown that good estimation performance require a number of sensors which is too high in most realistic (i.e., non academical) implementations.

The corresponding reduced graph representation, essentially made of several groups of parallel edges, consists of 46 cells. We provide a stylized version of it in Figure 10.

(Please place Figure 9 about here.)

The matrix of splitting ratios could be estimated from the data, but we make here a simpler assumption and assume that vehicles split according to the following rule: vehicles on slow lane cells remain on slow lane or turn into fast lane with a 70%-30% rule, and analogously for vehicles on the fast lane. If the next section has three parallel cells (cells 22-23-24 and cells 66-67-68), vehicles spit uniformly in such three cells.

We did not run an exhaustive search due to the relatively high dimension of the network and the consequent relevant computational load. Our comparison is instead with the locations of fixed loops installed by the Government Agency Centre national d'information routiére (CNIR) (Centre national d'information routiére (CNIR), 2014), which correspond to cells marked with a red dot in Figure 10. We shall show that the latter are positioned in a way that is in good accordance with the results of the Virtual Variance algorithm, even though our procedure allows for a slightly better design.

We run the Virtual Variance algorithm in three scenarios: 1) unconstrained scenario with total variance weight η = 0.2 and discrepancy weight κ = 20; 2) unconstrained scenario with η = 1 and κ = 20; 3) constrained scenario with number of sensors at most 10, initial η = 0.2, and κ = 20. We assume that σ 2 nom = 1 and that the cost per sensor is c = 1. In all cases, sensors are constrained either to be present in both lanes on a same section, or to be absent.

The results are summarized in Figure 10 and Table 2. The tables shows the number of sensors in the solution computed via the Virtual Variance algorithm, as well as the trace of the corresponding estimation error covariance V p (E m ) and the total cost V (E m ) = V p (E m ) + c|E m |. In Figure 10, cells found in the unconstrained scenario with η = 0.2 are marked with a green dot. As can be seen in Table 2, our algorithm requires one less sensors than the network deployed by CNIR and in addition the trace of the error covariance is smaller.

In the constrained scenario and in the unconstrained scenario with high η (which, as explained above, indirectly penalizes the number of sensors), the trace of the error covariance V p (E m ) increases, as expected. Furthermore, the chosen cells in the latter two cases are subsets of the cells chosen in the unconstrained case: in particular, in the constrained scenario all cells are kept except 13, 14, 49 and 50, and in the unconstrained scenario with high η the algorithm further discards cells 64 and 65.

(Please place Figure 10 about here.) (Please place Table 2 about here.)

Density reconstruction -experimental results

We provide here numerical results of the implementation of the data fusion algorithm for density reconstruction on data from the Rocade Sud.

On each sensing location and every T = 15 seconds, the system counts the number ϕ m e of vehicles that crossed the location, their average speed v m e , and the average occupancy o m e of the location. Since the latter is approximately proportional to the density of vehicles, so we shall assume that sensors can directly measure densities. In addition to fixed sensors, we use Floating Car Data provided by INRIX. Following INRIX schema, the Rocade has been further partitioned into FCD segments. One measurement of average speed is available on each FCD segment every 1 minutes. FCD segments partition the whole main line of the Rocade and include most onramps and offramps, but single lanes are not distinguished along the main line. FCD segments are represented in Figure 9 as rectangles encircling several sensing locations/cells.

For our experiments, we employ the sensor configuration obtained in the previous section via the Virtual Variance algorithm with η = 0.2. In particular, and in order to prove that the method shows good performance even with sparse equipment, we only use the sensors on the main line, which are, for reference, shown in light blue in Figure 9. Further, we don't use any information on flow or speed on the ramps.

Calibration of the Fundamental Diagram was performed via the algorithm described in Paragraph 3.1.1 and using the data from the GTL sensor network from April 10th, 2014, a working day (a Thu.) exhibiting very standard traffic pattern:

• very limited night time traffic;

• a peak of congestion in the morning (8:00 -10:00), triggered by vehicles exiting towards the city from the Rocade at the offramp of Eybens (cells 37/38) and spilling back until Meylan (cells 2/3), and a second, smaller peak of congestion triggered by vehicles entering in A480 at Rondeau (cells 66/67/68) but blocked by the high traffic on A480, and spilling back until around Libération (cells 61/62);

• a third, smaller, congestion triggered around Eybens around 14:00-15:00;

• in general, medium/heavy but fluid traffic from 10:00 to 16:00

• a second peak of congestion in the afternoon, again triggered by congestion at Rondeau at around 16:00, spilling back on the whole freeway in around 60 minutes, and lasting approximately two hours.

We considered an average vehicle length of ave = 5 meters, close to the value provided in [START_REF] Treiber | Traffic Flow Dynamics[END_REF]. This yields a jam density value of ρ jam = 200 veh/km.

(Please place Figure 11 about here.) As for the previous section, we consider the following simple rules to set the matrix of splitting ratios

• let e be a fast lane cell. Then 70% of vehicles continue on the fast lane cell and 30% turn into the slow lane cell; if in the following section there are three parallel cells, vehicles split uniformly;

• let e be a slow lane cell. If among the following cells there is not an offramp, then 70% of vehicles continue on the slow lane cell and 30% turn into the fast lane cell; if in the following section there are three parallel cells, vehicles split uniformly. Otherwise, 20% of the flow is directed towards the offramp, and the rest splits as previously specified;

• if e is an onramp cell and j is the following slow ramp cell, then R ej = 1.

• if e is a queue cell and j is the following onramp cell, then R ej = 1.

In words, vehicles split according to a 70%-30% lane-change rule in the cells on the main line, and at each offramp approximately 10% of vehicles exit from the freeway.

Implementation

To assess our method, we considered the whole month of April 2014 (except April 13th, a Sun., for which FCD measurements were not provided). A typical result of calibration of the Fundamental Diagram is illustrated in Figure 11,which shows in thick black the linear-quadratic Fundamental Diagram, in dashed thick black the corresponding standard linear Fundamental Digram in congestion regime, and as crosses the pairs (density, flow) measured on a day different from that used for calibration, in this case April 24th, 2014. As standard and well known, data in freeflow regime are in good accordance with the linear part, while data in congested regime are much more scattered and more difficult to fit. As it can be noticed, the standard bilinear Fundamental Diagram overestimates the flows in congested regime (the dashed think line is on average higher than the corresponding pairs (density, flow)), while the quadratic curve seems to better capture the average flow-density relation. Nonetheless, it is clear that the so found curve is only a very crude approximation of such a relation, which might be better captured using a stochastic description [START_REF] Jabari | A stochastic model of traffic flow: Gaussian approximation and estimation[END_REF]. Investigation of the latter possibility will be the focus of future research.

The proposed algorithm was implemented in Matlab on a non dedicated commercial laptop with 2.1 GHz i7-4600U CPU and 8 GB RAM. Optimization problems, required both in offline and online steps, were solved using standard Matlab functions as well as the modelling and optimization system CVX (Grant andBoyd, 2013, 2008). The time required for calibration of the Fundamental Diagrams is between 30 and 40 seconds for each cell, while reconstruction of all the samples for a whole day requires less 10 minutes, averaging 100 ms per sample.

(Please place Figure 12 about here.) (Please place Figure 13 about here.) Typical results are reported in Figures 1213. For validation purposes only, density and flow measurements from all GTL fixed sensors are considered ground truth. As such, the left panels show the evolution of the "true" measured density and flow in all the cells on the main line, over the whole day, for each day. On the x-axis, the 46 sensing locations along the main line (numbers correspond to the labels in Figure 9), on the y-axis, the 5760 time slots over the whole validation day (one slot every T = 15 seconds). The chosen colors range from green (low density/flow) to yellow (medium/critical density/flow) to red (high density/flow). In the right panels, we show, using the same legend, the results of density and flow reconstruction. As it can be observed, the estimation algorithm is able to represent the four congestion events described in the previous section in a reasonably good way, given the limited amount of information employed; in particular, observe that the two small congestions at Rondeau during the morning and at Eybens during early afternoon, when present, are both detected. Notice that the resulting estimate remains good and reasonably close to the real profile, despite the absence of flow measurements on ramps, which in the real system have a strong impact on the traffic dynamics, especially in peak hours. We further observe that the observer is able to capture different regimes in parallel cells. An example of such a phenomenon is the succession of congested/non congested cells 28-29-32-33 at 16:00 (Saint-Martin-d'Héeres -Saint-Martin-d'Héeres entrance), which correspond to slow/fast lanes.

(Please place The performance of the algorithm is quantitatively illustrated via the absolute error between measured and estimated flows and densities (8). The results are summarized in Table 3 andTable 4, which report the maximum absolute error δ on the 75%, 90% and 95% of the pairs (cell, time slots), for all cells on the main line and all samples during a day, for all the considered days, and for both densities and flows. In Table 3 time slots correspond to the original 15 seconds time slots. In Table 4 we consider time slots of 5 minutes and we average in time all quantities. For the original slots, the table shows an average error of less than 11 veh/km for the 75% of pairs and less than 23 veh/km for the 90% of the pairs. In 5% of the (cell, time) pairs the error is higher but still less than around 35 veh/km. The magnitude of the errors is influenced by the high sampling rate of the system, which, as it is well known, increases the noise power. This is confirmed by a comparison with Table 4, for the averaged system, which shows errors on average 40% less for flows and 25 -30% less for densities.

We conclude with a comparison with an an oracle that knows

• the exact outflow f out e (t) for every cell e and every sample time t (namely, for each 15 seconds time slot);

• whether the cell is in freeflow or in congestion, for every cell e and every sample time t.

The oracle then runs the same estimation algorithm as described above. Notice that by exploiting its available information, the oracle is able to propagate the flows in the network with zero error (i.e., ∆f out = 0) and always knows whether a cell's density is low (if the cell is in freeflow) or high (if the cell is congested). Consequently, the only source of noise that affects the oracle is the model uncertainty which is inherent in the Fundamental Diagram. The average error obtained by the oracle is 5.8 for the 75% of the pairs (cell, time) and 12.5 for the 90% of the pairs, even with the high amount of additional (and precise) information available to the oracle. Indeed, this confirms that estimation in traffic systems is a difficult task, and that errors of absolute magnitude around 10-20 veh/km can be acceptable, as they capture the qualitative trend features of the traffic system -such as low density, medium density, high density. On the other side, estimate of flows are rather low, being less than 1 ∼ 2 vehicles for the 75% of the pairs (cell, time), and less than 2 ∼ 3 for the 90% of the pairs. The biggest difference between estimated and measured densities can be observed at Eybens exit (cells 37-38), where the estimated flow constantly predicts a higher density than the measured one. The explanation is straightforward: Recall that we do not use ramp data in order to show the prowess of our method even employing a very small number of sensors. However, as mentioned in the description of the data, the exit of Eybens is a critical point whose ramp is selected by a high fraction of vehicles to exit the Rocade towards the town; in addition, the corresponding cells belong to a long FCD segment running from Saint-Martin-d'Hères (cells 32-33) to Eybens entrance (cells 41-42), which provides, during peak time, just one set of rather low speed measurements, which do not distinguish between the stretch of road before Eybens exit (congested and at low speed) against that after Eybens exit (uncongested and at high speed). Due to the so obtained low speed measurement, the algorithm tends to estimate a high number of vehicles along the whole segment, instead of two different regimes before and after the ramp. Analogously, we observe a mismatch between measured and estimated flow in the first and last sections, which are due to unobserved flows from onramp and to offramps. A second discrepancy is the smoothness of the reconstructed density and flow, as compared with the more scattered measurements. The latter is due to the high measurement rate, which during stop-and-go phenomena results in measurements which rapidly oscillate between stopped vehicles and low or medium speed. On the converse, the optimization based flow reconstruction and the first order mass conservation law for densities have a low pass effect therefore producing smoother, more regular patterns. Further research direction will investigate the possibility to detect stop-and-go phenomena and reproduce, at least qualitatively, the resulting irregular patterns.

Conclusions

This paper addressed the problems of fusion of heterogeneous sources of information for density estimation in Road Transportation Networks and optimal sensor placement via a heuristic that we called Virtual Variance algorithm. A gradient descent procedure for the calibration of the Fundamental Diagram is also discussed. Efficacy of the proposed solutions is illustrated on a regular grid and on the real world scenario of the Rocade Sud in Grenoble. Future research directions include and are not limited to estimation of statistical properties of measurement noises from real data and development of stochastic models for the relation between flows, speed and densities, optimization of the observer's parameters for minimization of mean-square reconstruction error, calibration of the matrix of splitting ratios, tuning of the weights in the Virtual Variance algorithm, and extension of the optimal sensor placement strategy to maximize density reconstruction performance. a row-stochastic matrix. The average of the 100 instances of the cost V (E m ) (namely, the estimator variance) is m cost = 6.93 with variance σ cost = 0.28246. The number of chosen cells varies between seven and nine. In Figure 14 we plot the distribution of cell usage. The figure is to be interpreted as follows: Four cells are essentially always used. Six cells are used between half and two-thirds of the times. Four other cells are used sporadically, the rest are never used. This suggests that there is a bulk of around ten cells that is particularly informative for this grid-like structure with splitting ratios specified in the manuscript. Whether this is a generic feature of the Virtual Variance algorithm, or a consequence of the very regular pattern of the chosen grid will be investigated in future research.

10

  Stylized representation of the main line part of the Rocade Sud. White ovals represents junctions of the graph. The selected positions of 17 fixed sensors by CNIR are marked by red dots, those of the 16 sensors chosen by the Virtual Variance algorithm by green dots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Calibration of the Fundamental Diagram on the cell Eybens exit -slow lane. The linear-quadratic Fundamental Diagram, calibrated using data from April 10th, 2014, is shown in thick line. The dashed thick line represents the corresponding linear Fundamental Diagram in congested regime. Each cross is a (flow, density) pair measured on April 24th, 2014, one for each time slot of T = 15 seconds during the whole day. Flow is the number of counted vehicles crossing the sensing location, density is the measured density of vehicles, during the time slot. . . . . . . . . 12 Numerical results of density estimation on all cells: measured densities (left panel) and estimated densities right panel). Day: 24-04-2014 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Numerical results of flow estimation on all cells: measured densities (left panel) and estimated densities right panel). Day: 24-04-2014 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Distribution of cell usage resulting from the virtual variance algorithm applied to 100 different instances of perturbed matrix of splitting ratios. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  e a e ρ 2 + b e ρ + c e , ρ c e < ρ ≤ ρ jam where • ρ c e is the critical density. It partitions the set of densities [0, ρ jam ] into the freeflow low-density region [0, ρ c

  0 is the freeflow speed on cell e; the value C e = v ff e ρ c e is the capacity of the section or road, namely, the maximum number of vehicles that can flow through it during a period T ; • we assume that the Fundamental Diagram congested region is a quadratic function of the density. The following relations among the parameters a e , b e and c e hold for consistency      a e (ρ c e ) 2 + b e ρ c e + c e = v ff e ρ c e a e ρ jam e 2 + b e ρ jam e + c e = 0 a e ≥ 0 Remark 2. Recall that we defined ϕ e (•) to be the flow of vehicles through a section during a sample time T . As such, the units of the speed v ff e are km/T . Remark 3. A standard choice for the Fundamental Diagram is the a triangular Fundamental Diagram, for which in the congested region ϕ e (ρ e ) = ω e (ρ jam e -ρ e )

  f out e (t).receives the measurements {v m e (t)} e∈E when available, or holds the last speed measurements received;-For each cell e, computes the two possible densities ρ 1 e (freeflow) and ρ 2 e (congested) corresponding to flow f out e (t) assuming that the local flow ϕ e = ϕ e (ρ e ) is exactly determined by the Fundamental Diagram;-For each cell e, computes the two velocities v e (ρ 1

Figure 1 :Figure 2 :

 12 Figure 1: Flowchart of the paper's contributions: Optimal Sensor Placement for sensor network design, Offline Fundamental Diagram Calibration, and Online Density Estimation.Curvy blocks represent information provided to the three procedures.

Figure 3 :

 3 Figure 3: Roads are partitioned in cells and FCD segments. The depicted FCD segment (in dot-dashed line) includes, among others, the three cells e, j and k, so that s(e) = s(j) = s(k).

Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 45678 Figure 4: Comparison of speed measurements from sensors and Floating Car Data at location Taillat on Grenoble Rocade Sud, 24/04/2014.

Figure 9 :Figure 10 :Figure 11 :Figure 12 :Figure 13 :Figure 14 :

 91011121314 Figure 9: The experimental setting: the town of Grenoble and the Rocade Sud and a stylized version of the freeway. The positions along the main line of the 22 sections of the main line in which sensors have been placed is shown as red pin. The positions of the 68 fixed sensing locations are shown in the stylized map. Each sensing location also corresponds to a cell. Light blue circles denote fixed sensors that are selected by the Virtual Variance algorithm and are used in the implementation of the Density Reconstruction algorithm. Each rectangle represents one FCD segment, often providing average speed measurements over more than one cell.

  p is a candidate vector of cumulative flows. Then, if p ≥ |R i | the system has a unique solution which is the true vector of flows. Conversely, if p < |R i |, the system is undetermined; • Noisy measurement scenario: assume p = |E m | noisy measurements at sensors s(1), . . . , s(p) are available and assume to adopt the Best Linear Unbiased Estimator presented above to estimate the flows. If p

  Table 3 about here.) (Please place Table 4 about here.)

Table 1 :

 1 cells equipped with sensors V p (E m ) performance of sensors in E m List of symbols.

	Graph notation and modeling parameters
	V		set of junctions		E	set of cells
	R i		set of onramps		R o	set of offramps
	E -e		set of cells preceding cell e E + e	set of cells following cell e
	T		sampling time		t	time index
	e R ek splitting ratio from e to k length of e	ave L	average vehicular length non-onramp sublaplacian
	λ e		inflow into onramp e		v limit e	speed limit on e
	d e (•) demand function on e		s e (•)	supply function on e
	ϕ e		Fundamental Diagram on e ρ jam	jam density
	v ff e		freeflow speed on e		ρ c e	critical density on e
	C e		flow capacity on e		ω e	wave speed on e
	Variables and measurements
	ρ e		density on e		o e	occupancy on e
	f in e		inflow in e		f out e	outflow from e
	v e		speed on e		f e	cumulative outflow from e
	ϕ m e v m s	flow measurement on e speed meas. on segment s v m ρ m e e	density measurement on e speed measurement on e
	ρe		density estimate on e		ρe	density pseudo-meas. on e
	f in	e inflow estimate on e		f out	e	outflow estimate on e
	Algorithms parameters
	Γ		density observer gain		γ	data fitting parameter
	E c		cost of a sensor		Σ ee	inverse of Virtual Variance on e
	η		total variance parameter	κ	homogeneity parameter
			Scenario	η	optimal # sensors	V p (E m ) V (E m )
			Fix			17	3.8161 20.8161
			Unconstrained	0.2		14	3.6867 17.6867
			Unconstrained	1		10	5.1732 15.1732
			Constrained, # ≤ 12 0.32		12	4.4972 16.4972

m

Table 2 :

 2 Results of the four considered scenarios.

While we employ SI units for simplicity and for coherence with the data in the Grenoble Traffic Lab, the presentation of our results would obviously be unchanged if other systems of measurements, such as the imperial system, were used.

Here we use that for x ∈ R N it holds true ||x|| 1 ≤ √ N ||x|| 2 .

Acknowledgements

The authors gratefully thank Rene Fritz, Mike Corlett, and INRIX Europe for providing the Floating Car Data that led to the results presented in this work. We also thank Roland Dollet and Sylvain Nachef for introducing us to the topic of Density Reconstruction in Road Transportation Systems.

Appendix

An equivalent formulation for (16) and its solution

The following Lemma formulates a problem that is equivalent to ( 16). The proof is straightforward.

Lemma 1. The solution (K z , q z ) to (16) is given by q z = 0 and K z the solution of

Proof of minimum number of sensors Recall that L = L on L nn is the lower part of the matrix L = I -R T once we relabel the cells in such a way that the onramps are the first

) in which the roles of cells and junctions are reversed, and in particular in which

The following result is adapted from [START_REF] Lovisari | Stability of monotone dynamical flow networks[END_REF].

Lemma 2. Let G = (V, E) be a graph and J ∈ R |V|×|V| be a sublaplacian of G. Then all the eigenvalues of J have negative real part except possibly eigenvalues in 0. Moreover, if S is the set of nodes v for which u J vu < 0, then J is invertible if for every u there exists a directed path in G from u to a node v ∈ S.

In the case under analysis, we take J = L T nn and S to be the set of cells directly following an onramp, namely S = {e ∈ E : ∃j ∈ R i : R je > 0}. Recall that by assumption for every cell e ∈ R i there exists an origin j ∈ R i and a path from j to e, so there must also be a k ∈ S and a path from k to e (at most being k = e), so the assumptions of Lemma 2 are satisfied. This establishes that L nn is invertible, and therefore that L is a full row-rank matrix with rank |E \ R i |. Since the number of columns of L is |E|, it follows that its kernel has dimension r = |R i |, and thus rank {V } = |R i |.

(Please place Figure 14 about here.) Sensitivity of the Virtual Variance algorithm with respect to splitting ratios

In this section we present numerical results suggesting that the sensitivity of the Virtual Variance algorithm, namely, the sensitivity of sensor positioning as a function of small variations in the splitting ratios, might be low in some scenarios. We have considered the regular grid presented in Section 5.1.1 and run 100 times the Virtual Variance algorithm with η = 2, each time adding zero mean gaussian noise with variance σ = 0.1 to each nonzero entry of the nominal matrix of splitting ratios, and then renormalizing the latter to obtain again