Gramian-Based Optimization for the Analysis and Control of Traffic Networks

Gramian-Based Optimization for the Analysis and Control of Traffic Networks

Gianluca Bianchin,  and Fabio Pasqualetti,  This material is based upon work supported in part by ARO award 71603NSYIP, and in part by NSF award CNS1646641. The authors are with the Department of Mechanical Engineering, University of California, Riverside, {gianluca, fabiopas}

This paper proposes a simplified version of classical models for urban transportation networks, and studies the problem of controlling intersections with the goal of optimizing network-wide congestion. Differently from traditional approaches to control traffic signaling, a simplified framework allows for a more tractable analysis of the network overall dynamics, and enables the design of critical parameters while considering network-wide measures of efficiency. Motivated by the increasing availability of real-time high-resolution traffic data, we cast an optimization problem that formalizes the goal of minimizing the overall network congestion by optimally controlling the durations of green lights at intersections. Our formulation allows us to relate congestion objectives with the problem of optimizing a metric of controllability of an associated dynamical network. We then provide a technique to efficiently solve the optimization by parallelizing the computation among a group of distributed agents. Lastly, we assess the benefits of the proposed modeling and optimization framework through microscopic simulations on typical traffic commute scenarios for the area of Manhattan. The optimization framework proposed in this study is made available online on a Sumo microscopic simulator based interface [1].

I Introduction

Effective control of transportation systems is at the core of the smart city paradigm, and has the potential for improving efficiency and reliability of urban mobility. Modern urban transportation architectures comprise two fundamental components: traffic intersections and interconnecting roads. Intersections connect and regulate conflicting traffic flows among adjacent roads, and their effective control can sensibly improve travel time and prevent congestion. Congestion is the result of networks operating close to their capacity, and often leads to degraded throughput and increased travel time.

The increasing availability of sensors for vehicle detection and flow estimation, combined with modern communication capabilities (e.g. vehicle-to-vehicle (V2V) and vehicle-to-infrastructure (V2I) communication), have inspired the development of infrastructure control algorithms that are adaptive [2], that is, policies that adjust the operation of the system based on the current traffic conditions. Nevertheless, the remarkable complexity of modern urban transportation infrastructures has recently promoted the diffusion of control policies that are distributed, that is, algorithms that adapt the operation of individual network components based on partial knowledge of the current network state and dynamics. The lack of a global network model, capable of capturing the interactions between spatially-distributed components and capable of modeling all the relevant network dynamics, often results in suboptimal performance [3]. In this paper, we propose a simplified model to capture the time and spatial relationships between traffic flows in urban traffic networks in regimes of free-flow. The model represents a tradeoff between accuracy and tractability, and sets out as a tractable framework to study efficiency and reliability of this class of dynamical systems. The proposed framework is employed in this work for the control of green split times at the signalized intersections.

Related Work: The design of feedback policies for the control of urban infrastructures is an intensively studied topic, and the proposed techniques can mainly be divided into two categories: routing policies and intersections control. Routing policies use a combination of turning preferences and speed limits in order to optimize congestion objectives, and have been studied both in a centralized [4] and distributed [5] framework. Conversely, intersection control refers to the design of the scheduling of the (automated) intersections so that the flow through intersections is maximized, and can be achieved (i) by controlling the signaling sequence and offset, and/or (ii) by designing the durations of the signaling phases. The control of signals offset typically aims at tuning the synchronization of green lights between adjacent intersections in order to produce green-wave effects [6, 7], and consists of solving a group of optimization problems that take into account certain subparts of the infrastructure, while minimizing metrics such as the number of stops experienced by the vehicles. In contrast, the durations of green time at intersections impact the behavior of certain traffic flows within the network, and plays a significant role in the efficiency of large-scale dynamical networks [3].

