Dynamic Route Flow Estimation in Road Networks Using Data from Automatic Number of Plate Recognition Sensors

Abstract

The traffic flow on road networks is dynamic in nature. Hence, a model for dynamic traffic flow estimation should be a very useful tool for administrations to make decisions aimed at better management of traffic. In fact, these decisions may in turn improve people’s quality of life and help to implement good sustainable policies to reduce the external transportation costs (congestion, accidents, travel time, etc.). Therefore, this paper deals with the problem of estimating dynamic traffic flows in road networks by proposing a model which is continuous in the time variable and that assumes the first-in-first-out (FIFO) hypothesis. In addition, the data used as model inputs come from Automatic Number of Plate Recognition (ANPR) sensors. This powerful data permits not only to directly reconstruct the route followed by each registered vehicle but also to evaluate its travel time, which in turn is also used for the flow estimation. In addition, the fundamental variable of the model is the route flow, which is a great advantage since the rest of the flows can be obtained using the conservation laws. A synthetic network is used to illustrate the proposed method, and then it is applied to the well-known Nguyen-Dupuis and Eastern Massachusetts networks to prove its usefulness and feasibility. The results on all the tested networks are very positive and the estimated flows reproduce the simulated real flows fairly well.

1. Introduction

The traffic management in road networks is a current problem that administrations, academic and technical organizations are always facing due to the serious problems that it may cause, such as climatic pollution, congestion, increase in accidents, etc. Today, in order to carry out any management operation or mitigation of these effects, it is required to know the traffic behavior—which is sometimes very difficult due to its dynamic nature and because it is the result of users’ decisions—which uses the infrastructure for its activities. Despite this, due to its computational implementation and operational simplicity, the models most used by policymakers, with a view to transportation planning and policies related with sustainability, try to reproduce the traffic conditions from a static or steady-state approach (Sheffi [1]; Ben-Akiva and Bierlaire [2]). Dynamic traffic estimation tools, of which the first applicable models were developed years ago (Daganzo [3] or Yperman [4]), have particularities both from theoretical and computational points of view that make them more complex and require to be defined and implemented in greater detail by using a greater number of resources; however, they permit the knowledge of space-time traffic conditions, being very useful to make some management decisions, even from a sustainability point of view (see Wang et al. [5] for a compressive and detailed review about environmentally sustainable applications with traffic simulation models). Some application examples about the use of dynamic traffic models for this purpose are: (i) to help users in their decisions about departure times and route choice according to the existing traffic on the network to lower the congestion or the harmful emissions, (ii) to reproduce realistic flow phenomena such as shock wave and queuing (due to acceleration and stop-and-go conditions) that improve the accuracy of vehicular emission estimations and (iii) the design and management of environmentally sustainable transportation networks, among others.

Since the model proposed in this paper is related to dynamic traffic estimation, a deeper literature review is developed in the next section.

2. Literature Review and Contributions of This Paper

Dynamic traffic models try to estimate or reproduce the flow patterns, in time and space, over a network for given demand conditions. To achieve this, the approach that can be used by these models is two-fold: If the demand conditions are known, it would be enough to use a dynamic traffic assignment model (DTA) to obtain the estimation of flow patterns in the elements of the network as output. If the demand conditions are not known a priori, they can be estimated from observed traffic data by means of some technique, to later be used in the subsequent estimation of the flow patterns.

2.1. The Dynamic Traffic Assignment (DTA)

In terms of time and space treatment, DTA models can be mainly classified as discrete or continuous. According to some authors (Smith [6]; Daganzo [3,7]; Cantarella and Watling [8]; Iryo [9]), discrete models use nets to model space and time, trying to reproduce the dynamic behavior of the flow that occurs in a discrete time of a given period of the day. On the other hand, continuous dynamics approaches model the space and time as continuous variables. These models allow to reproduce many of the real traffic phenomena, such as variability over time of capacities and demand flows, queue formation and dissipation at bottlenecks, etc. However, since they are mainly based on hydraulic analogies or kinematic waves, very complex models are required to reproduce these situations (Lighthill and Whitham [10,11]; Newell [12,13,14]). In addition to the two previous groups of models, the literature includes others that allow the properties of both to be combined simultaneously (see Nogal [15] for more detail).

To solve the DTA either for a discrete or continuous approach, we need (i) an assignment strategy that models how users choose or select their departure time and routes, usually according to a choice model (Cantarella et al. [16]), and (ii) a network performance model, often referred to as Dynamic Network Loading (DNL), which allows to reproduce the flow propagation over the network and relationships between flows and travel delays given a route inflow obtained from a previous assignment strategy.

In turn, according to Han et al. [17], a DNL model has to have some of the following components: (i) some form of link and/or route dynamics, (ii) analytical relationships between flow/speed/density and link traversal time, (iii) flow propagation constraints, (iv) a model of junction dynamics and delays, (v) a model of route traversal time and (vi) appropriate initial conditions. In particular, the analytical relationships should try to capture the driver behavior through the following (Szeto and Lo [18]):

• First-in-first-out (FIFO) rule: vehicles that enter the link earlier will leave it sooner.
• Causality: speed and travel time of a vehicle on a link are only affected by the speed of vehicles ahead.
• Queue spillback: end of the queue spilling backwards in the network.

Regarding DNL models, the work of Dell’Orco [19] can be highlighted, who formulates a model with the mesoscopic approach for the flow propagation study on a network, with the subsequent object to join it with modelling of traffic pollution. Osorio et al. [20] propose an analytical model of dynamic queues that takes into account the rebound and captures the backward wave in congested conditions, with boundary conditions that allow dynamic traffic modelling in a network. Di Gangi et al. [21] propose a dynamic model with the mesoscopic approach combined with a strategy for traffic control, that allows to faithfully reproduce the queue spillback phenomena. Castillo et al. [22] make some corrections in previous DNL model assumptions to satisfy the FIFO rule and reproduce traffic flows more faithfully in a model with the macroscopic approach, adding a minimum delay to consider the queue of vehicles produced by congestion in subsequent links. Ge and Fukuda [23] propose a DNL model based on the reservoir method, within which the flow in links is modeled through partial differential equations based on introducing a vehicle tracking method as an innovation. Han et al. [17] propose a multistep model that integrates a Dynamic User Equilibrium (DUE) assignment model with simultaneous choice of route and departure time and a DNL model formulated as a system of differential algebraic equations (DAEs) that involves the variational theory for kinematic wave models. Their algorithm is tested in several traffic networks with different sizes.

About the assignment strategy that has to be solved prior to the network loading, the dynamic Origin-Destimation (O-D) matrix can be assigned based on several approaches, such as for example program-based or variational inequality formulations (Smith [24]; Janson [25]). Here, unlike the static ones, the cost results depend on the route and on the departure time selected by users. In any case, it can be said that these are basically an extension of the static models.

However, quantifying pre-fixed trips’ proportions (i.e., obtaining the dynamic route flows) through the dynamic O-D matrix by means of a route choice criterion could be difficult if there is insufficient time-dependent field information of the user’s behavior (which traditionally has resulted as a complicated task). Traditionally, transportation agencies could provide this information through time-dependent trip tables, derived from surveys or models, that describe trips and departure times within a relatively short time interval (Chiu et al. [26]).

2.2. Dynamic O-D Matrix Estimation from Field Data

From the 1980s to date, special attention has been paid to the possibility of estimating the previously mentioned dynamic O-D matrices from observed traffic data. Such data can come from different types of sensors located on the network, tolling stations, probe vehicles, etc. (Aulebert et al. [27]). Specifically, sensors and their associated technology can become an easy tool to deploy and can be very efficient in collecting this data, although variables such as economic cost, utility, ease of installation and type of data to collect determine the choice of the sensor (Bernás et al. [28]).

However, since the beginning, the most widely used data for these models have come from the counting sensor, intrusive or non-intrusive on the road network (Guerrero-Ibáñez et al. [29]), such as inductive loop detectors, magnetic detectors, piezoelectric detectors, video camera, among others. These types of sensors only quantify and classify the vehicles that interact indirectly with them, and from their data it is possible to quantify time-dependent traffic intensities in a section of the road network, which can be used to develop more general traffic and mobility studies (Pulugurtha and Nambisan [30]; Won [31]; Krivda et al. [32]; Wang et al. [33]). Generally, the techniques used for this O-D estimation are based on a joint approach of this step with a subsequent DNL model through a bi-level model (Bert [34]), which can involve the use of least squares optimization problems, entropy maximization or Kalman filters for the O-D estimation step (Balakrishna et al. [35]; Kim and Jayakrishnan [36]).

Despite this, the dynamic estimation of O-D matrices with data from counting sensors is complicated due to the under-specified nature of the problem. Even for static problems, the number of O-D pairs are usually much larger than the number of installed sensors, and hence, there is an infinite set of solutions satisfying the conservation laws (Castillo et al. [37]; Mínguez et al. [38]).

To solve these difficulties, the development of new and more advanced sensor technologies has made data collection more efficient and they provide more information than the simple vehicle count (Djukic et al. [39]). In this new field, Automatic Vehicle Identification (AVI) sensors come to collect data that allow for more accurate studies, because they collect a type of information that no other sensor can collect in addition to the simple fact of counting and classifying vehicles by groups. In particular, time stamps to quantify travel times, and characteristics or identification of each vehicle to reconstruct the route flows. One example of this type of sensor is the Automatic Number Plate Recognition (ANPR) system, which is becoming very useful for conducting real-time traffic and mobility studies on larger-scale urban and inter-urban road networks at a moderate cost (see Álvarez-Bazo et al. [40]). The usefulness of the data collected by these sensors has already been tested in studies aimed at the dynamic estimation of O-D matrices and adjustment of dynamic models (Dixon and Rilett [41]; Vaze et al. [42]; Robinson and Venter [43]; Liu et al. [44]), as well as estimating the evolution of traffic over time based on data collected by different sources (Hadavi et al. [45]; Song et al. [46]).

2.3. Contributions of This Paper

One of the main topics that has not yet been deeply studied is the capacity of the use of the data collected by the sensors for the direct observation and dynamic reconstruction of route flows. That is, from the data collected and with only the use of a DNL model, being able to estimate and reconstruct the variation in time of the intensities of the route flows. The interest of observing and estimating the dynamic flows in routes with AVI sensors is where the assignment strategy step of O-D flows to routes is being implicitly carried out, an advantage that could not be done with counting sensors.

