Models, Algorithms and Applications of DynasTIM Real-Time Traffic Simulation System

Abstract

Intelligent Transportation Systems (ITS) have the potential to improve traffic conditions and reduce travel delays. As a decision support software system for ITS, DynasTIM is based on the principle of dynamic traffic assignment and developed for real-time online simulation, prediction and optimization of dynamic traffic flows in urban or expressway networks. This paper introduces the models, algorithms and some typical applications of DynasTIM. The main contents include: the functional architecture; the application architecture of the system; dynamic OD (Origin-Destination) flows estimation method with novel formula for assignment matrix computation; mesoscopic traffic model using variable-length speed influence region and calibrating speed online based on connected vehicles data; and parallel SPSA algorithm based urban area signal optimization method. The functions of DynasTIM are implemented basically through three main modules: state estimation (ES), state prediction and control strategy optimization (PS&CSO), and guidance strategy optimization (GSO). The case study is aimed at the populated Futian Central Business District (CBD) road network in Shenzhen, China, which has an area of about 7 square kilometers. Based on the archived turning counts collected from 359 video traffic detection locations, DynasTIM was calibrated offline for this network, in order to validate the capability of simulating actual traffic conditions, and to set up basic conditions for testing signal optimization methods. The results show that the simulation output flows of DynasTIM have fairly good matching accuracy with the real surveillance flows in the field. Furthermore, for the CBD network with 38 signalized intersections, the signal optimization method is evaluated and better signal timing plans are found which can reduce about 13% average travel delay, compared with the signal plans currently implemented in the field.

1. Introduction

Intelligent Transportation Systems (ITS) have the potential to reduce traffic congestion in urban and intercity traffic networks. For ITS to be effective, it is necessary to optimize dynamic network traffic management scenarios, probably including coordinated signal controls, traffic guidance, travel incentives, dynamic congestion pricing and optimal combinations of these strategies. Additionally, finding optimal scenarios has strict time constraint for real-time response to the current and unfolding dynamic network conditions. However, determining the optimal solution directly is rarely feasible, and a more realistic way is to keep trying, evaluating and modifying scenarios until the optimal or satisfactory one is found.

Since actual traffic network systems are usually quite complex, real-time online traffic simulation is presently almost the only practicable and risk-free online evaluation means. For each iteration in the optimization process, the effectiveness of the current scenario can be evaluated by an online traffic simulator, and the results are returned to the optimization algorithm to direct the adjustment in the next iteration. The process is often iterated hundreds of times so that a satisfactory optimization result can be achieved, and should be completed within several minutes, generally. DynasTIM (Dynamic network assignment system for Traffic Information Management) is such a software system that can be applied to the real-time simulation, prediction and optimization of dynamic network traffic flows, therefore providing information and decision support for ITS. It was originally developed as a simulation based Dynamic Traffic Assignment (DTA) system, followed by adding some novel features including system architectures, area network signal optimization based on the parallel SPSA algorithm, and the mesoscopic traffic model with variable-length speed influence region and online speed calibration functionality via connected vehicles data.

The introduction to the system structure and case studies of DynasTIM was initially published in a conference paper by the author [1]. This paper is an expanded version of the conference paper, mainly adding the introduction to key models and algorithms of the system. The organization of this paper is as follows. In Section 2, related work is discussed. Section 3 describes the architecture of DynasTIM including state estimation, state prediction and strategy optimization modules. Then, the application architecture of DynasTIM is introduced in Section 4. What follows is a description of dynamic demand analysis in Section 5, and supply simulation in Section 6. Section 7 introduces the signal optimization method. The software design is presented in Section 8. Then, case studies are given in the city of Shenzhen, China, where DynasTIM was calibrated and the signal optimization algorithm was evaluated. Lastly, conclusions and some suggestions for future work are given.

2. Literature Review

For dynamic traffic assignment, numerous models and solution methods have been established since the advent of Intelligent Transportation Systems (ITS), which range from mathematical programming, optimal control, variational inequality to simulation-based models. More details can be found from the reviews by Peeta and Ziliaskopoulos [2], and Chiu et al. [3]. This paper focuses on the simulation-based DTA systems. Although they lack the theoretical insights compared to analytical DTA models, the main advantages of which, over an analytical one, are those that the traffic dynamics is usually more realistic, e.g., the FIFO rules can be violated with the consideration of the speed difference across lanes due to signal controls or the spillback in actual networks, and the modeling of personalized traffic guidance strategy and complex signal control logic is more convenient as well.

In general, simulation-based DTA systems include two key components: one used to calculate the time-varying path flows with a specified service level; and the other to simulate vehicle movements along paths and obtain time-dependent traffic conditions and measure of effectiveness (MOE) with specified path flows [4,5]. The former is called the “demand” side and is basically used to explain where vehicles in the network depart, where to go, when to leave, and which paths to choose. The latter is called the “supply” side and is used to analyze the capacities of network elements, which are time dependent and may vary with dynamic traffic signal controls, incidents, weather, etc.

The demand side estimates and predicts dynamic Origin-Destination (OD) flows in the network [6,7,8,9,10,11] and utilizes some behavioral rules or models, e.g., Wardrop principles [12], discrete choice models [13], to discretize the OD flows into a list of vehicles, each of which is assigned the corresponding origin, destination, travel path and departure time. Generally, the OD flows have to be estimated with some algorithms based on the historical information, and even the real-time surveillance data for online applications. The supply side loads the list of vehicles onto a mesoscopic or microscopic traffic simulator to simulate the complete movement trajectory of each vehicle from its origin to the destination along the given path, according to some traffic flow models, which basically differentiate various traffic simulation systems [14,15]. The trajectories of all vehicles in the network are merged, diverged, stacked or intertwined, some like ripples in a lake. As a result, the network traffic conditions are generated due to the cumulative effect of the trajectories (ripples), including flows, densities, speeds, queue lengths, travel times and the measures of network performance, etc. In particular, the assignment matrices describing the mapping relationship between dynamic OD flows and sensor flows are obtained by counting and aggregating the relevant simulation output data on the supply side, which are often the most important parameters for some dynamic OD flow estimation algorithms [7,8]. In short, a simulation-based DTA system reproduces a variety of real-life traffic dynamics characteristics, estimates and predicts network conditions, and evaluates network performance, through modeling the demand side, the supply side, and the complex interaction process between them. In the end, the resulting simulation outputs may be applied in various planning and real-time application situations.

For planning applications of a DTA system, which are generally performed offline using historical archived data to calibrate the system [16,17,18,19,20], it may be used to evaluate the impacts on network performance, of new transportation infrastructure and policies, incident management strategies, and evacuation and rescue plans for emergencies. In the situation of real–time online applications of a DTA system, which is implemented online and calibrated using real-time surveillance data [21,22], the real-time DTA system is utilized to estimate and predict network traffic states, and perform the online evaluation and optimization of signal control and travel guidance strategies, in order to provide real-time information and decision support for the operations of Traffic Management Centers (TMCs).

Originally, traffic simulation systems are mainly applied to the offline evaluation and optimization of operation scenarios. In recent years, with some more refined transportation planning projects carried out in some cities, the within-day, not just day-by-day or year-by-year, traffic dynamics must be considered to analyze the propagation characteristics of traffic congestion in a network at different times throughout the day, and great progress has also been made in computer hardware and software technologies. Therefore, simulation-based DTA systems are gradually applied to some projects to assist in the development of congestion mitigation strategies. So far, the simulation-based DTA systems that have been applied to actual networks include but are not limited to, Transmodeler [15], INTEGRATION [23], METROPOLIS [24], CONTRAM [25], DYNASMART [26,27,28,29,30], DynaMIT [31,32,33,34,35,36,37,38], VISTA [39], Dynameq [40,41,42,43], AIMSUN [44,45,46], DynusT [47,48] and DTALite [49]. With the rapid development of high performance computation techniques, traffic data collection means and communication technologies, it is highly possible that some of the simulation-based DTA systems will be applied in the real-time online situations and for large-scale traffic networks. In this paper, these systems are also called real-time online traffic simulation systems as well.