Widely-used distributed signaling control methods include SCOOT [8], RHODES [9], OPAC [10], and emerge as the most common techniques currently employed in major cities. The sub-optimal performance of the above methods [11] has motivated the development of max-pressure techniques [2]. Max-pressure methods use a discrete-time model where queues at intersections have unlimited queue lengths. Under this assumption, max pressure is proven to maximize the throughput by stabilizing the network. Centralized policies require higher modeling efforts but, in general, have better performance guarantees [12]. Among the centralized policies, the Traffic-Responsive Urban Control framework [3] has received considerable interest for its simplicity and good performance. Based on a store-and-forward modeling paradigm, the method consists of optimizing network queue lengths through a linear-quadratic regulator problem that uses a relaxation of the physical constraints to abide with the high complexity. Variations of these techniques to incorporate physical constraints have been studied in [13, 14]. The increased complexity of urban traffic networks has recently motivated the development of simplified (averaged) models to deal with the switching nature of the traffic signals [15]. However, the highly-nonlinear behavior of this class of dynamical systems still limits our capability to consider adequate optimization and prediction horizons [16], and the development of tractable models capable of capturing all the relevant network dynamics is still an open problem.

Contribution: Motivated by the considerable complexity of traditional models for urban transportation systems, in this paper we propose a simplified framework to capture the behavior of traffic networks that are operating close to the free-flow regime. In this model, each road is associated with multiple state-variables, representing the spatial evolution of traffic densities within the road. This assumption allows us to capture the non-uniform spatial displacement of traffic within each road, and to construct a simplified network model that results in a more-tractable framework for optimization.

We employ the proposed model to formalize the goal of optimally designing the durations of green splits at intersections, and we provide a connection with the optimization of a metric of controllability for a dynamical model associated with the traffic network. To the best of the authors’ knowledge, this work represents a novel, computationally-tractable, method to perform network-wide optimization of the green-splits durations at intersections. We provide conditions that guarantee network stability, and we characterize the performance of the system in relation to the overall network congestion. We use the concept of smoothed spectral abscissa [17] to numerically solve the optimization, and we demonstrate the benefits of our methods through of a microscopic simulator on the area of Manhattan. We characterize the complexity of our algorithms, and propose a method to parallelize the computation in order to be more-efficiently executed by a group of distributed cooperating agents. Our results and simulations suggest that the increased system performance and stability obtained by the methods justify the increment in complexity deriving from a network-wide model description.

Organization: The rest of this paper is organized as follows. Section II presents the problem setup, relates the underlying assumptions to previously-established traffic models, and formalizes the goal of minimizing an overall measure of network congestion by selecting the durations of the green split times at the intersections. Section III discusses our approach to solve the resulting nonlinear optimization problem. The solution method is at first presented from a centralized perspective, and subsequently adapted for implementation over a distributed architecture (Section IV). In the distributed paradigm, the optimization shall be solved by a group of cooperating agents, each with limited (local) knowledge of the physical network interconnection. Section V is devoted to numerical simulations that validate our assumptions and the solution technique. Section VI concludes the paper.

Ii Dynamical Model of Traffic Networks and Problem Formulation

We model urban traffic networks as a group of one-way roads interconnected through signalized intersections. Within each road, vehicles move at the free-flow velocity, while traffic flows are exchanged between adjacent roads by means of the signalized intersection connecting them. In this section, we discuss a concise dynamical model for traffic networks in regimes of free flow, that will be employed for the analysis.

Ii-a Model of Road and Traffic Flow

Let denote a traffic network with roads and intersections . Each element in the set models a one way road, whereas intersections regulate conflicting flows of traffic among adjacent roads (see Section II-B). We assume that exogenous inflows enter the network at (source) roads and, similarly, vehicles exit the network at (destination) roads , with . The following standard connectivity assumption ensures that vehicles are allowed to leave the network.

Assumption 1.

For every road there exists at least one path in from to a road .

We denote by the length of road , and we model each road by discretizing it into cells (see Fig. 1), where the parameter is a constant discretization step.

Fig. 1: Road discretization and corresponding network model. The portion of road comprised between spatial coordinates and , is referred to as cell , and variable denotes the cell density.