Among the few references, Rao et al. [47] carried out an analysis of within-day variations of O-D patterns from data collected with sensors that incorporate the ANPR system, which is capable of estimating the probability of the route followed by a vehicle among all the possible candidate routes within an area or traffic network. Their results deduced that with a sampling or observation of 60% of the flows, it is possible to accurately estimate the O-D patterns. Chen et al. [48] propose a model for estimating flows on routes from a combination of data collected by probe vehicles and AVI sensors, using a neural network model that despite offering good estimation results can be improved in certain aspects, with one being the dynamic route choice in congested networks.

With all this in mind, this paper proposes a method for estimating dynamic route flows in road networks to then obtain the rest of the flows by means of the application of a Continuous Dynamic Network Loading (CDNL) model. For this, a Generalized Least Squares (GLS) based method is used to estimate the time-dependent route flows from field data.

The step forward given in this work is that this field data is obtained from ANPR sensors which provide (1) powerful time-dependent information directly related with route flows and (2) users’ travel time information which allows the model to obtain better estimations.

The paper is organized as follows: Section 3 summarizes the variables used all throughout the paper and their meaning, and Section 4 introduces the two main tools used in the model. In Section 5, the proposed model is explained in detail to later be applied to several well-known networks in Section 6. Finally, some conclusions and future research are provided in Section 7.

3. Notation

$a$	Link.
$\mathit{A}$	Set of links.
${\mathit{A}}_{\mathit{w}}\left(\mathit{s}\right)$	Set of links associated to cumulative observed flow $W_s^a\left(t\right)$.
$D_a\left(t\right)$	Lower bound of link travel time for vehicles entering link $a$ at time $t$.
$D_a^0$	Travel time assuming the maximum allowed speed on link $a$.
$E_a\left(t\right)$	Cumulative flow associated with the exit of link $a$ by time $t$.
$e_a\left(t\right)$	Flow intensity associated with the exit of link $a$ by time $t$.
$E_r^a\left(t\right)$	Cumulative flow associated with the exit of link $a$ and route $r$ by time $t$.
$e_r^a\left(t\right)$	Flow intensity associated with the exit of link $a$ and route $r$ by time $t$.
$h_r^0\left(t\right)$	Prior or out-of-date route flow intensity at the origin of each route $r$ at time $t$.
$h_r\left(t\right)$	Route flow intensity at the origin of each route $r$ at time $t$.
$H_r^0\left(t\right)$	Prior or out-of-date cumulative flow intensity at the origin of each route $r$ at time $t$.
$H_r\left(t\right)$	Cumulative flow intensity at the origin of each route $r$ at time $t$.
$ite{r}_{CDNL}$	Iteration of the CDNL approach.
$ite{r}_{global}$	Iteration of the global approach.
$\mathit{N}$	Set of nodes.
$Q_{ak}$	Queue dissipation time.
$r$	Route.
$R$	Set of routes.
$\mathit{S}\left(a\right)$	Set of all links downstream of link $a$.
$t_a^{out}$	Exit time of the vehicle that enters the link $a$ at time $t$.
$t_{ak}^{out}$	Link exit time of a user who enters link $a$ at time $t_k$.
$t_k$	Set of discrete times to be considered for $k=1,\dots ,m$.
$tol$	Tolerance.
$T_{veh}$	Registered travel time of the vehicle.
$W_s^a\left(t\right)$	Cumulative observed flow in set $s$ recorded with the sensor installed in link $a$ up to time $t$.
$x_a\left(t\right)$	Number of vehicles (link volume) in link $a$ at time $t$.
$x_a^{max}$	Maximum number of vehicles (link capacity) in link $a$ leading to a travel time ${\beta }_a\left(1+{\gamma }_a\right)$.
${\beta }_a$	Saturation factor of link $a$.
${\gamma }_a$	Saturation exponent of link $a$.
${\delta }_a$	Saturation factor ahead of link $a$.
${\epsilon }_{CDNL}$	Control error of the CDNL approach.
${\epsilon }_{global}$	Control error of the global approach.
${\theta }_r^a\left(t\right)$	Departure time of a vehicle that follows route $r$ and left the link $a$ at time $t$.
${\lambda }_r^s$	Route-set of scanned links incidence matrix.
$\rho$	Relaxation factor.
${\tau }_{a_n}^{-1}\left(t\right)$	Entry time of a vehicle that exists in link $a$ at time $t$.
${\tau }_a^{iter}\left(t\right)$	Exit time of a vehicle from link $a$ associated with entering it at time $t$ at iteration, $iter$.
${\phi }_s^{sr}\left(t\right)$	Observed flow proportion in set $s$ that used sub-route $sr$ during the period of time from $t-\Delta t$ to $t.$
${\chi }_a^{sr}$	Link-sub-route incidence matrix.

4. Background

The proposed model is designed to apply a Continuous Dynamic Network Loading (CDNL) model based on the one proposed by Castillo et al. [22], and use data from ANPR sensors to estimate the route traffic flows (see, for example, Castillo et al. [37], Cerrone et al. [49] or Álvarez-Bazo et al. [40]). Both are explained in the following sections.

The CDNL Model

This model needs the following information:

• The network topology.
• The route cumulative flow, $H_r\left(t\right)$, or alternatively, the flow intensity, $h_r\left(t\right),$ at the origin of each route, $r$.
• The density-cost function of each link, $a$:

  $$D_a\left(t\right)=f\left({\overrightarrow{\xi }}_a,{\overrightarrow{v}}_a\left(t\right)\right)$$

  (1)

  where ${\overrightarrow{\xi }}_a$ is the parameter vector and ${\overrightarrow{v}}_a\left(t\right)$ is the variable. A BPR (Bureau of Public Roads)-based function [50] would be a possibility to use, as proposed by Castillo et al. [22].

To illustrate concepts, let us use the simple example shown in

Table 1. Network, O-D pairs and routes used as an illustrative example.

Network	O-D	$\mathbf{Route}\mathit{r}$	$\mathbf{Link}\mathit{a}$
	1–4	1	1	5	8	10
		2	2	8	10
		3	3	9	10
		4	3	6	8	10
		5	4	7	9	10
		6	4	7	6	8	10
	2–4	7	5	8	10
	3–4	8	7	6	8	10
		9	7	9	10

. It is a network $\left(\mathit{N},\mathit{A}\right)$, where $\mathit{N}$ is a set of 7 nodes, with 3 origin nodes (1, 2 and 3) and 1 destination node (4), and $\mathit{A}$ is a set of 10 links. The topology of the network leads to 3 O-D pairs, as indicated in Table 1. The table also shows a set $\mathit{R}$ of 9 routes, whose route flow intensities are assumed known. The route 6 of this illustrative example is used as an analysis case to facilitate the understanding of this model. Specifically, the upper curve of Figure 1 represents the traffic intensity over time at the origin of this route. The model allows to know how the origin route flow wave is propagated through all links of the route and provides the traffic flow intensity at the exit of each link, $a$. Note that this figure and all the figures throughout the paper represent an extended day of 30 h because some of the simulated users finish their trips after the hour 24. Actually, from 24 to 30 h, the figures represent the traffic flow from 00:00 to 6:00 of the next day, which should be mixed with those new users, but for the sake of simplicity, we have preferred to model just the users of one day.

In particular, this route is made up by links 4, 7, 6, 8 and 10. In Figure 1, the upper plot corresponds to the origin node traffic intensity, and the other five plots indicate how this wave moves due to time and deforms due to congestion. This model can be classified as dynamic and it provides results for each instant through the approximation of link exit–entry time by monotone cubic-spline functions. Therefore, this model allows to determine the position of any vehicle at any time by only knowing the moment when it began its journey through the network and vice versa. With this in mind, consider, for example, a vehicle that started the route 6 at 13:58 h. The green areas represent the number of vehicles of the corresponding route 6 which left each link before the studied vehicle. All areas below these curves are identical at all times because the FIFO condition is satisfied in this model.

Nevertheless, the shape of the curves evolves with time depending on the link congestion degree. For example, the wave of link 6 is narrower and longer. This implies that it presents a lower traffic intensity and a longer travel time due to higher congestion degree. This fact evidences the model complexity and quality.

Obviously, the model considers all routes containing link $a$ to evaluate the potential traffic problems that may arise from flow confluence. Thus, the traffic flow intensity at the exit of link $a$ is the result of the traffic intensities due to all routes containing link a. Specifically for this analysis case, the flow intensities at the exit of links 4, 7, 6 and 8 are given by adding the traffic intensities of the 3, 8, 4 and 6 routes respectively, that use each of the mentioned links.

Figure 2 illustrates the flow intensity curves for all links of route 6, where the different colors refer to the different route components. Note that the sum of route flows takes into account the delays associated with each route. Therefore, the knowledge of both the position and the origin of any vehicle is guaranteed using this continuous dynamic loading model.

Note that one of the cores of this continuous dynamic network loading traffic model is that it needs to know all routes’ flow intensities at their origins. Taking this into account, the aim of this paper is to develop a method to reconstruct these $h_r\left(t\right)$ functions from field data. The novel contribution of this paper is that, differently from the majority of papers existing in the literature, we have used the powerful data collected from ANPR sensors.

The interest in using this data for traffic flow estimation has increased due to its considerably larger amount of resultant information compared with other standard methods (see, for example, Castillo et al. [37,51,52]).

To take advantage of all the potential of the ANPR data, the set of sensors have to be optimally located at a given traffic network to reconstruct both route flows and travel times (see, for example, Mínguez et al. [38], Cerrone et al. [49] or Sánchez-Cambronero et al. [53]). The sensors automatically recognize the plate number of the passing vehicle through a scanned link and record the passing times (Álvarez-Bazo et al. [40]). In the following, we briefly explain how to use the ANPR data by using a simple example.

Consider the network topology of the illustrative example and assume a set of sensors located at the set $\mathit{S}\mathit{L}=\left\{3,5,7,10\right\}$ (marked in red in Table 1). Then, intersecting each route in $\mathit{R}$ with set $\mathit{S}\mathit{L}$, one can build the “scanning map” (see Castillo et al. [52]), to obtain the Sets of Combinations of Scanned Links ($\mathit{SCSL}$), where flow is expected to be directly monitored (i.e., observed). For this particular case, the components of this set are as follows:

• $s_1$ = {5, 10}, where the vehicles that transverse links 5 and 10 can be computed;
• $s_2$ = {10}, where the vehicles that use only in link 10 can be computed;
• $s_3$ = {3, 10}, where the vehicles to be scanned in links 3 and 10 can be computed;
• $s_4$ = {7, 10}, where the vehicles to be scanned only in links 7 and 10 can be computed.

Then, with the route information of Table 1, the relations between these observed flows and the desired route flows can be derived by means of the matrix $\Lambda$, where ${\mathsf{\lambda }}_r^s=1$ if the flow of route $r$ is contained in the observed flow related to the set $s$ of $\mathit{SCSL}$, and 0 otherwise, i.e.,

$$W_s^a\left(t\right)={\displaystyle \sum }_r{\lambda }_r^s{E}_r^a\left(t\right)$$

(2)

where $W_s^a\left(t\right)$ is the cumulative observed flow in set $s$ recorded with the sensor installed in link $a$ up to time $t$. For example, let $W_{s_1}^{a_5}\left(t\right)$ be the cumulative observed flow in $s_1$ recorded with the sensor installed in link 5, which is the sum of the cumulative flows $E_{r_1}^{a_5}\left(t\right)$ and $E_{r_7}^{a_5}\left(t\right)$ of routes $r_1$ and $r_7$, respectively. Similarly, $W_{s_2}^{a_{10}}\left(t\right)$, i.e., the observed flow in $s_2,$ will directly be the cumulative flow $E_{r_2}^{a_5}\left(t\right)$ of route $r_2.$

Table 2. Scanning map with the Set of Combinations of Scanned Links ($\mathit{SCSL}$) and the components of $\Lambda$ matrix for the illustrative example.

Scanning Map					Components of Λ Matrix
$\mathit{S}\mathit{C}\mathit{S}\mathit{L}$	$\mathbf{Links}\mathbf{in}\mathit{S}\mathit{L}$				$\mathbf{Routes}\mathbf{in}\mathit{R}$
	${\mathit{a}}_{\mathbf{3}}$	${\mathit{a}}_{\mathbf{5}}$	${\mathit{a}}_{\mathbf{7}}$	${\mathit{a}}_{\mathbf{10}}$	${\mathit{r}}_{\mathbf{1}}$	${\mathit{r}}_{\mathbf{2}}$	${\mathit{r}}_{\mathbf{3}}$	${\mathit{r}}_{\mathbf{4}}$	${\mathit{r}}_{\mathbf{5}}$	${\mathit{r}}_{\mathbf{6}}$	${\mathit{r}}_{\mathbf{7}}$	${\mathit{r}}_{\mathbf{8}}$	${\mathit{r}}_{\mathbf{9}}$
$s_1$		X		X	1	0	0	0	0	0	1	0	0
$s_2$				X	0	1	0	0	0	0	0	0	0
$s_3$	X			X	0	0	1	1	0	0	0	0	0
$s_4$			X	X	0	0	0	0	1	1	0	1	1

shows this information for the rest of the relations and Figure 3 shows a graphical representation of all the cumulative observed flows, $W_s^{\mathrm{a}}\left(t\right)$.

Note that one $W_s^a\left(t\right)$ for each link of each $s$ can be obtained. However, for computational purposes (see Section 4), it is convenient to define a new set of links, ${\mathit{A}}_{\mathit{w}}\left(s\right)\subset \mathit{SL}$, such that only one flow be computed for each $s.$ In our case (see Figure 3), these links are $a_5$ for $s_1$, $a_{10}$ for $s_2$, $a_3$ for $s_3$ and $a_7$ for $s_4$.

So far, most of the existing models for traffic flow estimation that use plate scanning data are related to static flow estimation. However, since this set of sensors also supplies very valuable travel time information, such as those included in

Table 3. Extract of the registered data by the sensors located in $\mathit{SL}=\left\{3,5,7,10\right\}$ (the plate numbers are simulated and cannot be related to any vehicle).

Item	Plate #	Time Stamps in Each Link				Travel Time (Tveh)	SCSL
		3	5	7	10
…	…	…	…	…	…	…	…
13,548	1226AQU			09:59:59	14:10:09	04:10:11	$s_4$
13,549	3332AQO	10:00:04			11:22:11	01:22:08	$s_3$
13,550	7539UII		10:00:05		11:23:04	01:22:59	$s_1$
13,551	4095EEA			10:00:05	11:22:12	01:22:08	$s_4$
13,552	7558IEE				10:00:07		$s_2$
13,553	9208QAU			10:00:08	14:10:19	04:10:11	$s_4$
13,554	6179QZO			10:00:08	11:22:16	01:22:08	$s_4$
13,555	8266IAU		10:00:11		11:23:10	01:22:59	$s_1$
13,556	2349QQE	10:00:11			14:10:23	04:10:11	$s_3$
13,557	0195UIZ			10:00:17	14:10:29	04:10:12	$s_4$
13,558	1840UEI			10:00:18	11:22:26	01:22:07	$s_4$
13,559	5891QIQ	10:00:19			11:22:26	01:22:07	$s_3$
13,560	6282OUA			10:00:19	14:10:31	04:10:12	$s_4$
13,561	2272QIZ				10:00:21		$s_2$
13,562	7420UAQ			10:00:21	11:22:28	01:22:07	$s_4$
13,563	4068AIA		10:00:23		11:23:23	01:22:59	$s_1$
13,564	6073UQU		10:00:23		11:23:23	01:22:59	$s_1$
13,565	7390UEA			10:00:26	14:10:38	04:10:12	$s_4$
13,566	2716UAI	10:00:28			14:10:40	04:10:12	$s_3$
13,567	0959IAU			10:00:32	11:22:39	01:22:07	$s_4$
13,568	2977IAI			10:00:33	11:22:40	01:22:07	$s_4$
…	…	…	…	…	…	…	…

, we propose to also use it for dynamic flow estimation. In addition to the information described above, vehicles’ sub-route travel time can be derived since a single user can be identified in different scanned links at different times (Sánchez-Cambronero et al. [54]; Gentili and Mirchandani [55]).

For example, the vehicle with plate number 3332AQO has been registered when it left the link 3 at time 10:00:04 and the link 10 at time 11:22:11, meaning that the travel time ($T_{veh}$) for this particular vehicle from one stamp to the other was of 1 h, 22 min and 8 s. The vehicle with plate number 2349QQE was also registered in the same $\mathit{SCSL}{s}_3=\left\{3,10\right\}$, but the time stamps of this vehicle indicate that its travel time was 4 h, 10 min and 11 s, which means that despite the sum in the same observed flow (i.e., $W_{s_3}^{a_3}\left(t\right)=E_{r_3}^{a_3}\left(t\right)+E_{r_4}^{a_3}\left(t\right)$), it is probable that they do not share the same route, hence more information is needed to derive the cumulative route flows. In particular, one can write:

$${\displaystyle \sum }_r{\sigma }_{r,a}^{sr}\left(E_r^a\left(t\right)-E_r^a\left(t-\mathsf{\Delta }t\right)\right)={\displaystyle \sum }_s{\phi }_s^{sr}\left(t\right)\left(W_s^a\left(t\right)-W_s^a\left(t-\mathsf{\Delta }t\right)\right)$$

(3)

where ${\sigma }_{r,a}^{sr}$ = 1 if the sub-route $sr$ is included in route $r,$ and 0 otherwise. ${\phi }_s^{sr}\left(t\right)$ is, for the interval $\mathsf{\Delta }t$ before time $t,$ the proportion of the observed flow in set $s$ that have used sub-route $sr$. We have to note here that we consider a sub-route as a sequence of consecutive links $\left(a_1,a_{2,}\dots ,a_n\right)$ included in at least one route of $\mathit{R}$, such that $a_1$ and $a_n$ are included in $\mathit{SL}$.

Table 4. Sub-route links and components of $\Gamma$ and $\Phi$ matrices.

					Components of Γ									Components of Φ Matrix
					$\mathbf{Routes}\mathbf{in}\mathit{R}$									$\mathit{S}\mathit{C}\mathit{S}\mathit{L}$
Sub-Route	Links				${\mathit{r}}_1$	${\mathit{r}}_2$	${\mathit{r}}_3$	${\mathit{r}}_4$	${\mathit{r}}_5$	${\mathit{r}}_6$	${\mathit{r}}_7$	${\mathit{r}}_8$	${\mathit{r}}_9$	${\mathit{s}}_1$	${\mathit{s}}_2$	${\mathit{s}}_3$	${\mathit{s}}_4$
$s{r}_1$	3	9	10		0	0	1	0	0	0	0	0	0	0	0	${\phi }_{s_3}^{s{r}_2}\left(t\right)$	0
$s{r}_2$	3	6	8	10	0	0	0	1	0	0	0	0	0	0	0	${\phi }_{s_3}^{s{r}_3}\left(t\right)$	0
$s{r}_3$	7	9	10		0	0	0	0	1	0	0	0	1	0	0	0	${\phi }_{s_4}^{s{r}_4}\left(t\right)$
$s{r}_4$	7	6	8	10	0	0	0	0	0	1	0	1	0	0	0	0	${\phi }_{s_4}^{s{r}_5}\left(t\right)$

shows matrices $\Gamma$ and $\Phi$, and allows to conclude that with this information, the estimation of almost all routes can be greatly improved. In particular, since travel times of sub-routes $s{r}_2$ and $s{r}_3$ are very different, Equation (3) allows to have an accurate estimation to distinguish the users of route 3 from those of route 4.

We propose to obtain the parameters of matrix $\Phi$ as a function of $P\left(sr|T_{veh}\right),$ i.e., the probability of each user to choose the sub-route $sr$ knowing that its travel time was $T_{veh}$. Thus, using the Bayes’ theorem, we have:

$$P\left(sr|T_{veh}\right)=\frac{P\left(sr,t\right)P(T_{veh}|sr)}{{{\displaystyle \sum }}_{i}P\left(s{r}_i,t\right)P(T_{veh}|s{r}_i)}$$