To this day, most representative real-time online traffic simulation systems may be DynaMIT and DynaSmart, for which the key models and algorithms concerning dynamic OD matrix estimation, offline/online calibration are extensively studied, and a large number of field tests are implemented. However, it is still necessary for these systems to: improve the real-time traffic estimation and prediction accuracy; to model travel choice behavior more realistically in the background of large amounts of travelers’ preference and trajectory data available; to establish more efficient and effective OD flows estimation methods for large scale networks; and to apply simulation based evaluations to traffic control strategy optimization, etc. For the optimization of signal control plans of urban areas, Carolina et al. [50,51,52,53] propose a simulation based optimization method, which results in well-performing signal plans and is highly computationally efficient, through the integration of a derivative-free trust region algorithm and a metamodel that combines a microscopic simulator with an analytical queueing network model. For this method, there are still some questions that are not clear enough, such as: whether it can be used for the real-time optimization of urban area signal controls and how large scale; whether the queuing network model allows for the change of network conditions due to traveler’s response to informed information; and whether the queueing network model is capable of describing physical spillback queues, especially for Manhattan-like networks, since the modeling of some basic traffic dynamics may be important for the effective evaluation of signal control strategies. This paper also tries to answer some of the above mentioned questions.

3. Architecture of DynasTIM

In Figure 1, the overall structure of DynasTIM is illustrated. The main functionality modules are state estimation (ES), state prediction and control strategy optimization (PS&CSO) and guidance strategy optimization (GSO). The basic principles of state estimation (ES) and state prediction (PS) are similar to those of DynaMIT [32]. The ES and PS processes both apply the demand simulation and the supply simulation to implement the estimation and prediction of network conditions. The major functional differences between DynasTIM and other simulation based DTA systems, e.g., DynaMIT and DynaSmart-X, include the signal optimization and the guidance strategy optimization.

3.1. State Estimation and State Prediction

In the stage of the state estimation, DynasTIM, based on historical and real-time input data, makes use of the “demand simulation” module to estimate the dynamic OD flows for the current time interval (estimation interval) and, through the travel choice behavior models, disaggregate the OD flows into a list of travelers with different origins, destinations, departure times, travel paths and personal preferences. Then, the system loads the “travelers” onto the “supply simulator” to implement the fast simulation for the dynamic propagation process of network flows, so that the estimation of entire network conditions can be obtained from the simulation outputs which include flows, speeds, densities, queue lengths, link travel times, OD travel times, network performance indicators, etc. If there are considerable deviations between the estimated conditions and the corresponding real-life surveillance data, the dynamic OD flows and the model parameters are calibrated online by an iteration process, where the link impedances and the assignment matrices reflecting the mapping relationships between dynamic OD flows and surveillance flows at sensor locations are updated iteratively, through the Method of Successive Average (MSA). The iteration process does not end until the simulation outputs are consistent with the real-life surveillance data or the number of the iterations exceeds the upper limit.

For the supply simulator, also known as the mesoscopic traffic simulator, the macroscopic speed-density (k-v) relationship is used to calculate vehicle speeds for uninterrupted flow facilities, and the queue model and the capacity model to describe the queuing behavior of vehicles approaching to bottlenecks (e.g., intersections and lane drops). Compared to the microscopic traffic simulator, in which the detailed vehicle movements, e.g., car following, lane changing and overtaking, are modeled, the mesoscopic traffic simulator with higher simulation efficiency is capable of the real-time online applications and easier to create the road network and calibrate fewer model parameters. Compared to the macroscopic traffic simulator, the mesoscopic one simulates and traces the movement of each vehicle in the network, and depicts the response of travelers with different attributes (social economic attribute, values of time, trip purpose, information source, etc.), to diverse traffic information, this capability of which is especially important for the modeling of traffic guidance strategies that already have substantial effects on network conditions nowadays.

With the state estimation results as the start point, the “demand simulation” module in the stage of state prediction uses an autoregressive model to predict dynamic OD flows, which are then loaded onto the mesoscopic traffic simulator to simulate vehicle movements in order to generate short-term prediction of network conditions. Likewise, the state prediction is an iterative process and the MSA is used to update the link impedances as well. Theoretically, the iteration process produces a network assignment result closer to the SUE (Stochastic User Equilibrium) [54]; in practice, it generates the link impedances that are more consistent with actual travel experiences and applied to the generation of proactive guidance strategies, to avoid some unexpected “overreaction” (see, e.g., Ben-Akiva et al. [32] for more details) of travelers to the predicted network conditions.

In fact, state estimation is to make the traffic simulation system reproduce the current traffic conditions as accurately as possible, while state prediction is to predict the future traffic conditions based on the current situation. The two share the demand simulation and supply simulation modules. The difference is that the main inputs of state estimation are the current dynamic OD flows obtained by dynamic OD estimation algorithm based on historical data and real-time surveillance data, while the main inputs of state prediction are the predicted dynamic OD flows, which can be obtained based on the autoregressive model.

3.2. Strategy Optimization

The main purpose of state prediction is to achieve the active traffic demand and supply management. Therefore, after the completion of the state prediction, the system optimizes the signal control strategy proactively to improve the network supply capacity. When the optimized control strategy is confirmed by TMC to be implemented in the field, DynasTIM will make another online simulation evaluation based on this control strategy to generate the updated network state prediction, as the start point of the subsequent “GSO” process. Usually, the control strategy is optimized every 5~60 min and just optimized for one time in each outermost “ES→ PS&CSO→ GSO” run cycle. Additionally, the list of travelers generated in the state prediction process, which are obtained by discretizing the predicted dynamic OD matrix, will be stored as one of the input data for the GSO process.

The GSO process helps travelers optimize their departure times, trip modes and travel paths by providing them with the pre-trip and en-route travel information based on the predicted traffic conditions, thereby proactively influencing the network demand. What is different from DynaMIT is, DynasTIM, after completing the state prediction and control strategy optimization, does not wait for the next estimation interval of surveillance data available for the following run cycle, but immediately begins receiving real-time detection data, and the data aggregation interval is usually 30 to 60 s. These may include connected vehicles data for the online correction of the calculated speed of vehicles in the mesoscopic traffic simulator, incident data for the online calibration of the lane group capacities, the OD information (e.g., origin, destination, and departure time) of each reserved trip for the online calibration of dynamic travel demand, etc. Presently, the GSO process does not necessarily predict the dynamic OD flows, but instead directly uses the predicted traveler list generated in the PS&CSO process as input. The process can update the state prediction results generated in the previous PS&CSO process, similarly through the demand-and-supply simulation, particularly taking the impact of incidents into consideration to get more timely response to capacity changes in the network.

In the future, with the large-scale deployment of intelligent and connected vehicles, it is possible that all vehicle trips can be scheduled by TMCs, or the destination of each vehicle in the network will be known by the data analysis. As a result, the impact of each trip on network conditions, perhaps referred to as “ripple effect”, can be predicted by the GSO process. Moreover, through the optimization of each network user’s travel choice, this process is beneficial to the improvement of overall network performance caused by the cumulative ripple effects. Generally, the guidance strategy is updated every 30~60 s and can be generated according to the Stochastic User Equilibrium (SUE) principle [54]. The GSO process has been running till the surveillance data for the next state estimation interval are available. Then, DynasTIM proceeds to the next “ES→ PS&CSO→ GSO” run cycle.

3.3. Rolling Horizon Implementation

Similar to DynaMIT [32], DynasTIM also uses a sliding window framework to enhance the qualities of the network state prediction and predictive-based strategy optimization, and the basic concept is depicted in Figure 2. At T = 8:00, execution cycle 1 begins. The system estimates the current network conditions at 8:00 through the integration of historical information and real-time surveillance data collected during the most recent 5 min. Based on the estimation result, a prediction process is then performed to generate the predicted network conditions as the basis for the following signal optimization. After completing the optimization process, DynasTIM proceeds to accept real-time surveillance data and optimize guidance strategies till the time is 8:05, when the surveillance data for cycle 2 are available. At T = 8:05, a new execution cycle begins.