We denote by the traffic density111The terms density and occupancy will be used interchangeably in the remainder of this paper. associated with the -th cell of road , . We assume that inflows of vehicles enter the road in correspondence of its upstream cell (i.e. ); accordingly, outflows leave the road in correspondence of its downstream cell (i.e. ). We model the relation between traffic flows and cell density by assuming that flows of vehicles move in regimes of free flow with constants velocity , which represents the average speed along link . Then, the dynamics of the road state are described by:


Differently from conventional dynamical models for urban links (e.g. [18]), the space discretization technique (1) allows us to capture the fact that the density of vehicles may not be uniform along the road. Notably, this feature plays a key role in modeling the road outflows in correspondence of signalized intersections (see Section II-B).

Fig. 2: Fundamental diagram describing the static speed-density relation. The almost-flat behavior in regimes of free-flow and heavy congestion allows us to approximate , as the speed of the of flow is not sensibly affected by changes in density.
Remark 1.

(Equivalence Between (1) and Hydrodynamic Models) The dynamical model (1) derives from the mass-conservation continuity equation [19] in certain traffic regimes, as we explain next. Let the function denote the (continuous) density of vehicles within road at the spatial coordinate and time . Let denote the (continuous) flow of vehicles along the road, and let traffic densities and flows follow the hydrodynamic relation

We first complement the above continuity equation with the Lighthill-Whitham-Richards static relation, , in which traffic flows instantaneously change with the density. Then, we include the speed-density fundamental relation , where represents the speed of the traffic flow (see Fig. 2), to obtain

Solutions to the above relation are kinematic waves [20] moving at speed . We consider regimes of free flow where the speed of the kinematic wave can be approximated as . As illustrated in Fig. 2, this approximation is accurate in regimes of free flow or congestion, characterized by . Therefore, we let denote the average speed of the flow along the link, and consider the approximated continuity equation

We then discretize in space the above linear continuity equation, by defining the discrete spatial coordinate

and by replacing the partial derivative with respect to with the difference quotient

This discretization leads to the dynamical model (1) after introducing the boundaries conditions, represented by inflows and outflows , and by replacing the spatially discretized density with the functions of time .

Ii-B Model of Intersection and Interconnection Flow

Signalized intersections alternate the right-of-way of vehicles to coordinate and secure conflicting flows between adjacent roads. Every signalized intersection , , is modeled as a set , consisting of all allowed movements between the intersecting roads. For road , let denote the (unique) intersection at the road upstream; similarly, let denote the (unique) intersection at the road downstream. We model the effect of signalized intersections through a set of green split functions that assume boolean values (green phase) or (red phase), and let the road inflows be


where denotes the intersection transmission rate. It is worth noting that the notation represents the transmission rate from road to and, similarly, denotes the green split function that controls traffic flows from and directed to . In (2), we incorporated exogenous inflows and outflows to roads (flows that are not originated or merge to modeled intersections or roads) into and respectively, by defining the input term and output term . We note that if and only if , and if and only if .

Fig. 3: Typical set of phases at a four-ways intersection.
Fig. 4: Green splits for the four ways intersection in Fig. 3.
Example 1.

(Example of Intersections and Scheduling Functions) Consider the four-ways intersection illustrated in Fig. 3. The intersection is modeled through the set of allowed movements:

Intersections are often partitioned into a group of phases, where each phase represents a set of movements that can occur simultaneously across the intersection. For intersection , a typical set of phases is

The green split function is then defined by alternating the available phases within a cycle time , that is, for all ,

where are the switching times. See Fig. 4 for a graphical illustration.

Next, we make critical use of the underlying road-discretization technique (1), and model transmission rates as functions proportional to the occupancy of the cell at downstream of each road, that is,


where is a parameter that models the speed of the outflow, and includes the average routing ratio of vehicles entering road from . Although (3) constitutes an approximation of the intersection transmission rate, which often contains saturation terms (e.g. see [21] for a discussion), the combination of (3) with the proposed space-discretization technique (1) performs well in practice (see Section V for numerical validation of (3)), at least in the considered regimes.

Remark 2.