(4)

On one hand, $P\left(sr, t\right)$ is the prior probability of the vehicle (which enters to the sub-route by time t) to use the sub-route $sr$. We propose to obtain this probability using the logit formula with the prior information, i.e.,

$$P\left(sr,t\right)=\frac{ex{p}^{-\theta C_{sr}\left(t\right)}}{{{\displaystyle \sum }}_{i}ex{p}^{-\theta C_{s{r}_i}\left(t\right)}}$$

(5)

where $C_{sr}\left(t\right)$ is assumed to be the prior travel time of sub-route $sr$ by time $t$ (i.e., the user’s entry time to the sub-route). Note that to obtain this travel time, the continuous dynamic loading model used in this paper is very useful since it allows to reproduce the travel time of each single user no matter the link entry time nor the used route.

On the other hand, $P(T_{veh}|sr)$ is the likelihood that the vehicle whose travel time was $T_{veh}$ followed sub-route $sr.$ We assume that the greater the difference between the measured travel time ($T_{veh})$ and the estimated time for each user, the more likely it is that it had not used the sub-route. Therefore, the probability that the vehicle whose travel time was $T_{veh}$ did not follow the sub-route $sr$ is:

$$P(\overline{T_{veh}|sr})=\frac{\left|T_{veh}-{{\displaystyle \sum }}_a{\chi }_a^{sr}{D}_a\left(t_{veh}^a\right)\right|}{{{\displaystyle \sum }}_i\left|T_{veh}-{{\displaystyle \sum }}_a{\chi }_a^i{D}_a\left(t_{veh}^a\right)\right|}$$

(6)

where $T_{veh}$ is the registered travel time of the vehicle and ${\chi }_a^{sr}$ is 1 if link $a$ is included in sub-route $sr$, and 0 otherwise. Hence, ${\displaystyle \sum }_a{\chi }_a^{sr}{D}_a\left(t_{veh}^a\right)$ is the estimated travel time using any expression of (1), or as we will see in the next section, a spline interpolation. Then, since a vehicle only can choose one sub-route:

$$P\left(T_{veh}|s{r}_q\right)={\displaystyle \sum }_{i\ne q}P(\overline{T_{veh}|s{r}_i}).$$

(7)

The rationale behind Expression (4) is that we must assign each vehicle to the sub-route which estimated travel time as closer to the actual measured time using the ANPR data. Then, ${\phi }_s^{sr}\left(t\right)$ can be obtained as:

$${\phi }_s^{sr}\left(t\right)=\frac{{{\displaystyle \sum }}_{veh|t_{veh}\in \mathsf{\Delta }{t}_n}P\left(sr|T_{veh}\right)}{W_s^a\left(t\right)-W_s^a\left(t-\mathsf{\Delta }t\right)}$$

(8)

where $W_s^a\left(t\right)-W_s^a\left(t-\mathsf{\Delta }t\right)$ is the number of observed vehicles in the set $s$ during the period $t-\mathsf{\Delta }t$.

To sum up, the use of ANPR data as shown in Figure 3 implies knowing only some $W_s^{\mathrm{a}}\left(t\right)$ cumulative flow (or the corresponding flow intensity if we derive $W_s^{\mathrm{a}}\left(t\right)$) at the end of the scanned links. However, the basic data for traditional dynamic traffic assignment, as proposed by Castillo et al. [22], needs the route flow intensities at their origins. For this reason, the core of the proposed model is to reconstruct the route flow intensity at the origin of each route in such a way that once those route flows are assigned to the network, the results represent the already observed flow as much as possible.

5. The Proposed Model for Dynamic Route Flow Estimation Using Data from ANPR Sensors

In this section, the dynamic route flow estimation is developed considering all the background described above and based on some assumptions.

5.1. Model Assumptions

Assumption 1.

The proposed model is proposed to be a FIFO (first-in-first-out) consistent model. This means that users who enter the link earlier will leave it sooner. This assumption is very common in most dynamic traffic models since it avoids some computing problems, and in macroscopic modelling, FIFO is considered a better approximation to reality than allowing arbitrary deviations from FIFO (Peeta and Ziliaskopoulos [56]; Carey et al. [57]).

Assumption 2.

The link transversal time$D_a\left(t\right)$of a vehicle that enters link$a$at instant$t$needed to propagate the route flow, is assumed to be a function of the number of vehicles$x_a\left(t\right)$(link volume) in link$a$at instant$t$and, to reproduce the shock wave due to congestion, it also depends on the immediate downstream link volumes. In particular, the BPR function has been adapted in this paper as follows (actually,$D_a\left(t\right)$is only given when no queue exists. To deal with queue and respect the FIFO rule, some modifications have to be done; in particular, adding a minimum delay equivalent to considering a queue of vehicles to congestion if it exists (see the algorithm)):

$$D_a\left(t\right)=D_a^0\left[1+{\beta }_a{\left(\frac{x_a\left(t\right)}{x_a^{max}}\right)}^{{\gamma }_a}+{\delta }_a\underset{b\in \mathit{S}\left(a\right)}{\max }{\left(\frac{x_b\left(t\right)}{x_b^{max}}\right)}^{{\gamma }_b}\right]$$

(9)

where$D_a^0$is the travel time assuming the maximum allowed speed,${\beta }_a$and${\gamma }_a$are the saturation parameters,${\delta }_a$is a parameter used to take into account the congestion ahead of link$a$and$\mathit{S}\left(a\right)$is the set of all links downstream of link$a$. Regarding the number of vehicles on link$a$at time$t$, $x_a\left(t\right)$, it is calculated using the cumulative flow at the end of the link$a$until time$t$$(E_a\left(t\right))$and the cumulative flow at the end of link$a$until the time$t_a^{out}$ ($E_a\left(t_a^{out}\right)$), which is the exit time of the vehicle that enters the link at time $t$ (that can be computed using (9)). That is:

$$x_a\left(t\right)=E_a\left(t_a^{out}\right)-E_a\left(t\right)$$

(10)

A shown in Figure 1, this number of vehicles can be obtained looking backwards in the cumulative route flow $H_r\left(t\right)$ at the origin of each route $r$ that includes link $a$, i.e.,

$$E_a\left(t\right)={\displaystyle \sum }_r{E}_r^a\left(t\right)={\displaystyle \sum }_r{H}_r\left({\theta }_r^a\left(t\right)\right)$$

(11)

where ${\theta }_r^a\left(t\right)$ is the departure time of a vehicle that follows route $r$ and left the link $a$ at time $t$. Then,

$${\theta }_r^a\left(t\right)={\tau }_{a_n}^{-1}\left({\tau }_{a_{n-1}}^{-1}\left(\dots {\tau }_{a_1}^{-1}\left(t\right)\right)\right)$$

(12)

where $a_1,a_{2,}\dots a_n$ are the links of route $r$ and ${\tau }_a^{-1}$(t) is the entry time of a vehicle that exists at link $a$ at time $t$. As in Castillo et al. [22], we propose to use monotone cubic splines approximations to ${\tau }_a^{-1}$ as a pair of discrete entry–exit times ($t, t_a^{out})$, i.e.,

$${\tau }_a^{-1}\left(t\right)\approx splin{e}_a\left(t;\left(t_a^{out}, t\right)\right)$$

(13)

Note that knowing these splines and the observed $W_s^a\left(t\right),$ it is possible to reconstruct the route flow intensity $h_r\left(t\right)$ at the origin of each route using the proposed algorithm explained in the next section. In addition, the use of monotone cubic splines in this paper is justified because of two reasons:

• The link travel time has to be a continuous function to obtain the departure time of each identified vehicle in the ANPR sensor. Using a discretized function, this cannot be done.
• To guarantee that the FIFO rule holds, the link travel time function has to be monotonously increasing, and this spline guarantees that (see the Matlab manual for pchip function).

Note that since ${\tau }_a^{-1}\left(t\right)$ is a function of $x_a\left(t\right)$ and $x_a\left(t\right)$ is a function of ${\tau }_a^{-1}\left(t\right)$, an iterative process is needed to correctly load all the route flows through the network.

Assumption 3.

The ANPR sensors located in links of set$\mathit{SL}$are going to be installed at the end of the links. This assumption implies that the data registered by the sensors, which allows to obtain the observed flow$W_s^a\left(t\right),$is compatible with all those stated in Assumption 2. In particular, to reconstruct the route flow intensity$h_r\left(t\right)$at the origin of each route using Equation (13).

Assumption 4.

The exit time of a link coincides with the entry time of the next one of the same route. This implies that the time to transverse a node is negligible compared with the link travel time.

Assumption 5.

Any error in coding the vehicles’ plates may affect the quality of the traffic flow estimation results and also obtaining travel times in the network. However, dealing with this problem is a complex task that deserves a deeper analysis, and it falls beyond the scope of this work. For that reason, we will assume error-free data in the proposed model. In any case, the high number of reasons that errors could happen can be classified in two groups: those that come from the data collection phase and those that happen due to the data processing phase. The strategy to address these errors is different. In particular:

• During the data collection, one can find a high number of problems which may suggest that the used OCR (Optical Character Recognition) technique leads to a wrong or incomplete plate. Among these, we can name a few: high speed of the vehicles, a vehicle hiding the plate of another, blurry image, image with high brightness due to the sun position, the vehicle is far from the camera so the plate cannot be identified, the angle of the vehicle with respect to the sensor made the plate unreadable, etc. These problems can be mitigated using a sensor of high quality, installing replicated sensors (i.e., several sensors which register the same traffic flow) and/or with a detailed study of the sensor position in the scanned link or avoiding links where those problems are expected. In short, trying to minimize the probability of error in the data-collection phase (see, for example, Castillo et al. [58] and Álvarez-Bazo et al. [40], where the authors worked in this direction).
• Once the data are collected, if the plate number is wrong, the user cannot be correctly identified passing for the links of its$\mathit{SCSL}$(the same user may appear in several rows of a dataset, such as the one shown in Table 3). Therefore, this user will be assigned to a wrong$\mathit{SCSL}$(or even to several$\mathit{SCSL}$) and hence, computed as a user of another route (see, for example, Sánchez-Cambronero et al. [59] or Zhu et al. [60]). To deal with this problem, specific models are needed were data besides the license plate (such as, for example, the vehicle’s manufacturer, model, color, etc.) can be useful to assign each user to the correct$\mathit{SCSL}$.