4. Application Architecture of DynasTIM

Figure 3 shows the connections between DynasTIM and other application modules. DynasTIM is deployed at the Traffic Management Center (TMC) and integrates historical information with real-time input (surveillance data, dynamic traffic control scenarios, etc.), to achieve the following functions:

• Real-time estimates and predictions of complete network traffic states, including flows, densities, speeds, queue lengths, travel times for all of network links, and dynamic OD flows, etc.;
• The optimization of traffic signal control and travel guidance strategies based on predicted traffic states to generate proactive traffic management measures.

The historical database stores the network topology and geometric characteristic data, signal control scenarios, sensor flows, speeds, occupancies, link travel times and offline calibrated parameters including mesoscopic traffic model parameters, travel choice behavior model parameters, dynamic OD flows, etc. The real-time surveillance data are collected from probe vehicles, video detectors, V2I devices, loop detectors and so forth. In addition, the incident data including the location, severity and expected duration are also important real-time input information. The historical and real-time data are fed to the front-end data fusion module for preprocessing, including the filtering to reduce input errors, the classification to obtain different types of data (such as speeds, densities, flows, travel time etc.), and the aggregation to get different time scales of data, in order to satisfy various requirements of the components of DynasTIM. For example, the aggregation interval of sensor flows is usually 5~15 min for dynamic OD flow estimation, while it is often 30~60 s for Floating Car Data (FCD) used to correct the calculated speed of vehicles in the mesocopic traffic simulator of the system.

For the signal control scenario optimization, generally there are two modes to choose, including the scenario selection mode like SCATS [55] and the scenario generation mode like SCOOT [56]. The former evaluates the performance for each of the existing signal control scenarios and selects the optimal one to recommend to the TMC. The latter can globally optimize the timing parameters for all signal controllers in the network.

With respect to the guidance strategy optimization, DynasTIM generates the predicted traffic conditions to be sent to the cloud platform for back-end data fusion. By analyzing the travelers’ preference data, the cloud platform can provide personalized travel information service for them, such as the recommendations for the optimal departure time, trip mode and travel path with the highest reliability or the shortest travel time, etc.

5. Demand

5.1. Offline Estimation of Dynamic OD Flows

Dynamic OD flows are the basic input data of DynasTIM. For the actual road network, it is currently difficult to obtain historical dynamic OD flows directly. Thus, the Recursive Lease Square (RLS) algorithm [57] with less historical information required is adopted to estimate dynamic OD flows based on sensor flows, link travel times and apriori OD flows. If Floating Car Data (FCD) are sufficient, link travel times can be calculated by analyzing the FCD; otherwise, they are obtained directly from the simulation outputs. The RLS algorithm usually requires at least one interval of apriori OD flows to initiate the estimation process, however, even on the assumption that all apriori OD flows are zero, the algorithm can still obtain the equivalent dynamic OD flows, which generate same network assignment results (e.g., link flows) as those of actual OD flows. The Generalized Least Squares (GLS) formulation for RLS algorithm is:

$$\underset{x_h}{\min } \left[{\left(x_h-x_h^a\right)}^T{W}_h^{-1}\left(x_h-x_h^a\right)+{\left(\partial y_h-a_h^h{x}_h\right)}^T{R}_h^{-1}\left(\partial y_h-a_h^h{x}_h\right)\right]$$

(1)

subject to

x_h ≥ 0

(2)

where x_h is the vector of OD flows to be estimated for current estimation time interval h, $x_h^a={\widehat{x}}_{h-1}$ takes the OD flow estimate for time interval (h − 1) as the apriori value for time interval h, $\partial y_h=y_h-{\displaystyle \sum }_{p=h-p^{\prime }}^{h-1}{a}_h^p{\widehat{x}}_p$, y_h is the vector of measured traffic counts at detectors during time interval h, ${\widehat{x}}_p$ is the vector of estimated OD flows for previous time interval p, $a_h^p$ is an assignment matrix that represents the contribution ratios of x_p or ${\widehat{x}}_p$ on y_h, p’ is the corresponding number of time intervals for the longest OD trip in the network, W_h is the variance-covariance matrix for $x_h^a$ and R_h is the variance-covariance matrix of the measurement error vector, which both are set to unit matrix in the absence of any prior knowledge.

For the RLS algorithm, the key point is the calculation of the assignment matrices $a_h^p$ (p = h, h − 1, …, h−p). In this paper, an analytical computation method for an assignment matrix is proposed, each element of which is calculated by the following formulation:

$$a_{lh}^{rp}={{\displaystyle \sum }}_{k\in K_r}{\phi }_{lh}^{kp}{\rho }_k^{rp}$$

(3)

where, $a_{lh}^{rp}$ at the l^th row and the r^th column of matrix $a_h^p$ represents the contribution ratio of the r^th OD flow during departure time interval p on the measured traffic count at detector l during time interval h; ${\rho }_k^{rp}$ is the choice probability of travelers, corresponding to OD pair r and departure time interval p, for path k from K_r which is the path set for this OD pair; ${\phi }_{lh}^{kp}$ represents how much percentage of the k^th path flow departing the origin during time interval p, contributes to the measured traffic count at detector l during time interval h. ${\rho }_k^{rp}$ is usually calculated by some path choice model, e.g., Path-size Logit model [13], and ${\phi }_{lh}^{kp}$ is computed by the following expression.

$${\phi }_{lh}^{kp}=ℒ\left({\mathsf{\Gamma }}_h{\displaystyle \cap }{\Delta }_l^{kp}\right)∕ℒ\left({\Delta }_l^{kp}\right)$$

(4)

where Γ_h = [ (h-1)T, hT] is the time period corresponding to interval h, T is the duration of arbitrary estimation time interval or departure interval; ${\Delta }_l^{kp}=\left[t_{1l}^{kp}, t_{2l}^{kp}\right]$ is the time window in which all the travelers for OD pair r and departure time interval p pass through sensor l, $t_{1l}^{kp}$ and $t_{2l}^{kp}$ are the arrival times at sensor l for travelers departing the origin at the beginning and at the end of departure interval p respectively and moving along path k; $ℒ$ (•) = sup (•) − inf (•), is the function to calculate the duration of the time period denoted by “•”, and sup (•) and inf (•) are separately the operators to get the least upper bound and the greatest lower bound of the time period denoted by “•”.

Proof. 

Assume for the travelers corresponding to OD pair r, who depart the origin during time interval p and move along path k, their arrivals at sensor l satisfy Poisson distributions P(λ), and there exists no overtaking on path k. Then,

$$\mathsf{\lambda }=d_k^{rp}/ℒ\left({\Delta }_l^{kp}\right),$$

where $d_k^{rp}$ is the count of travelers corresponding to OD pair r and departure interval p that choose path k, and their arrival times at the sensor l must be in the time window ${\Delta }_l^{kp}$ as defined above. As shown in Figure 4, the contribution of $d_k^{rp}$ to the measured count at detector l during time interval h just comes from the arrivals in the time window that is the intersection of ${\Delta }_l^{kp}$ and interval h. Since the arrival process satisfies Poisson distribution P(λ), the contributed traffic count is $\mathsf{\lambda }\cdot ℒ\left({\mathsf{\Gamma }}_h{\displaystyle \cap }{\Delta }_l^{kp}\right)$. Thus,

$${\phi }_{lh}^{kp}=\mathsf{\lambda }\cdot ℒ\left({\mathsf{\Gamma }}_h{\displaystyle \cap }{\Delta }_l^{kp}\right)/d_k^{rp}=ℒ\left({\mathsf{\Gamma }}_h{\displaystyle \cap }{\Delta }_l^{kp}\right)/ℒ\left({\Delta }_l^{kp}\right).$$

□