(Turning Rates and Conservation of Flows at Intersections) The existence of independent lanes for different turning preferences in proximity of intersections often leads to the decomposition , where the function represents the constant average routing ratio of vehicles entering road from . The conservation of flows at intersections can be guaranteed by letting .

We emphasize that the proposed model represents a simplified framework with respect to traditional, more-established, models. In fact, the assumptions we make are realistic for systems that are operating close to their free-flow regime, and do not take into account collateral phenomena, such as back-propagation in regimes of congestion, which would require ad-hoc adjustments in the model.

Ii-C Switching and Time-Invariant Traffic Network Dynamics

Individual road dynamics can be combined into a joint network model that captures the interactions among all modeled routes and intersections. To this aim, we adopt an approach similar to [14], and assume that exogenous outflows are proportional to the number of vehicles in the road, that is, , . By combining Equations (1), (2) and (3), we obtain


where , is the overall number of modeled network cells, derives from (2), and

where denotes the -th canonical vector of size .

We note that the matrix in (4) is typically sparse, because not all roads are adjacent in the network interconnection, and its sparsity pattern varies over time as determined by the green splits . Thus, the network model (4) represents a linear switching system, where the switching signals are the green split functions.

Fig. 5: Network model associated with a traffic network composed of intersections and roads. Each road is associated with a set of states that represent the density of the cells within the roads.
Example 2.

(Example of Traffic Network) Consider the network illustrated in Fig. 5, with and . The network comprises four destination roads ( for all , and otherwise), and four source roads (, with only if ). Let and for all . Then, the matrices in (4) read as

for all . Notice that for all times if for all .

Next, we make the typical assumption that scheduling functions are periodic, with period . That is, for all , , and for all times :

Let denote the set of time instants when a scheduling function changes its value, that is,

Notice that the network matrix A in (4) remains constant between any consecutive time instants and . We denote each constant matrix by , and refer to it as to the -th network mode. Further, let , with and , denote the duration of the -th network mode. We employ a state-space averaging technique [22] and define a linear-time invariant approximation of the switching network model (4):


where , and , . We note that the averaging technique preserves the sparsity pattern of the network, that is, if and only if for some .

In general, the approximation of the behavior of the switching system (4) with the average dynamics (5) is accurate if the operating period is short in comparison to the underlying system dynamics. Under suitable technical assumptions, the deviation of average models with respect to the network instantaneous state has been characterized in e.g. [15]. In [15], the authors formally provide a bound on the accuracy of the approximation, where the bound becomes tighter for decreasing values of , and increasing values of road lengths. Although a formal characterization of the performance of average model is often an extremely challenging task, and the formal results in [15] hold for a single roads in the presence of a single signalized intersections, the accuracy of averaging techniques has been shown to be satisfactory in similar applications (see e.g. [16]). A numerical validation of the averaging technique, and its validity in relation to the cycle-time , is discussed in Section V (see Fig. 7).

Ii-D Problem Formulation

In this paper, we consider the dynamical model (5) and focus on the problem of designing the durations of the green split functions so that a measure of network congestion is optimized. As illustrated through the relation (see (5)), the average model allows us to design the durations of the network modes, rather than their exact sequence. This approach motivates the adoption of a two-stage optimization process. First, the durations of the modes is optimized by considering a joint model that captures the dynamics of the entire interconnection. Second, offset control techniques (see e.g. [7]) can be employed to decide the specific sequence of phases, given the durations of the splits and by considering local interconnection models. In this paper, we focus on the first stage. To formalize our optimization problem, we denote by the vector of the queue lengths originated by the signalized intersections, and model as the occupancy of the cells at roads downstream, that is,


We assume the network is initially at a certain initial state , and focus on the problem of optimally designing the mode durations that minimize the -norm of the vector of queue lengths . Namely, we consider the following dynamical optimization problem:

subject to (7a)