5.2. The Proposed Algorithm

Figure 4 shows the complete flowchart of the proposed algorithm, which will be explained in detail in this section.

INPUT: The network topology, the set of intervals to interpolate the splines $t_k \left(k=1,2,\dots , m\right)$, the free-flow travel time $D_a^0$, the link parameter $x_a^{max}$, the link saturation parameters ${\beta }_a$, ${\gamma }_a$ and ${\delta }_a$, the tolerance to stop the process, $tol$, the maximum allowed iterations $itermax$, the set of scanned links $\mathit{SL}$, the cumulative observed flow $W_s^a\left(t\right)$ and prior or out-of-date route flow intensity $h_r^0\left(t\right)$.

OUTPUT: All the estimated dynamic traffic flow information: link travel time by means of the set of spline functions ${\tau }_a^{-1}\left(t\right)$, route flow intensity $h_r\left(t\right)$, cumulative route flow $H_r\left(t\right),$ link flow intensity $e_a\left(t\right)$, link cumulative flow $E_a\left(t\right)$, etc.

Step 0: Initialization.

To start the algorithm, we propose to load a prior or out-of-date route flow intensity $h_r^0\left(t\right)$ which is going to be updated in each iteration of the algorithm. Therefore, we obtain the cumulative flow $H_r^0\left(t\right)$ by using Equation (14):

$$H_r^0\left(t\right)={{\displaystyle \int }}_0^t{h}_r^0\left(t\right)dt$$

(14)

The essential tool to propagate the cumulative flow is to know $D_a\left(t\right)$ to then interpolate the basic points $t_k,$ to obtain ${\tau }_a^{-1}\left(t\right).$ Since initially the number of vehicles $x_a\left(t\right)$ to compute Equation (9) is not known, we need to make an initial guess so as to later iterate until convergence. So, let the initial set of points for the spline function (13) be $\left(t_{ak}^{out},t_k\right)$, where $t_{ak}^{out}=t_k+D_a^0$. Also, set:

• ${\epsilon }_{global}=2tol$ to control the differences in the results between two consecutive iterations of the global approach.
• $ite{r}_{global}=1$ to start the iteration process of the global approach.
• ${\epsilon }_{CDNL}=2tol$ to control the differences in the results between two consecutive iterations of the CDNL approach.
• $ite{r}_{CDNL}=1$ to start the iteration process of the CDNL approach.

Step 1: Obtain ${\tau }_a^{-1}\left(t\right) \forall a\in \mathit{A}$ by applying the Continuous Dynamic Network Loading (CDNL) model.

In this step, the current estimation of the cumulative route flow $H_r^0\left(t\right)$ is loaded to the network. Therefore, an iterative process is used as follows:

Step 1.1: Check the continuity or stop the CDNL iterative process. While ${\epsilon }_{CDNL}>tol$ and $ite{r}_{CDNL}$< $itermax$, go to Step 1.2. Otherwise, go to Step 2.

Step 1.2: Calculate $D_a\left(t_k\right)$ in all links of the network. The aim of this step is to obtain the lower bound of each link travel time in each discrete time $t_k$ using Equation (9). Then, for $\mathrm{a}=1$ to $\mathrm{a}={\mathrm{n}}_{\mathrm{links}}$, obtain:

$$D_a\left(t_k\right)=D_a^0\left[1+{\beta }_a{\left(\frac{x_a\left(t_k\right)}{x_a^{max}}\right)}^{{\gamma }_a}+{\delta }_a\underset{b\in \mathit{S}\left(a\right)}{\max }{\left(\frac{x_b\left(t_k\right)}{x_b^{max}}\right)}^{{\gamma }_b}\right],$$

(15)

where

$$x_a\left(t_k\right)=E_a\left(t_{ak}^{out}\right)-E_a\left(t_k\right)$$

(16)

and using the current spline function (obtained in the current $ite{r}_{CDNL}$) and Equations (11) and (12):

$$\begin{array}{c}{E}_a\left(t_k\right)={\displaystyle \sum _r}{E}_r^a\left(t_k\right)\\ ={\displaystyle \sum _r}{H}_r^0\left({\theta }_r^a\left(t_k\right)\right)| {\theta }_r^a\left(t_k\right)={\tau }_{a_n}^{-1}\left({\tau }_{a_{n-1}}^{-1}\left(\dots {\tau }_{a_1}^{-1}\left(t_k\right)\right)\right)\end{array}$$

(17)

$$\begin{array}{c}{E}_a\left(t_{ak}^{out}\right)={\displaystyle \sum _r}{E}_r^a\left(t_{ak}^{out}\right)\\ ={\displaystyle \sum _r}{H}_r^0\left({\theta }_r^a\left(t_{ak}^{out}\right)\right)| {\theta }_r^a\left(t_{ak}^{out}\right)={\tau }_{a_n}^{-1}\left({\tau }_{a_{n-1}}^{-1}\left(\dots {\tau }_{a_1}^{-1}\left(t_{ak}^{out}\right)\right)\right)\end{array}$$

(18)

Step 1.3: Obtain the new $t_{ak}^{out}$ taking into account the existence of queues. Since $D_a\left(t_k\right)$ is a lower bound of the link travel time of a user that enters to link $a$ at time $t$, we have that $t_{ak}^{out}\ge t_k+D_a\left(t_k\right)$. However, if a queue exists in the link, we must take into account the queue dissipation time $Q_{ak}$. To calculate if a queue exists or not between $t_{ak}^{out}$ and $t_{ak-1}^{out}$, we have to know the link flow that used the link in that period of time, i.e., $E_a(t_{ak}^{out})-E_a(t_{ak-1}^{out})$. Assuming that if $x_a^{max}$ vehicles flow in a time of $D_a^0(1+{\beta }_a$) (i.e., when $x_a\left(t\right)/x_a^{max}=1$), then, using the proportionality rule,

$$Q_{ak}=[E_a(t_{ak}^{out})-E_a(t_{ak-1}^{out})]\frac{D_a^0(1+{\beta }_a)}{x_a^{max}}.$$

(19)

Therefore, to be consistent with the FIFO rule (Assumption 1), we propose to obtain the $t_{ak}^{out}$ values to be used to update the splines (13) as follows:

$$t_{ak}^{out}=max\{t_k+D_a\left(t_k\right),t_{ak-1}^{out}+Q_{ak}\}$$

(20)

Step 1.4: Update the splines using the basic exit–entry times $(t_{ak}^{out}, t_k)$. The splines are updated using Equation (13) to calculate the current ${\epsilon }_{CDNL}$ by means of the data provided by the updated splines and those obtained in the previous iteration, i.e.,

$${\epsilon }_{DNL}={\displaystyle \sum }_{a,t}\left|{\tau }_a^{iter}\left(t_k\right)-{\tau }_a^{iter-1}\left(t_k\right)\right|$$

(21)

Do $ite{r}_{CDNL}=ite{r}_{CDNL}+1$ and go to Step 1.1.

Step 2: Update ${\widehat{E}}_a^r\left(t_k\right)$ using the splines obtained in the CDNL model of Step 1.

Since one of the basic inputs of the model to be solved in Step 3 is a reference cumulative flow associated with the exit of scanned link $a\in {\mathit{A}}_{\mathit{w}}\left(\mathit{s}\right)$ and route $r$ by each time $t_k$, in this step, the current version of the links’ travel times (the splines) are used to obtain this flow as follows:

$${\widehat{E}}_r^a\left(t_k\right)=H_r^0\left({\theta }_r^a(t_k)\right).$$

(22)

We have to note here that not all ${\widehat{E}}_a^r\left(t_k\right)$ have to be calculated, but only those such that link $a$ belongs to set ${\mathit{A}}_{\mathit{w}}\left(\mathit{s}\right),$ which is a subset of $\mathit{SL}$. Go to step 3.

Step 3: Reconstruct the cumulative route flow $H_r\left(t\right)$ using the observed data.

It is well-known that the most widely used methods of O-D trip matrices estimation from measured traffic data (link counts, ANPR data, etc.) are based on mathematical programming techniques, such as least squares methods or entropy-based methods, among others. In this paper, we propose to use a least squares-based model to reconstruct the cumulative route flow at the beginning of each route, again using the splines but this time from $E_r^a\left(t_k\right)$ flows.

Step 3.1: Solve a least squares-based optimization approach. To obtain $E_r^a\left(t_k\right)$ from the observed flow $W_s^a\left(t_k\right)$, solve the following optimization problem:

$$\underset{E}{\min } z={\displaystyle \sum }_{r,t_k,a}{\left(\frac{{\widehat{E}}_r^a\left(t_k\right)-E_r^a\left(t_k\right)}{G_r^a\left(t_k\right)}\right)}^2$$

(23)

Subject to:

$$W_s^a\left(t_k\right)={\displaystyle \sum }_r{\lambda }_r^s{E}_r^a\left(t_k\right) \forall s; \forall t_k; \forall a\in {\mathit{A}}_{\mathit{w}}\left(\mathit{s}\right)$$

(24)

$$\begin{array}{c}{\displaystyle \sum _r}{\sigma }_{r,a}^{sr}\left(E_r^a\left(t_k\right)-E_r^a\left(t_{k-1}\right)\right)={\displaystyle \sum _r}{\phi }_s^{sr}\left(t_k\right)\left(W_s^a\left(t_k\right)-W_s^a\left(t_{k-1}\right)\right) \\ \forall sr; \forall t_k|k\rangle 1; \forall a\in {\mathit{A}}_{\mathit{w}}\left(\mathit{s}\right) \end{array}$$

(25)

$$E_r^a\left(t_k\right)\ge E_r^a\left(t_{k-1}\right) \forall s; \forall t_k|k\rangle 1; \forall a\in {\mathit{A}}_{\mathit{w}}\left(\mathit{s}\right)$$

(26)