In general, $t_{1l}^{kp}$ and $t_{2l}^{kp}$ are obtained through historical link travel times which come from Floating Car Data (FCD) or corresponding simulation outputs generated by traffic simulators. As Formulations (1) and (2) is a typical Nonnegative Least Squares (NNLS) optimization problem, it can be solved by some tool software like Matlab or the corresponding algorithm (e.g., Van Benthem [58]) program coded by C/C++ languages like that of DynasTIM.

5.2. Real-Time Estimation and Prediction of Dynamic OD Flows

Although the NNLS algorithm is employed in the offline OD flow estimation problem, it lacks sufficient computation efficiency for the real-time dynamic OD flow estimation of large-scale networks. Inspired from Ashok and Ben-Akiva [7] and Bierlaire and Crittin [59], the deviations between actual and historical OD flows rather than the OD flows themselves are estimated, which eliminates the nonnegative constraint and transfers the OD estimation problem into an unconstrained Sparse Generalized Least Square (SparseGLS) optimization problem solved by highly effective LSQR algorithm [60]. The sparsity comes from the fact that each OD flow has contribution only to the traffic counts at the sensors, which are located on the corresponding OD paths and only account for a small portion of all sensors in the network. The SparseGLS formulation is given by:

$$\underset{\delta x_h}{\min } \left[{\left(\delta x_h-\delta x_h^a\right)}^T{\mathcal{W}}_h^{-1}\left(\delta x_h-\delta x_h^a\right)+{\left(\delta y_h-a_h^h\delta x_h\right)}^T{ℛ}_h^{-1}\left(\delta y_h-a_h^h\delta x_h\right)\right]$$

(5)

where

$$\delta x_h^a={{\displaystyle \sum }}_{p=h-q^{\prime }}^{h-1}{f}_{h+1}^p{\widehat{\delta x}}_p,$$

(6)

$$\mathbf{\delta }{y}_h=y_h-{{\displaystyle \sum }}_{p=h-p^{\prime }}^h{a}_h^p{x}_p^H-{{\displaystyle \sum }}_{p=h-p^{\prime }}^{h-1}{a}_h^p{\widehat{\delta x}}_p,$$

(7)

the parameter vector to be estimated is denoted by $\mathbf{\delta }{x}_h=x_h-x_h^H$, ${\widehat{\mathbf{\delta }x}}_p$ is the estimated parameter vector for previous time interval p, $x_p^H$ denotes the historical OD flows estimated by RLS algorithm for interval p, $f_{h+1}^p$ is the same auto-regressive process factor matrix on the OD flow deviations as Ashok and Ben-Akiva [7], which represents the contribution of ${\widehat{\mathbf{\delta }x}}_p$ to $\mathbf{\delta }{x}_h^a$, q′ is the number of previous time intervals that influence $\delta x_h^a$, ${\mathcal{W}}_h$ is the variance-covariance matrix for $\mathbf{\delta }{x}_h^a$ and ${\mathcal{R}}_h$ is the variance-covariance matrix of the measurement error vector which are both set to unit matrix if there are no sufficient historical data.

The computation formulation for the assignment matrices $a_h^p$ is:

$$a_{lh}^{rp}=\{\begin{array}{c}{\displaystyle \sum }_{k\in K_r}{\phi }_{lh}^{kp}{\rho }_k^{rp}, f_{rp}^H<{ℱ}^c\\ {\widehat{y}}_{lh}^{rp}/{\widehat{x}}_{rp}, f_{rp}^H\ge {ℱ}^c\end{array}$$

(8)

where ${ℱ}^c$ is the critical OD flow basically determined by computational experiments, $f_{rp}^H$ is the historical OD flow for OD pair r and departure interval p, ${\widehat{y}}_{lh}^{rp}$ is the number of travelers (vehicles) belonging to OD pair r and departure interval p that pass through sensor l during time interval h, ${\widehat{x}}_{rp}$ is the OD flow for OD pair r and departure interval p. Through the simulation based dynamic OD loading process, both ${\widehat{y}}_{lh}^{rp}$ and ${\widehat{x}}_{rp}$ are obtained from the traffic simulator to calculate the assignment matrices (AMs), which are called simulated AMs in this paper, in comparison with the analytical AMs calculated by Equations (3) and (4). However, for rather low OD flows, the simulated AMs are less reliable, because the path choice probabilities are simulated and each of the AMs is thus just one realization of the random simulation experiment, which may result in some incorrect zero elements in the simulated AMs, as no vehicles are assigned to the corresponding paths between the OD pairs with low OD flow. Therefore, for the OD pairs with low historical OD flows less than ${ℱ}^c$, the corresponding OD flows are estimated by the analytical AMs which are computed by the less fluctuating link travel times and the calculated path choice probabilities.

Although the SparseGLS algorithm is mainly used for the real-time estimation of dynamic OD flows, the algorithm is also capable of improving the offline estimation accuracy for OD flows obtained by the RLS algorithm as well. Usually, given the effective historical dynamic OD flow, the simulated AMs for the larger OD flows are more accurate compared to the analytical AMs that make the approximate assumption of vehicle arrivals at sensors satisfying the Poisson distribution, thus more accurate OD estimation results are obtained.

The prediction of dynamic OD flows is implemented using an autoregressive process that is same as the Kalidas Ashok’s approach [7]. The formulation is:

$$x_{h+1}=x_{h+1}^H+\mathbf{\delta }{x}_{h+1}=x_{h+1}^H+{{\displaystyle \sum }}_{p=h+1-q^{\prime }}^h{f}_{h+1}^p{\widehat{\mathbf{\delta }x}}_p$$

(9)

where $f_{h+1}^p$ is often calculated by linear regression models for each OD pair based on historical OD flows, and it is usually a diagonal matrix, as correlation across OD pairs are ignored in most cases.

5.3. Paths

The set of paths for all OD pairs is generated offline. Furthermore, through analyzing the FCD and License Plate Recognition (LPR) data, the vehicle travel trajectories between some OD pairs become known, then the feasible OD path set can be set up and the corresponding path choice probabilities can be calculated. For DynasTIM, these paths are referred to as “detection paths” and stored in a structured text file. For example, a line in the file is shown as follows:

13 0 23-7-5-3-1 0.7 0.4 0.2

where “13” is the ID of the origin node (O) contained in an origin zone; “0” is the ID of the destination node (D) contained in a destination zone; “23-7-5-3-1” is the ordered list of IDs for each link on the path (from O to D direction); the last three fields are the corresponding path choice probabilities for different types of vehicles, including cars, trucks, and HOVs.

For other OD pairs without sufficient FCD or LPR data to analyze the paths, the OD path set are established by performing the following steps:

• Use the k-shortest path algorithm [61] to generate the initial OD path set P;
• Add random perturbations to the impedances of all links in the network, simultaneously;
• Call the k-shortest path algorithm to generate new OD paths for the updated link impedances, and push the paths into the path set P;
• Repeat steps 2–3 multi times until the path set P is satisfactory.

The path choice probabilities can be calculated by Path-size Logit model [13]. However, the OD paths are pre-specified for buses or VIP fleets with scheduled routes.

When loading the network data, the system also reads all the OD paths into the computer memory, and the paths conform to the nesting requirement in order to save the memory. It means that if L₁-L₂-L₃-…-L_n is a feasible path from link L₁ to a destination node D (the downstream node of link L_n), then L₂-L₃-…-L_n must be a feasible path from link L₂ to D. For a specified destination node, all feasible paths from any link in the network to this node are assigned corresponding integer identifications (IDs). For example, a feasible path identified by the “CurrentPathID” can be represented as the following tuple:

〈CurrentLinkID, NextLinkID, NextPathID〉

where the start link on the path is identified by the “CurrentLinkID”, and “NextPathID” is the identifier of the path which must be followed out of the downstream node of the link labeled by the “NextLinkID”. The “NextPathID” is set to −1 if the downstream node of the link identified by “NextLinkID” is just the destination node. For a given destination node, the ordered list of links on any feasible path can be enumerated recursively until the “NextPathID” in the path tuple is equal to −1. Thus, the memory requirements of storing all OD path topology data are relatively modest, because each path can be depicted only by three integer IDs.

6. Supply

In this paper, supply means the network and the signal control system. A mesoscopic traffic simulator is designed to realize the supply simulation, which utilizes different classes of models to capture the complicated traffic dynamics characteristics in the network. These models include the network description, the anisotropic mesoscopic traffic model for uninterrupted flow facilities, the deterministic queuing model reflecting the effect of bottlenecks and the capacity model associated with designed road capacity constraints, signal controls, spillback and incidents.

For a linear road section with consistent width and number of lanes, which is in general defined as a segment (details for a segment is in Section 6.1), all vehicles on the segment can be divided into two parts: one is the “moving part” where different vehicles may move at different speeds calculated by individual speed-density function, which is detailed in Section 6.4; the other is the “physical queue part” where vehicles move at a relatively low and uniform speed calculated by the same deterministic queuing model as that of DynaMIT [32]. The queue part is formed usually due to the input capacity limit (spillback) of the downstream segment or the output capacity constraints of the current segment caused by the signal control or the incident. In some cases, the use of queue model may bring the additional computation advantage. For example, in the highly congested network, most of the vehicles are queuing up, using the simple queue model instead of the speed-density model (exponential function form) can therefore greatly reduce the computation overhead of mesoscopic simulation process to improve the real-time computation efficiency. Generally, the mesoscopic traffic simulator does not model the detailed movement process of vehicles inside intersections.

6.1. Network Representation

The basic elements of the road network include nodes, links, segments, lane groups, lanes and zones, etc., and the containment relationships between the elements are shown in Figure 5. The nodes correspond to intersections, traffic generation or attraction areas, and network boundaries. Typically, a zone may contain one or more nodes. A link represents a one-way road between two adjacent nodes, and can contain one or more segments. The definition of the segment in DynaMIT [32] is adopted by DynasTIM. Nodes and links together constitute the network topology description (directed graph model), as illustrated in Figure 6. The author has developed a Graphical Road traffic Network Editor (GRNE), shown in Figure 7, to facilitate the creation of these network elements. GRNE can directly read the road network file in OpenStreetMap format (XML file), automatically generate nodes, links, segments, and lane connectors. After that, one can further manually add new network elements (e.g., sensors, signals, zones, etc.) and set corresponding attributes of network elements.

A link includes one or more segments, and a segment is a road section with exactly the same geometry features (width, curvature, slope, etc.) and traffic attributes (number of lane groups, number of lanes, speed-density model parameters, etc.). Each segment can contain multiple lanes, and these lanes are divided into lane groups, of which the definition can refer to HCM2000 [62]. Each lane group has an output capacity limit at its downstream end. Most of the segments are defined in advance, however, some segments can be dynamically generated or removed during the runtime of the system. In this paper, it is assumed that an incident is always positioned at the most downstream end of the segment. When an incident occurs in the simulation process, the segment can be divided into two at the incident location and the output capacity of the upstream segment is updated to reflect the incident impact; when the incident disappears, the two segments are merged automatically to save computation costs.

6.2. Capacities

The capacity model is used to calculate the maximum number of vehicles that can leave the current lane group and enter the downstream segment (lane group) during a given time interval. Vehicles that exceed the capacity limit will remain in the current lane group to form a physical queue.

$$N_p^t=\min \left\{c_p^t, n_{p+1}{l}_{p+1}\left(k_{jam}^{p+1}-k_{p+1}^{t-1}\right)\right\},$$

(10)

where

p:

Subscript denoting a lane group. The index p increases along with vehicle traveling direction on the same link;

$N_p^t$:

Maximum number of vehicles leaving lane group p and enter lane group (p + 1) during time interval t;

$c_p^t$:

Output capacity of lane group p during time interval t;

n_p+1:

Lane number for the parent segment of lane group (p + 1);

l_p+1:

Length of the parent segment of lane group (p + 1);

$k_{jam}^{p+1}$:

Jam density for the parent segment of lane group (p + 1);

$k_{p+1}^{t-1}$:

Density for the parent segment of lane group (p + 1) during time interval (t − 1).

As a note, the parent segment of a lane group means the segment contains the lane group. For example, a segment can accommodate two lane groups, the first contains all left-turn lanes, and the other contains all of the through lanes and the right turn lanes.

6.3. Queuing part

For the queuing part, the deterministic queuing model of DynaMIT [32] is adopted in this paper. The travel delay of i^th vehicle is calculated by the formula in Equation (11):

i/c

(11)

where c is the output capacity of the corresponding lane group in the segment, that is the maximum number of vehicles leaving the queue per unit time (e.g., one second). If a vehicle arrives at the queue tail at time t, the position of it can be calculated as follows:

q(t) = q₀ + l (ct − m)

(12)

where q(t) is the position of the queue tail at time t; q₀ is the position at time 0; l = 1/k_jam is the average passenger car length, and k_jam is the jam density; ct is the total number of vehicles leaving the queue during a period of time t. As shown in Figure 8, m is the number of vehicles moving between the considered vehicle and the tail vehicle of the queue at time t = 0, that is, the number of vehicles reaching the queue before the considered one.

The applicable condition of the queue model is 0 < q(t) < L, where L is the segment length. Obviously, q(t) < 0 is not possible because the spillbacks to the upstream segments are modeled explicitly; q(t) ≥ L shows that the considered vehicle does not catch the queue, that is, the queue has disappeared when the vehicle arrives at the downstream end of the segment.

6.4. Moving Part

For each vehicle in the moving part of a segment, the basic principle of speed calculation is similar to that of the anisotropic mesoscopic model developed by Chiu et al. [14], that is, the prevailing speed of the vehicle is determined by using a macroscopic speed-density relationship based on the density in its immediate downstream speed influence region (SIR). However, in this paper, the variable-length SIR (vSIR) rather the fixed-length one is utilized to calculate the vehicle speed (as shown in Figure 9), which differs from the model of Chiu et al. Actually, when the vehicle speed is higher, a driver needs to observe the road conditions in the longer distance ahead, to adjust the current speed of the vehicle, and the SIR length is also related to the response time of the driver to the front traffic. Thus, the speed of each vehicle in the moving part of a segment can be calculated as follows.

$$k_i^{t-1}=\mathit{min}\left\{k_{jam}, \frac{N_i^{t-1}}{nl\left(v_i^{t-1}\right)}\right\}$$

(13)

$$l\left(v_i^{t-1}\right)=\mu v_i^{t-1}\times resp\_time$$

(14)

$$v_i^t=\mathit{min}\left\{v_{lim}^t, v_f{\left[1-{\left(\frac{max\left(0, k_i^{t-1}-k_0\right)}{k_{jam}-k_0}\right)}^{\beta }\right]}^{\alpha }\right\}$$

(15)

where

i:

Subscript denoting a vehicle;

t:

Superscript denoting a time interval;

$k_i^{t-1}$:

vSIR density for vehicle i during time interval (t − 1);

k_jam:

jam density, parameter to be calibrated;

k₀:

maximum density that sustains free-flow speed, parameter to be calibrated;

$N_i^{t-1}$:

Number of vehicles present in vSIR for vehicle i during time interval (t − 1), excluding vehicle I;

n:

Number of lanes present in vSIR;

$v_i^t,v_i^{t-1}$:

Speeds of vehicle i during time intervals t and (t – 1), respectively;

$l\left(v_i^{t-1}\right)$:

Length of vSIR for vehicle i during time interval (t – 1);

resp_time:

Average response time of a driver to the traffic conditions ahead, generally taken as 2~3 s or to be calibrated;

μ:

Parameter to be calibrated;

v_f:

Free-flow speed, parameter to be calibrated;

α, β:

Parameters to be calibrated;

$v_{lim}^t$:

Segment speed limit for the specified vehicle type during time interval t.

Compared with the fixed length SIR model, the AMS-vSIR model with variable-length SIR can naturally deal with the congestion effect. When a road segment is congested, the vehicle speed is low, so the length of the vSIR naturally shortens by Equation (14), vice versa. Therefore, unlike AMS with fixed SIR length, it is not necessary for the AMS-vSIR to calibrate the SIR length separately for different types of segments (such as freeways, ramps, arterials, etc.) or for different traffic conditions.

In some cities in the world, such as Beijing and Shenzhen, large-scale Floating Car Data (FCD) collection systems have been established. The upload cycle of the position and speed data for each probe vehicle can be pretty short and even less than 10 s. These probe vehicles and the intelligent and connected vehicles to be deployed on a large scale in the future can provide rich data source for the online calibration of the mesoscopic traffic simulator. In DynasTIM, a segment may be divided into multiple cells with equal length. During time interval t, if the reliable average speed estimate for a cell is obtained from the corresponding probe vehicles or connected vehicles, then the speeds of all simulated vehicles located within the cell can be calibrated online by the formula in Equation (16):

$$v_{i,t}^{adj}=v_i^t+\Delta v=v_i^t+\alpha \left(v_{k,t}^{pro{b}_{cell}}-v_i^t\right)$$

(16)

where

$v_{i,t}^{adj}$:

Adjusted speed of vehicle i during time interval t, by FCD or connected vehicle data;

k:

Index denoting a cell in a segment, where vehicle i is located during time interval t;

$v_{k,t}^{prob\_cell}$:

Average speed of cell k, calculated by the corresponding probe vehicle speeds;

α:

Smoothing factor to be calibrated, α∈ [0, 1], and α = 0 if FCD and connected vehicle data are not available for this cell.

As a result, the mesoscopic traffic simulator is able to utilize the adjusted speeds to advance the vehicles in the cells with adequate FCD or connected vehicle data available. This speed adjustment process in Equation (16) is also an on-line calibration approach that may be applied to both ES and GSO processes.

6.5. Simulation Process

In each simulation step, firstly the traffic dynamics parameters are calculated, which include segment input capacities, lane group output capacities, vSIR density and speed for each vehicle in the network, etc. For each lane group, the output counter C_out during the current simulation step is initialized to its output capacity. Then, all the vehicles are moved to their new positions according to the current vehicle speeds and the size of the simulation step. Usually, the step size ranges from 1 to 10 s, the selection of which should be the compromise between the modeling accuracy of traffic dynamics and real-time computational efficiency. The movement of each vehicle in a simulation step is detailed in Figure 10, by which it is convenient to develop simulation programs to process vehicle moving and calculate the newest position of the vehicle in the network at the end of the current simulation step. After the position of each simulated vehicle on each link in the network is updated, the whole simulation process is thus advanced by one step.

If the vehicle is a bus with scheduled route, then it does not have to change the en-route path. If there is some spillback occurred on the immediate downstream link of a vehicle along its current path, then it may change the en-route path to circumvent the congested area. By the modeling approach, some unrealistic deadlocks in the mesoscopic simulator can be avoided, and it seems that people are often doing so to bypass the anticipated traffic jam ahead in the real-life network. If a vehicle moves out of the current link during the simulation step, then it enters the next link on its path and proceeds to move until the total running time T is equal to the step size.

7. Signal Optimization

DynasTIM is able to model signal control scenarios and enable offline or real-time online optimization of the signal timing parameters, which can be converted to the output capacity constraints of lane groups at approaches to an intersection. Usually, the lane group capacity at a signalized intersection is computed as in the following equation:

$$c_i=s_i\frac{g_i}{C}$$

(17)

where s_i and g_i are saturation flow rate and effective green time for lane group i respectively, and C is cycle time. Saturation flow rate is a function of a host of parameters, such as lane width, approach grade, right and left turns, etc. Effective green and cycle times may be computed from signal phasing charts. For unsignalized intersections, lane groups capacities are computed as per the recommendation of HCM 2000 [62].

After completing the calibration of the online traffic simulation system with historical information and real-time surveillance data, the assumption is that the simulator reflects the traffic dynamics in the corresponding real-life network, and has almost the same current and short-term traffic conditions as those of the network. Therefore, the signal control strategies which are evaluated as the most effective by the online simulation system, can be transplanted into the real-life network to optimize the traffic conditions in the field.

The signal control strategies are optimized according to the target described in Wardrop’s second principle [12], that is the sum of the travel times of all travelers in the network is minimum. However, the study in this paper is to approach this target by optimizing signal control strategies rather than the DTA. For the optimization problem, minimizing the total travel time is equivalent to minimizing the total travel delay, and the optimization parameters are the green splits of all signal controllers in the network. Additionally, all signals have common and fixed cycle in the optimization process.

7.1. Optimization Model

The following definitions are adopted:

i:

Subscript denoting a signalized intersection;

j:

Subscript denoting a signal phase;

h:

Subscript denoting a time interval;

q:

Width of prediction window consisting of q time intervals;

w:

Width of optimization window consisting of w time intervals, and 0 ≤ w ≤ q;

d:

d = {d_h, d_{h + 1}, …, d_h + w}, the vector of travel delays for each link in the network from time interval h till time interval (h + w), which are also the measurements of network performance;

d_h:

Vector of travel delays for each link in the network during time interval h which are also the measurements of network performance during time interval h;

h (•):

Simulator function mapping the input parameters to the measurements of network performance (no analytical expression);

D:

Dynamic OD flows for the network, to be calibrated;

β:

Path choice model parameters of the simulator, to be calibrated;

S:

Supply model parameters of the simulator, to be calibrated, including speed-density model parameters of segments, saturation flow rates of lane groups;

nⁱ:

Subscript denoting the number of phases for signalized intersection i;

p_ijh:

Green split for phase j for intersection i during time interval h;

p_ih:

$p_{ih}={\left[p_{i1h} p_{i2h}\cdots p_{i{n}^{i}h}\right]}^T$, vector of green splits for each phase for signalized intersection i during time interval h (superscript T represents vector transpose);

P_h:

P_h = {p_ih}, vector of green splits for each phase for each signal during time interval h;

P:

P = {P_h, P_{h + 1}, …, P_h + w}, vector of green splits for each phase for each signal from time interval h to time interval (h + w);

P_min:

Vector of minimal green splits for each phase for each signal from time interval h to time interval (h + w);

P_max:

Vector of maximum green splits for each phase for each signal from time interval h to time interval (h + w).

With the above definitions, the objective function and constraints are expressed in the following equations:

$$\underset{P}{\min } \Vert \mathbf{d}\Vert$$

(18)

subject to

d = h (D, β, S, P)

(19)

P_min ≤ P ≤ P_max.

(20)

The simulation-based signal optimization problem is solved by the following P-SPSA algorithm.

7.2. Solution Approach

The signal optimization is based on the Parallel Simultaneous Perturbation Stochastic Approximation (P-SPSA) algorithm and simulation based function evaluation to optimize a given network performance indicator. It is designed for stochastic optimization problems and can be applied in the cases that inputs are corrupted by noises, which is especially suitable for the simulation-based signal optimization. In this paper, the parallelization of the gradient calculation is realized and the computational efficiency can be improved multi times, which is particularly crucial for real-time application situations.

In general, the stochastic optimization algorithm produces a sequence of parameter estimators that make the gradient of the objective function gradually converge to zero, where the updates of parameters at iteration k are as follows:

$${\widehat{\theta }}_{k+1}={\widehat{\theta }}_k-a_k\widehat{g}\left({\widehat{\theta }}_k\right)$$

(21)

where ${\widehat{\mathrm{\theta }}}_k$ is the M-dimensional parameter vector to be optimized at the beginning of iteration k; $\widehat{g}\left({\widehat{\mathrm{\theta }}}_k\right)$ is the current estimator of the gradient; $a_k$ is the gain of the step size at iteration k and calculated as follows:

$$a_k=\frac{a}{{\left(A+k\right)}^{\alpha }}$$

(22)

where a, A and α are algorithm parameters.

The SPSA algorithm is proposed by Spall [63,64] and provides huge savings in per iteration cost, through calculating the stochastic gradient with just two function evaluations (independent of the optimization problem size M), and can hence be applied to large-scale and even real-time online optimization situations. The gradient calculation at iteration k is as follows:

$$\widehat{g}\left({\widehat{\theta }}_k\right)=\frac{y\left({\widehat{\theta }}_k+c_k{\Delta }_k\right)-y\left({\widehat{\theta }}_k-c_k{\Delta }_k\right)}{2c_k}\left[\begin{array}{c}\begin{array}{c}{\Delta }_{k1}^{-1}\\ {\Delta }_{k2}^{-1}\\ ⋮\end{array}\\ {\Delta }_{kM}^{-1}\end{array}\right]$$

(23)

where Δ_k is a M-dimensional random perturbation vector, and each component of it is an independent Bernoulli random variable taking the values ±1 with probability 1/2; c_k is usually a small positive number that decreases as k increases; y(•) denotes the objective function value at a given parameter vector expressed by “•”. All the components of $\widehat{g}\left({\widehat{\theta }}_k\right)$ have the common numerator, which indicates the simultaneous perturbation acting on all components of ${\widehat{\theta }}_k$. Therefore, the SPSA algorithm has a M-fold saving in computation efficiency per iteration, compared to the component-wise perturbations and corresponding objective function evaluations in the Finite Difference Stochastic Approximation (FDSA) algorithm. To obtain more stable gradient approximations, especially when objective function evaluations contain random noise, the gradient replications are usually implemented to get the average gradient at each iteration [16,64]. For example, if the number of gradient replications (grad_rep) is three, the number of function evaluations is thus six, each of which can be assigned to an independent computation process to implement parallel function evaluations. This is also the original idea of the P-SPSA algorithm proposed in this paper. In comparison with the standard SPSA algorithm, the speed-up ratio of the P-SPSA algorithm is six theoretically, if the overhead for the process synchronization and data transmission is ignored. The P-SPSA algorithm steps are as follows:

• Initialize the current iteration step k = 0 and the algorithm parameters a = 2.0, A = 50, c = 1.9, α = 0.602, γ = 0.101, according to some suggestions in the literatures [16,64]; also initialize the parameter vector ${\widehat{\theta }}_k={\widehat{\theta }}_0$ that is the column vector of green splits for all signals in this paper.
• Set the number of gradient replications “grad_rep” in order to get the average gradient estimate for ${\widehat{\theta }}_k$ at each iteration step, which is usually set to 3 in our research.
• Increase k = k + 1. Calculate the gains a_k by Equation (22) and c_k = c/k^γ.
• Generate grad_rep M-dimensional random perturbation vectors ${\Delta }_k^p$ (p = 1, 2, …, grad_rep), each component of which is drawn from a Bernoulli ±1 distribution with probability of 1/2 for each ±1 sample.
• Use traffic simulator (e.g., DynasTIM) to perform the parallel objective function evaluations at (grad_rep × 2) points: ${\widehat{\theta }}_k^{p+}={\widehat{\theta }}_k+c_k{\Delta }_k^p$ and ${\widehat{\theta }}_k^{p-}={\widehat{\theta }}_k-c_k{\Delta }_k^p$, where p = 1, 2, …, grad_rep. Note that, before evaluating the objective functions, the constraints in Equation (20) are imposed on each point to ensure that they are within or on the parameter boundaries. Figure 11 shows more implementation details of this step in the MPI parallel computing environments, where np is the total number of computation processes in the MPI parallel domain, tn_p is the number of parameter vectors to be evaluated on each process, and stIdx and endIdx, numbered from 0 through (grad_rep × 2 − 1), are both the indices of parameter vectors in the list of parameter vectors.
• Use Equation (23) to calculate the stochastic gradient vectors $\widehat{g}\left({\widehat{\theta }}_k^p\right)$, p = 1, 2, …, grad_rep, and obtain the average gradient vector as follows:

  $$\widehat{g}\left({\widehat{\mathbf{\theta }}}_k\right)=\frac{{{\displaystyle \sum }}_{p=1}^{grad\_rep}\widehat{g}\left({\widehat{\mathbf{\theta }}}_k^p\right)}{grad\_rep}$$

  (24)
• Update the parameter estimate through Equation (21), which is adjusted again to avoid violations of the boundary constraints in Equation (20).
• Return to Step 3 until convergence or the maximum number of iterations is reached.

8. Software Design

DynasTIM is a complex traffic software system including plenty of models and algorithms. Figure 12 depicts its main functional modules, including network, path, supply, demand, calibration and optimization. The software system has a functional layered architecture and is written by C++ language. It supports the MPI based cluster computing environment, where the whole network can be split into multiple sub networks and loaded by different MPI processes to achieve parallel simulation. In this way, the computation efficiency can be enhanced. Additionally, in each SPSA iteration process, multiple simulation based function evaluations are often needed to obtain more reliable average gradient approximation, which can also be assigned to different MPI processes to perform parallel evaluation. Using a parallel computing technique to improve operating efficiency is crucial for the applications with strict time limit, such as real-time signal control strategy optimization. In addition, for the management of large number of historical and real-time surveillance data and simulation outputs, the software adopts the relational database management systems including MySQL, SQL Server and Access; and for the network data, path data, etc., because of the obvious hierarchical features, it uses structured text files (e.g., XML files) to store them.

9. Case Study

9.1. Network and Data

The network and data come from a cooperative project between Shenzhen Urban Transport Planning Center (SUTPC) and the author. The main purpose of this project is to study whether the existing traffic simulation technology can accurately reproduce the real traffic states for complex urban road networks with rich and accurate flow detection data.

This case study is based on a populated network from the Futian Central Business District (CBD) in Shenzhen, China, which is depicted in Figure 13. There are plenty of financial, commercial, cultural, service buildings, hotels and apartments etc., in the network. The study network comprises urban expressways, auxiliary roads, arteries, minor roads and branch roads. The computer representation of the network consists of 795 nodes connected by 1301 directed links (shown in Figure 14) in an area of around 7 square kilometers, and the numbers of different network elements are shown in

Table 1. Numbers of different network elements.

Zones	Nodes	Links	Segments	Lane Groups	Lanes	Signals	Sensors
134	795	1301	1426	1726	3172	38	359

.

Surveillance data in the case study, provided by SUTPC, include the generation and attraction volumes of main traffic analysis zones from one week during August 2016 for 24 h, and traffic counts from one weekday during November 2016 between the hours of 7:30 a.m. and 9:30 a.m. The traffic counts were obtained from 359 video detection locations, where the turning flows were recorded. All of the data were aggregated in 15 min intervals. There are eight estimation time intervals in total, each of which is 15 min long and ranges from 7:30 a.m. to 9:30 a.m. The network includes 38 signalized intersections with pre-timed signal control plans, which are also provided by SUTPC. In addition, SUTPC provided the k-v model parameters and the saturation flow data for different types of roads.

9.2. Calibration

In this paper, the Normalized Root Mean Square Error (NRMSE) statistic is used to evaluate the quality of the calibration results, which represents the discrepancy between the observed and the simulated traffic counts, and is computed by the formula in Equation (25).

$$NRMSE=\frac{\sqrt{({{\displaystyle \sum }}_{i=1}^S{\left(y_i-{\widehat{y}}_i\right)}^2)/S}}{\sqrt{({{\displaystyle \sum }}_{i=1}^S{y}_i^2)/S}}=\frac{\sqrt{{{\displaystyle \sum }}_{i=1}^S{\left(y_i-{\widehat{y}}_i\right)}^2}}{\sqrt{{{\displaystyle \sum }}_{i=1}^S{y}_i^2}}$$

(25)

where

$y_i$:

vector of observed traffic counts in the field;

${\widehat{y}}_i$:

vector of simulated traffic counts generated by the system;

S:

dimension of the two vectors.

A lower value of NRMSE indicates that the discrepancy between the simulation outputs and the observed counterparts is smaller.

Additionally, the mesoscopic traffic model parameter values are provided by SUTPC, which are based on their traffic simulation projects and calibrated roughly for different types of roads. The parameters for vSIR length calculation are set as follows: resp_time = 2.5 s, and μ = 2.0. The set of OD paths is generated by calling the k-shortest path (KSP) algorithm after adding random perturbations to the link impedances. The Path-size Logit model, which accounts for overlapping routes, is used to model the route choice behaviors in the network.

Based on the generation and attraction volumes of zones, the gravity model is applied in the calculation of the OD flows just for the first estimation interval (7:30–7:45), and a total of 8398 OD pairs are found. What follows is to use the RLS based dynamic OD estimation algorithm to obtain the apriori dynamic OD flows for all of the estimation intervals. Since most of the apriori OD flows are fairly small, the OD pairs with the largest OD flow less than three are removed in all estimation intervals, and finally a total of 1319 major OD pairs are obtained. Then, based on the apriori dynamic OD flows estimated by the RLS algorithm, the SparseGLS algorithm and the simulation-based assignment matrix calculation method are used to estimate the dynamic OD flows again, in order to improve the estimation accuracy further. The NRMSE results of two different OD estimation algorithms are given in

Table 2. Estimation results.

Algorithm	No. of Observations	NRMSE
Dynamic OD flow estimation with RLS algorithm	2872	0.1854
Dynamic OD flow estimation with SparseGLS algorithm	2872	0.1388

, and Figure 15 depicts the estimation results for the comparison between the simulation outputs with corresponding real-life count measurements during the morning peak from 7:30 to 9:30 a.m. The x-axis represents the detected traffic counts (including turning direction counts), and the y-axis represents the estimated traffic counts (the simulation outputs). The circles exactly located on the 45° line indicate that the estimated counts by the system are perfectly same as the corresponding measurements in the field. Compared with the RLS algorithm, the SparseGLS one increases the estimation accuracy of the traffic counts by about 25%. As the surveillance data including FCD get richer in the future, the estimation accuracy of mesoscopic model parameters, OD path set, and the route choice model parameters, will be improved further, thus, the better calibration results can be acquired. It should be noted that the flow estimation error in this case study is obviously lower than the experimental result of Transmodeler [15], which represents the state of the art of traffic simulation, and the estimation error of the latter should be more than 30%.

9.3. Urban Area Signal Optimization

Currently, the predetermined time-of-day signal control plans are applied in the Futian CBD network, where the two-hour morning peak is usually divided into 2–3 time intervals, and the signal plans are preset for each of them respectively. These plans have been continuously optimized by the Shenzhen TMC for many years, so as to be more adaptable to the dynamic traffic flows. One question to be answered is whether the existing signal plans can be further optimized and how much travel delay can be reduced. Therefore, the SPSA based algorithm established in this paper is utilized to optimize the timing plans of the 38 signal controllers in the network simultaneously. In this case study, each signalized intersection uses only one timing plan with the same signal cycle of 150 s, and the optimization is only for the green splits of all signals in the network. The setting of the cycle time is the result of many simulation tests and taking into account the safety and acceptable waiting time of pedestrians and non-motor vehicles.

Figure 16 shows the SPSA based signal optimization with 40 iterations. It can be seen the objective function value fluctuates and does not always decrease at some iterations, since this is a stochastic optimization problem and the calculated gradient vector is a stochastic approximation. However, the signal timing parameters and the corresponding objective function value in each iteration can be stored into a data pool. After the completion of the optimization process, the signal timing parameters with minimal objective function value can be extracted from the data pool as the final signal scenario.

Table 3. Network MOEs before and after signal optimization.

	Current Signal Plans	Optimized Signal Plans	Improvement (%)
No. of vehicles entering network	162,434	164,441	1.24
No. of vehicles leaving network	144,375	148,287	2.71
Average vehicle travel delay (s)	37.1	32.2	13.21
Total travel time (h)	23,603.9	21,370.8	9.46
Total travel distance (km)	312,744	316,815	1.30
Maximum node delay (h)	1436.5	1221.6	14.96
Network average speed (km/h)	13.2	14.8	12.12

gives the simulation result comparison between the currently implemented signal timing plans in the field and the plans optimized by the signal optimization algorithm proposed in this paper. It shows that the optimized signal plans can effectively reduce the travel delay and increase the average vehicle speed for the entire network while maintaining the slight improvement of the network throughput.

In comparison with the signal optimization methods proposed by Carolina et al. [50,51,52,53], the P-SPSA based optimization method requires fewer simulation evaluations, so it has more potential for real-time application. It does not need a metamodel as a bridge between the traffic simulator and the optimization algorithm, thus reducing the modeling cost. Meanwhile, even for the large-scale signal optimization problem for the congested CBD road network in Shenzhen, it seems that the algorithm can still produce fairly stable and reliable optimization result.

9.4. Computational Performance

To investigate the computational performance of the algorithms presented in this paper, the experiments are implemented on a Thinkpad T470p laptop with Intel Core i5-7300HQ Processor (up to 3.50GHz), 8 GB memory, 500GB Hard Disk, and Windows 10 Pro 64 operating system. As shown in Figure 17, the running times of the algorithms are compared with the duration of the analysis time period (2 h morning peak from 7:30 to 9:30). Only the algorithms of the running time less than 2 h have the potential to be deployed in real time. In this paper, the RLS algorithm is applied in estimating the OD flows for 8398 OD pairs in the network, then, the major OD flows for 1319 OD pairs are estimated by the SparseGLS algorithm, the running time of which only account for around 14% that of RLS algorithm. The computational performance of the standard SPSA algorithm [16,64] and the P-SPSA algorithm is compared, both of which perform 40 iterations. The P-SPSA algorithm initiates three computation processes (assigned to three CPU cores) to accelerate the objective function evaluation and as a result saves about 64% of running time. Thus, the parallel speedup ratio is around 2.78 with three CPU cores in this case, and theoretically it should be close to 6.0 with six computation processes assigned to six CPU cores. Supported by the parallel computation environment, the P-SPSA algorithm appears to be capable of the real-time optimization of signal timing parameters for fairly large network like the Futian CBD in this case study.

10. Conclusions and Future Directions

A real-time online simulation system, DynasTIM, used for traffic estimation, prediction and optimization, is presented. Through the effective integration of demand-and-supply modeling and simulation, and optimization techniques, the system can provide online decision support for ITS to help improve the performance of signal control and travel guidance strategies based on the predicted traffic states and network performance evaluation. The case study in the city of Shenzhen, China, shows that the system can estimate the traffic counts fairly well, and the signal optimization method generates better signal timing plans compared with those currently performed in the actual network.

An important feature of DynasTIM is, after completion of the state prediction, it is not waiting for the next estimation interval of surveillance data available before moving to the next state estimation cycle, but immediately begins receiving the field data to implement online calibration again (recall that the first online calibration is in the ES process), and then continues to simulate vehicle movements in the network. In this way, the impact of incidents and sudden travel demand previously unforeseen on the network performance, can be responded timely to determine whether to update the guidance strategy correspondingly. As a result, for real-time online applications, DynasTIM is always running without idle time. Additionally, an analytical assignment matrix computation formulation is proposed and proved for dynamic OD flows estimation, and the simulation based signal optimization framework is established based on the high efficiency P-SPSA algorithm. Although a variety of traffic models have been developed to describe real-life traffic dynamics characteristics, if actual link travel speeds can be accurately detected (e.g., from probe vehicles), it seems more reasonable that these speed data are directly exploited to calculate vehicle displacements in the mesoscopic traffic simulator instead of relying on the mesoscopic traffic models to calculate vehicle speeds and corresponding displacements based on traffic densities. For this reason, a speed adjustment method is proposed for the vehicles in the moving part of a segment. Indeed, the mesoscopic traffic models is used to calculate vehicle speeds for the links without any detected speed data available or for traffic prediction process.

In future work, it will be valuable to collect more abundant field data to better calibrate the parameters of mesoscopic traffic models and route choice models, and to also help refine and improve the models further . Moreover, the online calibration of real-time DTA systems is still an important aspect that needs more research. In addition, the development of hierarchical adaptive signal control algorithms based on the simulation optimization framework is also an interesting research direction.