Loosely speaking, the optimization problem (7) seeks for an optimal set of split durations that minimize the -norm of the impulse-response of the system to the initial conditions . Thus, similarly to [23], our framework considers the “cool down” period. In this scenario, exogenous inflows and outflows are not known a priori, and the goal is to evacuate the network as fast as possible, where the final condition is an empty system. Under this assumption, and the solution to (7), shall be re-computed when the current traffic conditions change significantly. Finally, we note that constraint (7e) guarantees feasibility of the resulting solutions. In fact, the time-normalization in (7e) implies that for any solution resulting from (7), there exists (at least) one sequence of movements at each intersection, with overall cycle time , and that abide the selected set of durations.

Iii Design of Optimal Network Mode Durations

In this section, we propose a numerical method to determine solutions to the optimization problem (7). The approach we discuss is centralized, namely, it requires knowledge of the state for all network cells. An extension of the framework to fit a distributed implementations is then discussed in Section IV. At its core, the technique relies on rewriting the cost function of (7) in relation to the controllability Gramian of an appropriately-defined linear system, as we explain next.

Lemma 2.

(Controllability Gramian Cost Function) Let

The following minimization problem is equivalent to (7):

subject to

By incorporating (7a), (7b), and (7c) into the cost function of optimization problem (7), we can rewrite:

from which the claimed statement follows. ∎

Next, we note that the cost function of (7) is finite only if the choice of parameters leads to a dynamical matrix that is Hurwitz-stable. Requiring Hurwitz-stability of corresponds to imposing , where denotes the spectral abscissa of . The following result proves network stability under optimal sets of split durations.

Theorem 3.

(Stability of Optimal Solutions) Consider a traffic network , and assume that there exists a path from every node in to at least one node in . Let for all , and for some . Then,


From the structure of (4) and from the assumption follows that for all , while for all , . Moreover, all columns of have nonpositive sum. In particular, the columns corresponding to destination cells have strictly negative sum, that is, for all , and for all such that . To show , we use the fact that destination cells in have no departing edges, and re-order the states so that

where , , is the submatrix that describes the dynamics of the destination cells, , and . The fact immediately follows from (2). The stability of follows from the connectivity assumption in the original network, and from the analysis of grounded Laplacian matrices (see e.g. [24, Theorem 1]). ∎

For the solution of the optimization problem (2), we propose a method based on a generalization of the concept of spectral abscissa, namely, the smoothed spectral abscissa. For a dynamical system of the form (5)-(6), the smoothed spectral abscissa [17] is defined as the root of the implicit equation


where . It is worth noting that the root is unique [17], and for fixed and it is a function of both and , namely .

Remark 3.

(Properties of the Smoothed Spectral Abscissa) For any , the smoothed spectral abscissa is an upper bound to , and this bound becomes exact as . To see this, we first observe that the integral exists and is finite for any , as the function is bounded and convergent as . On the other hand, for any the function becomes unbounded for and the above integral is infinite. It follows that, the left-hand side of (9) is finite only if or, in other words, for any finite , satisfies .

We observe that, by letting in (9), the optimization problem (2) can equivalently be reformulated in terms of the smoothed spectral abscissa as follows:

subject to

where the parameter is now an optimization variable. In what follows, we denote with the value of the optimization parameters at optimality of (III). Problem (III) is a nonlinear optimization problem [17], because the optimization variables and are related by means of the nonlinear equation (9).

For the solution of (III), we propose an iterative two-stages numerical optimization process. In the first stage, we fix the value of and seek for a choice of that leads to a smoothed spectral abscissa that is identically zero. In other words, we let , and solve the following minimization problem:

subject to

It is worth noting that every that is solution to (III) and that satisfies , is a point in the feasible set of (III), which corresponds to a cost of congestion .

In the second stage, we perform a line-search over the parameter , where the value of is increased at every iteration until the minimizer is achieved. We note that the optimizer of (III) with is indeed , an optimal solution to (III). Thus, heuristically, the progressive -update step in the two stages optimization problem allows us to determine (local) solutions to (III).

The remainder of this section is devoted to the description of the two-stage optimization process. The benefit of solving (III) as opposed to (III) is that we can derive an expression for the gradient of with respect to the mode durations , as we illustrate next. In the remainder of this section, with a slight abuse of notation, we use the compact form and, for a matrix , we denote its vectorization by .

Lemma 4.