Note that Equations (24) and (25) are the already described Equations (2) and (3) but using the discrete times $t_k$. Equation (26) is an additional constraint to ensure that the estimated cumulative flow is monotonously increasing. Note also that defining Equations (24) and (25) only for $a\in {\mathit{A}}_{\mathit{w}}\left(\mathit{s}\right)$, both the number of variables and equations are reduced considerably. Note also that since some ${\widehat{E}}_r^a\left(t_k\right)$ may take values of 0 (mainly for the first values of $t_k$), we have defined $G_r^a\left(t_k\right)={\widehat{E}}_r^a\left(t_k\right)$ if ${\widehat{E}}_r^a\left(t_k\right)>0.01,$ and $1$ otherwise.

Step 3.2: Update the cumulative flow $H_r\left(t\right)$. Use the current version of the splines (${\tau }_a^{iter}\left(t\right)$) to obtain the departure times (${\theta }_r^a\left(t_k\right))$. With them and the solutions of the model (23)–(26), the cumulative flow at the origin of each route can be updated as follows:

$$H_r\left({\theta }_r^a\left(t_k\right)\right)=E_r^a\left(t_k\right)$$

(27)

were ${\theta }_r^a\left(t_k\right)$ is obtained, using Equation (12), i.e.,

$${\theta }_r^a\left(t_k\right)={\tau }_{a_n}^{-1}\left({\tau }_{a_{n-1}}^{-1}\left(\dots {\tau }_{a_1}^{-1}\left(t_k\right)\right)\right)$$

(28)

Since the dynamic loading traffic model used in this paper is a continuous model, we have used a spline function to interpolate the flow using the discrete pairs ($H_r\left({\theta }_r^a\left(t_k\right)\right), {\theta }_r^a\left(t_k\right))$. After that, go to Step 4.

Step 4: Check convergence.

This step has two basic sub-steps:

Step 4.1: Update ${\epsilon }_{global}$. The ${\epsilon }_{global}$ is updated as the absolute error between the current $H_r\left(t_k\right)$ and the used $H_r^0\left(t_k\right)$ to load the network in Step 1, i.e.,

$${\epsilon }_{global}={\displaystyle \sum }_{r,t}\left|H_r^0\left(t_k\right)-H_r\left(t_k\right)\right|$$

(29)

Go to step 4.2.

Step 4.2: End of the algorithm. If the convergence criteria have been met (${\epsilon }_{global}<tol$ or $ite{r}_{global}=itermax$), return the current splines as the final splines and $H_r\left(t\right)$ as the cumulative inflow rate at the origin of route $r$ at time $t$, which are the tools to evaluate all the flow functions. Otherwise, use the following update equation:

$$H_r^0\left(t_k\right)=\rho H_0\left(t_k\right)+\left(1-\rho \right)H_r^0\left(t_k\right)$$

(30)

and go to Step 1. Note that we propose to use the relaxation factor $\rho$ in Equation (30) to avoid oscillations so that the algorithm converges monotonically to the final solution.

6. Examples of Application

In this section, the proposed model is applied to a small illustrative example to better understand its performance and to two other more complicated examples to prove that it is appropriate and useful in real practice.

6.1. The Simple Illustrative Example

In this section, the illustrative network of Table 1 is used. It consists of 7 nodes, 10 links, 3 O-D pairs and 9 routes. Let us assume that the traffic intensity of each route is modelled as a linear combination of normal densities, where the total area down the curve is the total route flow (in fact, the shape of these intensity curves can be represented with sufficient precision, as proposed in Equation (31)). In any case, this function can be replaced by others without any problem, i.e.,

$$h_r\left(t\right)={\displaystyle \sum }_{i=1}^m{h}_{r_i}f{N}_{\left({\mu }_{r_i},{\sigma }_{r_i}\right)}\left(t\right)$$

(31)

where $m$ is the number of waves (normally two), $f{N}_{\left({\mu }_{r_i},{\sigma }_{r_i}\right)}\left(t\right)$ is the pdf (probability density function) of the normal distribution, $h_{r_i}$ is the daily flow associated with the wave $i$ and ${\mu }_{r_i}$ are the times associated with the maximum intensity and the corresponding ${\sigma }_{r_i}$ standard deviation, which measures the traffic spread of each wave. The left part of

Table 5. Parameters used to perform the illustrative example.

Parameters Used to Obtain the Route Flow Intensity							Parameters Used in the Density-Cost Function (9)
Route	${\mathit{h}}_{{\mathit{r}}_{\mathbf{1}}}$	${\mathit{\mu }}_{{\mathit{r}}_{\mathbf{1}}}$	${\mathit{\sigma }}_{{\mathit{r}}_{\mathbf{1}}}$	${\mathit{h}}_{{\mathit{r}}_{\mathbf{2}}}$	${\mathit{\mu }}_{{\mathit{r}}_{\mathbf{2}}}$	${\mathit{\sigma }}_{{\mathit{r}}_{\mathbf{2}}}$	Link	Length (km)	${\mathit{D}}_{\mathit{a}}^{\mathbf{0}}$	${\mathit{\beta }}_{\mathit{a}}$	${\mathit{\gamma }}_{\mathit{a}}$	${\mathit{\delta }}_{\mathit{a}}$	${\mathit{x}}_{\mathit{a}}^{\mathit{m}\mathit{a}\mathit{x}}$ (veh)
1	1000	9.50	2.00	1300	16.50	1.50	1	111.80	0.93	1.00	2.00	0.33	7826.24
2	1900	10.00	1.80	900	18.00	1.50	2	111.80	0.93	1.00	2.00	0.33	7826.24
3	1700	9.00	1.90	800	15.00	2.00	3	111.80	0.93	1.00	2.00	0.33	7826.24
4	1380	7.90	1.90	1500	16.60	1.70	4	111.80	0.93	1.00	2.00	0.33	7826.24
5	1600	7.50	1.90	1100	16.00	2.30	5	70.71	0.59	1.00	2.00	0.33	4949.75
6	1300	9.00	1.50	800	16.00	1.80	6	100.00	1.67	1.00	2.00	0.33	2500.00
7	800	9.00	1.70	1600	15.50	2.20	7	70.71	0.59	1.00	2.00	0.33	4949.75
8	1100	9.00	1.80	1180	16.50	2.00	8	111.80	0.93	1.00	2.00	0.33	7826.24
9	1700	8.50	1.80	1800	16.60	2.40	9	111.80	0.93	1.00	2.00	0.33	7826.24
							10	50.00	0.42	1.00	2.00	0.33	10,000.00

shows the parameters used to represent both route flow intensity $h_r^0\left(t\right)$ and cumulative route inflow $H_r^0\left(t\right)$ as the out-of-date information needed as input in the proposed algorithm. Figure 5 shows these curves. Finally, the right part of Table 5 shows the parameters used in the density-cost function (15).

Let us assume again that 4 ANPR sensors were installed in the network in set $\mathit{SL}=\left\{3,5,7,10\right\}$, and to reproduce the observed data, the following procedure was adopted:

• The assumed “true” origin route flow intensities were obtained using the following parameters to compute Equation (31):

  ∘

  $h_{r_i}=Uniform\left(800, 2500\right)$.

  ∘

  ${\mu }_{r_1}=Uniform\left(6.5, 7.5\right)$; ${\mu }_{r_2}=Uniform \left(13.00, 14.00\right)$.

  ∘

  ${\sigma }_{r_i}=Uniform\left(1, 2.5\right)$.

• The network was loaded using the model described in Step 1 of the proposed algorithm.
• Each scanned vehicle was simulated using the previous results to obtain the curves shown in Figure 3.

Finally, an extended day of 30 h has been imposed to be partitioned into 151 disjoint and exhaustive $t_k, t_{k+1}$, etc., to obtain intervals of 12 min (0.2 h). We have set $\rho =0.2$, $tol=0.01$ and $itermax=30$ to initialize the algorithm.

Figure 6 shows the model results of the route flow intensity curves at all the links of the illustrative network. To represent the contribution of each route to the total link flow intensity in a better way, we have accumulated the flow intensities of each route, so each color corresponds to a unique route. In particular, the filled colored areas and the continuous lines represent the assumed “true flows” that we try to reproduce from ANPR data. Furthermore, the different lines represent the estimations at different iterations of the algorithm (in this case, iterations 1, 5 and the final estimation) showing the good performance of the approach.

To better analyze the estimation quality of each route independently, Figure 7 presents the route flow evolution at several iterations for routes 2, 4, 7 and 8. It shows very good estimation results even in congested links.

If the results of this figure are compared with the field data provided in Figure 3, some interesting conclusions can be drawn. For example, the estimation of the dynamic flow of route 2 converges to the assumed real flow since the observed $W_{s_2}^{a_{10}} \left(t\right)$ is directly the cumulative route flow $E_{r_2}^{a_{10}}\left(t\right)$. On the contrary, despite being very good estimations, the rest of them cannot converge exactly to the real flows since the observed flow mixes the flow of several routes. Comparing the estimation of the dynamic route flow of routes 7 and 4, a better performance can be observed in the results of route 4 since the estimation of this one can be improved using the travel time information, i.e., the equations from the model (23)–(26) related to this flow are:

$$W_{s_3}^{a_3}\left(t\right)=E_{r_3}^{a_3}\left(t\right)+E_{r_4}^{a_3}\left(t\right)$$

(32)

$$E_{r_3}^{a_3}\left(t_k\right)-E_{r_4}^{a_3}\left(t_{k-1}\right)={\phi }_{s_3}^{s{r}_3}(W_{s_3}^{a_3}\left(t_k\right)-W_{s_3}^{a_3}\left(t_{k-1}\right))$$

(33)

However, the estimation of route 7 cannot be improved using the travel time information since the scanned vehicles of routes 1 and 7 share the same sub-route (see Table 4), and hence the prior information plays a very important role in the flow estimation. Finally, similar comments can be made when analyzing the results of traffic dynamic estimation of route 8: in this case, the observed field data mixes the flow of routes 5, 6, 8 and 9 (see again Figure 3) and the travel time data improve the results since routes 5 and 9 on one side and 6 and 9 on the other share the same sub-route (see Table 4).

Regarding the convergence of the algorithm, the use of the relaxation factor avoids oscillations in the results between iterations so that the final estimation can be achieved in a few iterations. Note also that the design of the algorithm allows to decrease the CPU (central processing unit) time to complete each iteration as the approach progresses. This is because we start from a free flow solution for adjusting the splines (Step 0), which is far from the final travel time achieved in Step 1, while in the final iterations, the used splines are very close to the results. The same happens with ${\widehat{E}}_r^a\left(t_k\right)$ of Step 3.1. This leads to a CPU time of 20 s for the first iteration and 10 s for the last one.