(Descent Direction) Let denote the unique root of (9) with . Let , and let . Then,

where and are the unique solution to the two Lyapunov equations


and denotes the identity matrix.


The expression for the partial derivative of the smoothed spectral abscissa with respect to can be obtained from the composite function

where follows immediately from (7d), and the expression for the derivative of with respect to is given in [17, Theorem 3.2]. ∎

We note that equations (4) always admit unique solution. To see this, we use the fact that is an upper bound to , and observe that is Hurwitz-stable for every .

A gradient descent method based on Lemma 4 is illustrated in Algorithm 1. Each iteration of the algorithm comprises the following steps. First, (lines ) a (possibly non feasible) descent direction is derived as illustrated in Lemma 4. Second, (line ) a gradient-projection technique [25] is used to enforce constraints (7d)-(7f). The update-step follows (line ).

Input: Matrix , vector , scalars ,
Output: solution to (7)
1 Initialize: , ,
2 while  do
3       repeat
4             Compute by solving (9);
5             Solve for and : ; ;
6             ;
7             ;
8             Compute projection matrix ;
9             ;
10             ;
11             ;
13      until ;
14      ;
16 end while
17return ;
Algorithm 1 Centralized solution to (7).

Algorithm 1 employs a fixed stepsize , and a terminating criterion (line ) based on the Karush-Kuhn-Tucker conditions for projection methods [25]. The -update step, which constitutes the outer while-loop (line ), is then performed at each iteration of the gradient descent phase, and the line-search is terminated when cannot be achieved. To prevent the algorithm from stopping at local minimas, the gradient descent algorithm can be repeated over multiple feasible initial conditions .

Remark 4.

(Complexity of Algorithm 1) The computational complexity of Algorithm 1 in relation to the network size can be derived as follows. First, the computation of the current value of the smoothed spectral abscissa can be achieved via a root-finding algorithm (such as the bisection algorithm) over the implicit function (9). The complexity of this operation is a logarithmic function of the desired accuracy, and it requires the computation of the trace of the Gramian matrix. Thus, for given accuracy, the complexity of this operation is . Second, Algorithm 1 requires to determine solutions to the pair of Lyapunov equations (4). Modern methods to solve Lyapunov equations [26] rely on the Schur decomposition of the matrix , whose complexity is . It is worth noting that, given the Schur decomposition , where is upper triangular and is unitary, a decomposition for follows immediately by shifting to . Therefore, a single decomposition is required at each iteration of the gradient descent and the complexity of Algorithm 1 is therefore .

Iv Distributed Gradient Descent

The centralized computation of assumes complete knowledge of the matrices , and requires to numerically solve the pair of Lyapunov equations (4). For large-scale traffic networks, such computation imposes a limitation in the dimension of the matrix and, consequently, on the number of signalized intersections that can be optimized simultaneously. Since the performance of the proposed optimization technique depends upon the possibility of modeling and optimizing large network interconnections, a limitation on the number of modeled roads and intersections constitutes a bottleneck toward the development of more efficient traffic infrastructures. A possible solution to address this complexity issue is to distribute the computation of the descent direction in Algorithm 1 among a group of agents, in a way that each agent is responsible for a subpart of the entire computation (e.g. see Fig. 6). Agents can represent geographically-distributed control centers or clusters in parallel computing, each responsible for the control of a group of intersections.

In order to distribute the computation of solutions to (4), in this section we focus on distributively solving Lyapunov equations of the form


where is an unknown matrix, , . Let be partitioned as

where , . We assume that each agent knows only. Note that are sparse matrices, and their sparsity pattern depends upon the subpart of infrastructure associated with that agent. In addition, we assume that neighboring agents are allowed to exchange information by means of a communication infrastructure. Let be the communication graph, where each vertex represents one agent, and represents the interconnection. The following lemma constitutes the key ingredient of the method we propose to distributively compute .

Lemma 5.

(Distributed Solutions to (13)) Consider a Lyapunov equation of the form (13), and let be Hurwitz-stable. The following statements are equivalent:

  1. solves (13);