6.2. The Nguyen-Dupuis Network

In this section, the well-known Nguyen-Dupuis network is used to illustrate the proposed method. The network of Figure 8 shows its 13 nodes and 38 links shadowing the origin-destination nodes and highlighting the links equipped with ANPR sensors ($\mathit{SL}=\{1, 2, 3, 5, 8, 9, 16, 20, 31, 34\}$). We have considered 50 routes, which are shown in

Table 6. Set of routes for the Nguyen-Dupuis network example.

Route	Links						Route	Links
1	1	11	14	18	20		26	6	38	24	9
2	2	35	14	18	20		27	5	33	28	23
3	2	36	20				28	6	38	23
4	1	11	14	19	31		29	5	32	17	16
5	1	11	15	29	31		30	5	33	27	16
6	1	12	25	29	31		31	7	11	14	18	20
7	1	12	26	37			32	8	25	29	30
8	2	35	14	19	31		33	8	25	29	32	18	20
9	2	35	15	29	31		34	8	25	29	31
10	1	11	14	18			35	8	26	37
11	2	35	14	18			36	7	11	14	18
12	2	36					37	8	25	29	32	18
13	3	21	17	13	9		38	21	17	13	9
14	3	21	17	16	34		39	21	17	16	34
15	3	22	34				40	22	34
16	3	21	17	13	10		41	21	17	13	10
17	3	21	19	33	28	23	42	21	19	33	28	23
18	4	33	28	23			43	21	17	16
19	3	21	17	16			44	22
20	3	22					45	35	14	18	20
21	5	32	17	13	9		46	36	20
22	5	32	17	16	34		47	35	14	19	31
23	5	33	27	13	9		48	35	15	29	31
24	5	33	27	16	34		49	35	14	18
25	5	33	28	24	9		50	36

. The parameters used for the density-cost function were: ${\beta }_a=1, {\delta }_a=1/3$ and ${\gamma }_a=3$. Assuming that the coordinates of the nodes of Figure 8 are in km, and to define $D_a^0$ and $x_a^{max}$ parameters, a free flow speed of 100 km/h and a maximum density of 80 vehicles/km respectively, were used. Finally, to simulate both the data collected by the sensors and the prior route flow information, linear combinations of two normal densities were used in a similar way as in Section 6.1.

Table 7. Scanning map and its corresponding equilibrium flow equations.

	$\mathbf{Links}\mathbf{in}\mathit{S}\mathit{L}$										Equilibrium Flow (2)
$\mathit{S}\mathit{C}\mathit{S}\mathit{L}$	${\mathit{a}}_{\mathbf{1}}$	${\mathit{a}}_{\mathbf{2}}$	${\mathit{a}}_{\mathbf{3}}$	${\mathit{a}}_{\mathbf{5}}$	${\mathit{a}}_{\mathbf{8}}$	${\mathit{a}}_{\mathbf{9}}$	${\mathit{a}}_{\mathbf{16}}$	${\mathit{a}}_{\mathbf{20}}$	${\mathit{a}}_{\mathbf{31}}$	${\mathit{a}}_{\mathbf{34}}$	${\mathit{W}}_{\mathit{s}}^{\mathit{a}}\mathbf{\left(}\mathit{t}\mathbf{\right)}={\displaystyle \mathbf{\sum }_{\mathit{r}}}{\mathit{\lambda }}_{\mathit{r}}^{\mathit{s}}{\mathit{E}}_{\mathit{r}}^{\mathit{a}}\mathbf{\left(}\mathit{t}\mathbf{\right)}$
$s_1$	X							X			$W_1^1\left(t\right)=E_1^1\left(t\right)$
$s_2$		X						X			$W_2^2\left(t\right)=E_2^2\left(t\right)+E_3^2\left(t\right)$
$s_3$	X								X		$W_3^1\left(t\right)=E_4^1\left(t\right)+E_5^1\left(t\right)+E_6^1\left(t\right)$
$s_4$	X										$W_4^1\left(t\right)=E_7^1\left(t\right)+E_{10}^1\left(t\right)$
$s_5$		X							X		$W_5^2\left(t\right)=E_8^2\left(t\right)+E_9^2\left(t\right)$
$s_6$		X									$W_6^2\left(t\right)=E_{11}^2\left(t\right)+E_{12}^2\left(t\right)$
$s_7$			X			X					$W_7^3\left(t\right)=E_{13}^3\left(t\right)$
$s_8$			X				X			X	$W_8^3\left(t\right)=E_{14}^3\left(t\right)$
$s_9$			X							X	$W_9^3\left(t\right)=E_{15}^3\left(t\right)$
$s_{10}$			X								$W_{10}^3\left(t\right)=E_{16}^3\left(t\right)+E_{17}^3\left(t\right)+E_{20}^3\left(t\right)$
$s_{11}$			X				X				$W_{11}^3\left(t\right)=E_{19}^3\left(t\right)$
$s_{12}$				X		X					$W_{12}^5\left(t\right)=E_{21}^5\left(t\right)+E_{23}^5\left(t\right)+E_{25}^5\left(t\right)$
$s_{13}$				X			X			X	$W_{13}^5\left(t\right)=E_{22}^5\left(t\right)+E_{24}^5\left(t\right)$
$s_{14}$						X					$W_{14}^9\left(t\right)=E_{26}^9\left(t\right)+E_{38}^9\left(t\right)$
$s_{15}$				X							$W_{15}^5\left(t\right)=E_{27}^5\left(t\right)$
$s_{16}$				X			X				$W_{16}^5\left(t\right)=E_{29}^5\left(t\right)+E_{30}^5\left(t\right)$
$s_{17}$								X			$W_{17}^{20}\left(t\right)=E_{31}^{20}\left(t\right)+E_{45}^{20}\left(t\right)+E_{46}^{20}\left(t\right)$
$s_{18}$					X						$W_{18}^8\left(t\right)=E_{32}^8\left(t\right)+E_{35}^8\left(t\right)+E_{37}^8\left(t\right)$
$s_{19}$					X			X			$W_{19}^8\left(t\right)=E_{33}^8\left(t\right)$
$s_{20}$					X				X		$W_{20}^8\left(t\right)=E_{34}^8\left(t\right)$
$s_{21}$							X			X	$W_{21}^{16}\left(t\right)=E_{39}^{16}\left(t\right)$
$s_{22}$										X	$W_{22}^{34}\left(t\right)=E_{40}^{34}\left(t\right)$
$s_{23}$							X				$W_{23}^{16}\left(t\right)=E_{43}^{16}\left(t\right)$
$s_{24}$									X		$W_{24}^{31}\left(t\right)=E_{47}^{31}\left(t\right)+E_{48}^{31}\left(t\right)$

shows the scanning map with the feasible combinations of scanned links, which allows to obtain the observed traffic flow and to derive the Equation (2).

Table 8. Sub-route links and its corresponding travel time Equation (25).

sr	Sub-Route Links					${\displaystyle \mathbf{\sum }_{\mathit{r}}}{\mathit{\sigma }}_{\mathit{r}\mathbf{,}\mathit{a}}^{\mathit{s}\mathit{r}}\left({\mathit{E}}_{\mathit{r}}^{\mathit{a}}\mathbf{\left(}{\mathit{t}}_{\mathit{k}}\mathbf{\right)}\mathbf{-}{\mathit{E}}_{\mathit{r}}^{\mathit{a}}\mathbf{\left(}{\mathit{t}}_{\mathit{k}\mathbf{-}\mathbf{1}}\mathbf{\right)}\right)\mathbf{=}{\displaystyle \mathbf{\sum }_{\mathit{r}}}{\mathit{\phi }}_{\mathit{s}}^{\mathit{s}\mathit{r}}\mathbf{\left(}{\mathit{t}}_{\mathit{k}}\mathbf{\right)}\left({\mathit{W}}_{\mathit{s}}^{\mathit{a}}\mathbf{\left(}{\mathit{t}}_{\mathit{k}}\mathbf{\right)}\mathbf{-}{\mathit{W}}_{\mathit{s}}^{\mathit{a}}\mathbf{\left(}{\mathit{t}}_{\mathit{k}\mathbf{-}\mathbf{1}}\mathbf{\right)}\right)$
$s{r}_1$	2	35	14	18	20	$\left(E_2^2\left(tk\right)-E_2^2\left(tk-1\right)\right)={\phi }_2^1\left(tk\right)\left(W_2^2\left(tk\right)-W_2^2\left(tk-1\right)\right)$
$s{r}_2$	2	36	20			$\left(E_3^2\left(tk\right)-E_3^2\left(tk-1\right)\right)={\phi }_2^2\left(tk\right)\left(W_2^2\left(tk\right)-W_2^2\left(tk-1\right)\right)$
$s{r}_3$	1	11	14	19	31	$\left(E_4^1\left(tk\right)-E_4^1\left(tk-1\right)\right)={\phi }_3^3 \left(tk\right)\left(W_3^1\left(tk\right)-W_3^1\left(tk-1\right)\right)$
$s{r}_4$	1	11	15	29	31	$\left(E_5^1\left(tk\right)-E_5^1\left(tk-1\right)\right)={\phi }_3^4\left(tk\right)\left(W_3^1\left(tk\right)-W_3^1\left(tk-1\right)\right)$
$r{s}_5$	1	12	25	29	31	$\left(E_6^1\left(tk\right)-E_6^1\left(tk-1\right)\right)={\phi }_3^5\left(tk\right)\left(W_3^1\left(tk\right)-W_3^1\left(tk-1\right)\right)$
$s{r}_6$	2	35	14	19	31	$\left(E_8^2\left(tk\right)-E_8^2\left(tk-1\right)\right)={\phi }_5^6\left(tk\right)\left(W_5^2\left(tk\right)-W_5^2\left(tk-1\right)\right)$
$s{r}_7$	2	35	15	29	31	$\left(E_9^2\left(tk\right)-E_9^2\left(tk-1\right)\right)={\phi }_5^7\left(tk\right)\left(W_5^2\left(tk\right)-W_5^2\left(tk-1\right)\right)$
$s{r}_8$	5	32	17	13	9	$\left(E_{21}^5\left(tk\right)-E_{21}^5\left(tk-1\right)\right)={\phi }_{12}^8\left(tk\right)\left(W_{12}^5\left(tk\right)-W_{12}^5\left(tk-1\right)\right)$
$s{r}_9$	5	33	27	13	9	$\left(E_{23}^5\left(tk\right)-E_{23}^5\left(tk-1\right)\right)={\phi }_{12}^9\left(tk\right)\left(W_{12}^5\left(tk\right)-W_{12}^5\left(tk-1\right)\right)$
$s{r}_{10}$	5	33	28	24	9	$\left(E_{25}^5\left(tk\right)-E_{25}^5\left(tk-1\right)\right)={\phi }_{12}^{10}\left(tk\right)\left(W_{12}^5\left(tk\right)-W_{12}^5\left(tk-1\right)\right)$
$s{r}_{11}$	5	32	17	16	34	$\left(E_{22}^5 \left(tk\right)-E_{22}^5\left(tk-1\right)\right)={\phi }_{13}^{11}\left(tk\right)\left(W_{13}^5 \left(tk\right)-W_{13}^5 \left(tk-1\right)\right)$
$s{r}_{12}$	5	33	27	16	34	$\left(E_{24}^5\left(tk\right)-E_{24}^5 \left(tk-1\right)\right)={\phi }_{13}^{12}\left(tk\right)\left(W_{13}^5 \left(tk\right)-W_{13}^5\left(tk-1\right)\right)$
$s{r}_{13}$	5	32	17	16		$\left(E_{29}^5\left(tk\right)-E_{29}^5\left(tk-1\right)\right)={\phi }_{16}^{13}\left(tk\right)\left(W_{16}^5\left(tk\right)-W_{16}^5\left(tk-1\right)\right)$
$s{r}_{14}$	5	33	27	16		$\left(E_{30}^5\left(tk\right)-E_{30}^5\left(tk-1\right)\right)={\phi }_{16}^{14}\left(tk\right)\left(W_{16}^5\left(tk\right)-W_{16}^5\left(tk-1\right)\right)$

shows the resulting sub-routes and the corresponding Equation (25), which allows to improve the information for estimating dynamic route flows, taking advantage of the measured travel times.

Figure 9 shows some results after the application of the proposed model on the Nguyen-Dupuis network. As expected from the analysis of Table 7, the flow estimation of route 1 (fist column of Figure 9) gives almost the assumed real flow since it can be derived directly from the observed flow $W_1^1\left(t\right).$ Note that despite having assumed a prior flow (dashed lines in the figure) very different from the real flow (continuous-filled curve), the model is able to reproduce the real shape (dashed-dot line) of the curve, even the important distortion due to the congestion produced by link 14 where the wave crests almost disappear and the flow intensity is quite constant from 11:00 to 18:00 h. Similar curves were derived for routes 2, 13, 14, 15, 19, 27, 33, 34, 39, 40 and 43 (see Table 7).

The influence of including the travel time information of the observed vehicles in the model that come from different routes, but they are included in the same observed flow $W_s^a\left(t\right)$, is patent, for example, in the estimation of the flow of route 6 (see the second column of Figure 9). Despite having observed users of route 6 mixed with users of routes 4 and 5 in $W_3^1\left(t\right)$, Equation (25), in this case, written on $s{r}_3, s{r}_4$ and $s{r}_5$ (see Table 8), allows the model to obtain very good route flow estimation. For this particular case, 14 sub-routes were obtained which allows to improve the dynamic flow estimation of 14 more routes (2, 3, 4, 5, 6, 8, 9, 21, 22, 23, 24, 24, 29 and 30).

When the observed $W_s^a\left(t\right)$ split has to be done based only on the prior information, the final estimation, although improved, is not as good as desired. For example, the split of $W_{24}^{31}\left(t\right)$ to estimate the dynamic flow of routes 47 and 48 has to be done based on a prior flow that was far from the real one. We have chosen the link 35 to illustrate the evolution of the dynamic route estimation $e_{35}^{47}\left(t\right)$ and $e_{35}^{48}\left(t\right)$ in the upper right corner of Figure 9, where the influence of prior flows in the final estimation is patent due to the lack of more information.

Finally, the estimation of the link flow intensity also depends on the quality of the route flow estimation. For example, for the dynamic link flow of link 1 (see the bottom right corner of Figure 9), the model perfectly adjusts the shape of the total link flow (upper purple curve) because a sensor has been installed on it. However, some estimations of the route flows that use the link do not adjust to the real one, because the route estimation of routes 4, 5, 6, 7 and 10 are not exactly the real ones. The estimation of dynamic flow for link 14, despite being a very good estimation, does not fit with the real one due to the complexity of involving more routes with no perfect estimation.

All the comments above show the great influence of the correct ANPR sensor location to achieve a better quality of dynamic route flow estimation, and that deserves deeper research in a future work.

The number of iterations needed to reach the final estimation was 10. Indeed, the different plots of Figure 9 show the results for iteration 3, which in general is already very close to the final one.

6.3. The Eastern Massachusetts Network

The last example to show the performance of the proposed model consists of a large-scale road network. Specifically, the network modelled for this example is the Eastern Massachusetts road network shown in Figure 10, used, for example, by Zhang et al. [61] or Liang et al. [62]. The network is made up of 258 links and 74 nodes. Both for the definition of the network topology, as well as the values related to capacities, free-flow costs and parameters associated to the density-cost functions of the links, the data available in an open database have been used (Stabler [63]). The traffic demand for this network has been modelled assuming that there are 731 O-D pairs in the network, and whose trips are distributed over a total of 1452 routes. We have used the model proposed in Mínguez et al. [38] to locate 100 ANPR sensors, which are also indicated in Figure 10 (note that a total of 168 sensors are needed to be able to differentiate the users of all routes).

The aim of this section is to test the proposed model in a real network and to check that the computation times are reasonable. Regarding this last issue and taking into account that we have used a non-optimized code and some links results with a saturation above 0.8, the number of iterations required was 15 and the total CPU time was 8 h. The flow estimation results’ trends corresponding to this network are very similar to those of previous sections, and some results are shown in Figure 11, where again, the colors refer to the different route components but due to the great number of routes involved, the legend in this case has not been included.

For example, since the users of route 63 can be identified, the model is able to reproduce the real shape of this intensity flow. Route 41 is one of the routes for which estimation is improved due to the travel time information, and therefore, although the real shape is not completely reproduced, the final results are fairly good. Again, when only prior information is available, the final dynamic route flow estimation has a clear influence of these flows, as can be seen, for example, in the estimation of routes 47 and 48 (see the upper right corner of Figure 11). Finally, the figure shows the link flow intensity estimation of link 119 with the corresponding path flow components. We have chosen this link because it is the link with higher congestion of the network.

7. Conclusions and Future Research

7.1. Conclusions

The main conclusions that can be drawn from this paper are as follows:

• The proposed model is based on the interesting results of Castillo et al. [22] (continuous dynamic estimation of travel times, traffic flows in links and routes, etc.). One of the steps forward in this paper is the use of field data from ANPR sensors to estimate the dynamic route flows. In addition, the examples of applications exposed in this paper have proven the good performance of the estimations.
• Because the link exit–entry time functions ${\tau }_a^{-1}\left(t\right)$ were approximated by monotone cubic-spline functions, the model, which satisfies the FIFO rule, is able to model the link travel time at any time, which allows to easily introduce the field data (i.e., data of each single scanned vehicle) in order to reconstruct the dynamic route flow and to evaluate its travel time. In any case, if the approximation is not enough, the number of $t_k$ basic points can be increased.
• The knowledge of the dynamic traffic intensities during the day is very useful for authorities and transportation engineers and planners to know the mobility guidelines (i.e., departure times, routes, speeds, etc.). For example, for a peak hour of a given link, the proposed model allows (using ${\tau }_a^{-1}\left(t\right))$ to identify the worst time to start each route. This information can be used to develop sustainable policies which, for example, persuade certain users to modify their mobility patterns in order to improve the global performance of the network (switch to a sustainable mode, create a new bus route, etc.).
• The model developed in this paper is not proposed to provide results in real-time. On the contrary, the aim of the proposed model is to estimate the flows of a network using some observed variables which do not need to be obtained with a low CPU time. In any case, the used code was not optimized, and a more careful programming may lead to a faster development.
• The model has been applied to medium-sized networks, which made it useful in the engineering practice. For example, regarding sustainable transportation applications, the proposed approach can be integrated in a model to quantify the vehicle emissions throughout the day, which is a topic of great interest in many recent works (see Wang et al. [5] for more detail). The results can also be used to provide the basic historical day-by-day data, which can be used to predict (in real time) traffic conditions in the near future (see Castillo et al. [64] as an example).

7.2. Future Research

This paper represents a starting point on a topic that deserves deeper research. In particular, we think that some aspects to further develop are as follows:

• The existing models to locate ANPR sensors (such as the one used to locate the sensors in Figure 10) do not consider the advantages of considering restriction (25). The results shown in this paper prove the advantages of taking into account the travel time information to improve the dynamic traffic flow estimation. Therefore, a sensor location model to take even more advantage of the data from ANPR is needed.
• The results of the GLS-based optimization model proposed in this paper (Step 3 of the algorithm) is greatly influenced by the prior information, which may be very different to the real flow performance. Therefore, it is necessary to explore other models capable of reconstructing the cumulative route flow.
• The examples of applications shown in this paper are based on simulated traffic flows, including the vehicles’ plates and passing times. However, it is true that since we have used the proposed Continuous Dynamic Network Loading to simulate those plates, a real experiment needs to demonstrate the real power of the model. In particular, the real data should be used to calibrate the parameters of the used density-flow function and to evaluate the errors in the field data collected by the ANPR sensors. In any case, since the real intensity flows usually look like those simulated in the examples, we expect that the model would reproduce the real situation fairly well.
• There is a proven interest in quantifying the emissions of the vehicles in a network. For that reason, it would be interesting to combine the proposed model with an emission model, such that the estimated flows can be introduced as input to quantify the vehicle emissions throughout a day.