Samitha Samaranayake
email: samitha@alum.mit.edu

Jack Reilly

Walid Krichene

Maria Laura

Delle Monache

Paola Goatin

Alexandre Bayen

Discrete-time system optimal dynamic traffic assignment (SO-DTA) with partial control for horizontal queuing networks

We consider the System Optimal Dynamic Traffic Assignment problem with Partial Control (SO-DTA-PC) for general networks with horizontal queuing. The goal of which is to optimally control any subset of the networks agents to minimize the total congestion of all agents in the network. We adopt a flow dynamics model that is a Godunov discretization of the Lighthill-Williams-Richards (LWR) partial differential equation with a triangular flux function and a corresponding multi-commodity junction solver. Full Lagrangian paths are assumed to be known for the controllable agents, while we only assume knowledge of the aggregate turn ratios for the non-controllable (selfish) agents. We solve the resulting finite horizon non-linear optimal control problem using the discrete adjoint method.

Introduction

Dynamic traffic assignment (DTA) is the process of allocating time-varying origin-destination (OD) based traffic demand to a set of paths on a road network. This problem has been studied extensively over the last 35 years, since the seminal works of Merchant and Nemhauser (1978a,b). There are two types of traffic assignment that are generally considered, the user equilibrium or Wardrop equilibrium allocation (UE-DTA), in which users minimize individual travel-time in a selfish manner, and the system optimal allocation (SO-DTA) where a central authority picks the route for each user and seeks to minimize the aggregate total travel-time over all users. These principles were first presented in the work of [START_REF] Wardrop | Some theoretical aspects of road traffic research[END_REF] in the context of static traffic assignment and expanded on by [START_REF] Beckman | Studies in the Economics of Transportation[END_REF]. See [START_REF] Peeta | Foundations of Dynamic Traffic Assignment : The Past , the Present and the Future[END_REF] for a broad overview of dynamic traffic assignment. User equilibrium (UE) traffic assignment can lead to inefficient network utilization, highlighted by Braess' Paradox [START_REF] Braess | Über ein Paradoxon aus der Verkehrsplanung[END_REF], where adding capacity to the network can actually result in longer travel times for all users. It has been shown that this paradox can occur in real road networks [START_REF] Kelly | Network routing[END_REF] and that it is hard to design networks that are immune to it [START_REF] Roughgarden | On the severity of Braess's paradox: designing networks for selfish users is hard[END_REF]. In fact, it can be shown that the price of anarchy (PoA) [START_REF] Koutsoupias | Worst-case equilibria[END_REF], the worst-case ratio of the system delay caused by the selfish behavior over the system optimal solution, may be arbitrarily large even in simple networks [START_REF] Roughgarden | Bounding the inefficiency of equilibria in nonatomic congestion games[END_REF][START_REF] Swamy | The effectiveness of Stackelberg strategies and tolls for network congestion games[END_REF]. System optimal (SO) traffic assignment on the other hand leads to optimal utilization of the network resources, but is hard to achieve in practice since the overriding objective for individual drivers in a road network is to minimize their own travel-time. It is well known that setting a toll on each road segment corresponding to the marginal delay of the demand moves the user equilibrium towards a SO allocation [START_REF] Roughgarden | The Price of Anarchy is Independent of the Network[END_REF][START_REF] Vickrey | Congestion theory and transport investment[END_REF]. However, imposing time-varying tolls on each road segment is impractical and tolling in general is difficult to implement in many settings due to both infrastructure and political considerations.

An alternative approach is to attempt to control a fraction of the drivers (via direct control or some incentive scheme) and assign routes via a central authority that tries to minimize system-wide total travel time. This has been studied in the context of Stackelberg routing games [START_REF] Korilis | Achieving network optima using Stackelberg routing strategies[END_REF][START_REF] Roughgarden | Stackelberg scheduling strategies[END_REF][START_REF] Swamy | The effectiveness of Stackelberg strategies and tolls for network congestion games[END_REF], where the goal of the central controller is to assign routes to a fraction of the demand using a strategy that minimizes the system wide cost, while anticipating selfish behavior from the demand that is not being controlled. Most of this work has been in the area of communication networks and assumes a non-decreasing latency function and vertical queues. However, these assumptions are generally not satisfied in road traffic networks, with horizontal queuing, because of congestion propagation and more complex latency functions, due to the physics of flows and driver behavior [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF][START_REF] Daganzo | The cell transmission model, part II: Network traffic[END_REF][START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF]. Therefore, the literature on partial control in traffic assignment is sparse and and usually makes simplifying assumptions. For example, [START_REF] Aswani | Game-theoretic routing of GPS-assisted vehicles for energy efficiency[END_REF] use vertical queues and non-decreasing latency functions, while [START_REF] Krichene | Stackelberg Routing on Parallel Networks with Horizontal Queues[END_REF] consider simple parallel networks. [START_REF] Ziliaskopoulos | A linear programming model for the single destination system optimum dynamic traffic assignment problem[END_REF] formulated the single destination SO-DTA problem (with full control) as a Linear Program (LP) under a LP relaxation that approximates the non-linear dynamics of the system. However, the SO-DTA problem with partial control can not be formulated as a convex problem, even in the case of a single destination, without violating the first-in-first-out (FIFO) condition [START_REF] Carey | Nonconvexity of the dynamic traffic assignment problem[END_REF], due to the multiple commodities (selfish and cooperative demand). Furthermore, solving the SO-DTA problem with an LP relaxation of the dynamics can lead to the holding of vehicles on links when the model allows for a larger flow. It has been argued that this holding can be achieved in practice via variable speed limit (VSL) signs [START_REF] Ziliaskopoulos | A linear programming model for the single destination system optimum dynamic traffic assignment problem[END_REF] and makes sense when the goal of the problem is to also solve for the optimal VSL values [START_REF] Gomes | Optimal freeway ramp metering using the asymmetric cell transmission model[END_REF], but is impractical to implement in most cases. Thus, there is a need for a more general solution that does not require VSL.

A further complication of DTA in practical settings is the unavailability of origin-destination (OD) data for the entire demand. Most DTA solutions assume that this data is available, although it can be challenging to obtain in practice. Therefore, we formulate the partial control problem in a manner that requires full OD information only for the demand that can be controlled by the central authority, and junction split ratios for the remaining demand (which are much easier to obtain via inductive loop detectors for example).

We formulate the system optimal dynamic traffic assignment problem with partial control (SO-DTA-PC), using a traffic dynamics model similar to the Cell Transmission Model (CTM) [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF][START_REF] Daganzo | The cell transmission model, part II: Network traffic[END_REF], which is a Godunov discretization of the Lighthill-Williams-Richards (LWR) partial differential equation (PDE) [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF] with a triangular flux function. The CTM is a horizontal queuing model and uses a latency function that gives a constant delay when the traffic density is below a certain threshold and progressively increases as the density increases beyond this threshold, and is well accepted in the transportation community as a good first order approximation of road traffic dynamics.

We propose solving the SO-DTA-PC problem with the non-convex traffic dynamics from [START_REF] Delle Monache | A PDE-ODE Model for a Junction with Ramp Buffer[END_REF] and limited OD data with complete split ratios as a non-linear optimal control problem. This formulation generalizes to multiple sources and multiple destinations. The next challenge is in finding efficient descent methods for this non-convex optimal control problem. There is a vast literature on optimization techniques for non-convex control problems (see [START_REF] Bertsekas | Nonlinear programming[END_REF] and the references therein) that can be utilized to solve this problem. While gradient based methods do not provide any guarantees of converging to the optimal solution in non-convex optimization problems, they can still be used to find local minima. One of the main computational challenges in this approach is the efficient computation of the gradient, since this computation must be repeated a large number of times. We show that the structure of our dynamical system allows for very efficient computation of the gradient via the discrete adjoint method [START_REF] Bayen | Adjoint-based control of a new Eulerian network model of air traffic flow[END_REF]Giles andPierce, 1997, 2000;[START_REF] Jameson | Aerodynamic shape optimization techniques based on control theory[END_REF][START_REF] Raffard | An adjoint-based parameter identification algorithm applied to planar cell polarity signaling[END_REF]. If the state vector is n dimensional and the control vector is m dimensional, direct computation of the gradient takes O(n 2 m) time. The adjoint method generally reduces the complexity to O(n 2 + nm), but the structure of our system allows for further reduction of the complexity to O(n + m) by exploiting the sparse nature of the forward system.

It should be noted that this work currently does not model the response of the selfish users (demand). It is clear that a change in the network state will result in a response from the selfish users as in a Stackelberg game. Finding the optimal control for a Stackelberg game is NP-Hard in the size of the network for the class of increasing latency functions even in the static case [START_REF] Roughgarden | Stackelberg scheduling strategies[END_REF] and it is common to use approximate strategies [START_REF] Roughgarden | Stackelberg scheduling strategies[END_REF][START_REF] Swamy | The effectiveness of Stackelberg strategies and tolls for network congestion games[END_REF]. We wish to extend this work in the future to model the selfish response.

The contributions of this article are as follows: 1) Formulating of the system optimal dynamic traffic assignment with partial control problem as a multi-commodity finite horizon optimal control problem, 2) Defining the appropriate multi-commodity junction model for the network dynamics with reasonable data requirements, 3) Solving the gradient of the system with O(n + m) time complexity for an n dimensional state space and m dimensional control vector using the discrete adjoint method, and 4) Experimental results for illustrating the benefits and applications of the technique. The rest of the article is organized as follows. In Section 2, we present the traffic dynamics model with its assumptions, junction solver and boundary conditions. Following that, we introduce the notion of controllable and non-controllable flows in the network, in Section 3. In Section 4, we define the forward model of the system with explicit solutions to the junction problems. Section 5 formulates the optimization problem, derives the adjoint system and describes the gradient decent method. This is followed with numerical results in Section 6. We conclude with some final remarks and future research directions in Section 7.

Traffic model

The aggregate traffic dynamics are modeled using a macroscopic traffic flow model based on the Lighthill-Williams-Richards (LWR) PDE [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF]. We use a multicommodity variant with buffers of the PDE model developed in [START_REF] Garavello | Traffic flow on networks[END_REF]. This model imposes strong boundary conditions at the entrances to the network, so that no vehicles are dropped due to congestion propagating outside the bounds of the network, an important consideration in the optimal control setting. We then use a Godunov discretization [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF] of the network PDE model as explained in [START_REF] Reilly | An efficient method for coordinated ramp metering using the discrete adjoint method[END_REF] to obtain an equivalent discrete model.

Network definitions

The road network is divided into cells, indexed by i ∈ A. We add a buffer cell at the entrances of the network, to be able to impose the boundary demands in a strong sense. Each junction, indexed by z ∈ J , connects a set of incoming links J in z to a set of outgoing link J out z . The total flow in the network is decomposed into a set of |C| commodities that correspond to different types of flow.

Definition 1 (Supply and demand). The supply of a cell i at time step k, denoted σ i (k), is maximal flow that can enter the cell, while the demand δ i (k) is the maximal flow that can leave the cell. By assumption, buffers have no supply and the sinks have no demand.

Definition 2 (Density). The density on a link i at time step k, denoted by ρ i (k), is the total number of vehicles on the link during that time step divided by the length of the link L i . The vehicles in the link could be from any of the |C| commodities in the network.

Definition 3 (Single commodity density). The density induced by a single commodity c on a link i at time step k, denoted by ρ i,c (k), is the total number of vehicles of commodity c on the link during that time step divided by the length of the link L i , and satisfies

ρ i (k) = c∈C ρ i,c (k) (1)
Definition 4 (Initial conditions). The initial conditions of the network are the densities of each commodity at each link at time step k = 0 and are denoted ρ i,c (0).

Definition 5 (Inflow and outflow). The inflow (resp. outflow) from a cell i at time step k, denoted

f in i (k) (resp. f out i (k))
, is the total flow leaving (resp. entering) the cell at time step k. By assumption, buffers have no inflow and sinks have no outflow.

Definition 6 (Single commodity inflow and outflow). The inflow (resp. outflow) from a cell i at time step k corresponding to commodity c, denoted f in i,c (k) (resp. f out i,c (k)), is the total flow of commodity c leaving (resp. entering) the cell at time step k.

f in i (k) = c∈C f in i,c (k) (2) f out i (k) = c∈C f out i,c (k) (3)
Definition 7 (State evolution). The state of the network at time step k is given by the density ρ i,c (k) of each commodity c at each cell i. The density evolution is governed by the following dynamics, which simply correspond to mass conservation.

ρ i,c (k) = ρ i,c (k -1) + ∆t L i f in i,c (k -1) -f out i,c (k -1) ∀i ∈ A \ (B ∪ S), ∀k ∈ [[1, T f ]], ∀c ∈ C (4) ρ i,c (k) = ρ i,c (k -1) + ∆t L i • f in i,c (k -1) ∀i ∈ S, ∀k ∈ [[1, T f ]], ∀c ∈ C (5)
with initial condition

ρ i,c (0) = ρ 0 i,c ∀i ∈ A \ S, ∀c ∈ C (6) ρ i,c (0) = 0 ∀i ∈ S, ∀c ∈ C (7)
Assumption 1. The flux function defining the relationship between density and flow is given by the triangular fundamental diagram shown in figure 1. This is a first order approximation of the empirical relationship between flow and density [START_REF] Dervisoglu | Automatic calibration of the fundamental diagram and empirical observations on capacity[END_REF].

f ρ v -w ρ jam F Figure 1: Triangular fundamental diagram.
Assumption 2. First-in first-out (FIFO) property. We assume that no vehicles entering a link at time step t will overtake vehicles that have already entered the link at some time step t ′ < t.

Discretization requirements

The network model is composed of links (each link representing a physical road) and nodes (each node representing a junction between some incoming and outgoing links). As with any numerical scheme, the accuracy of a macroscopic traffic model increases with the granularity of the network discretization. In addition, to ensure the convergence of the solution of the discretized model to the solution of the continuous LWR equation as the time step ∆t and space discretization L goes to zero, the network must satisfy the Courant-Friedrichs-Lewy (CFL) conditions, which are standard requirements in numerical analysis [START_REF] Godunov | A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics[END_REF][START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. See Appendix B for a more detained explanation of the CFL requirements.

Controllable and non-controllable flow

There are two types of flows that are transported in the network. Controllable flows that have origin destination requirements, but can be routed along any path in the network, and non-controllable flows that have fixed paths. These flows are modeled by distributing the total flow of the network into multiple commodities, as explained below.

Assumption 3. Path decomposition of controllable flow. We assume that the controllable flows from each origin destination pair are restricted to a small pre-determined subset of paths in the network.

Definition 8 (Non-controllable commodity). There is a single non-controllable commodity c n that represents all non-controllable flow in the network.

The paths of the flow corresponding to the non-controllable commodity are defined via the junction split ratios.

Definition 9 (Split ratio). The split ratio of a commodity c at cell i and time step k among the outgoing cells j ∈ Γ(i), denoted β ij,c (k), is the fraction of the commodity c flow out of cell i that is entering cell j at time step k, such that j∈Γ(i) β ij,c (k) = 1.

Definition 10 (Controllable commodities). The controllable commodities c c ∈ CC correspond to the controllable flow. There is a unique controllable commodity that corresponds to each path that the controllable flow can be routed along in the network. A controllable commodity is then equivalent to a tuple (origin, destination, path).

Definition 11 (Conservation of flow). The path of a controllable commodity is defined via the junction split ratios for this commodity.

β ij,c (k) = 1 if the path of commodity c includes cell i and cell j 0 otherwise (8)
Definition 12 (Origin-destination demand). The number of controllable vehicles that seek to travel from origin o ∈ B to destination s ∈ S at time step k is given by

D (o,s) (k) • ∆t. It is an exogenous input.
Assumption 4. Data requirements. We assume that the origin destination information for all the controllable flow and the aggregate path information for all the non-controllable flow is known.

While at first glance this might seem like a lot of information to gather, it is in fact reasonable to assume in road traffic networks. We assume that the controllable flows are vehicles that are cooperating with the traffic coordination system that is trying to route vehicles efficiently and therefore will share their origin destination information. The aggregate paths of the non-controllable flows can be obtained from historical traffic patterns. The caveat is that empirical split ratios also include the contribution of the controllable flows and therefore must be pre-processed to remove this contribution.

Remark 1. The junction split ratios for the non-controllable commodities may be time-dependent, while the junction split ratios for the controllable commodities are not, since a controllable commodity corresponds to a single path.

Definition 13 (Equivalence between controllable commodity and path). Any controllable commodity c ∈ CC is a tuple (o, s, p) where o ∈ B, s ∈ S and p is a path (i.e. a sequence of cells). We define the function Ω as follows:

Ω : CC → S × B (9) c → (o, s).
Ω -1 ({o, s}) is then the set of commodities corresponding to the flows from source o to destination s.

Definition 14 (Controllable flow control). A control u is an allocation of the controllable demand over the set Ω -1 ({o, s}) for each time step. Formally u is defined as:

u : CC × [[0, T -1]] → [0, 1] (10) (c, k) → γ c (k) γ c (k)
(k) • γ c (k) • ∆t.
Definition 15 (Physically feasible control). A physically feasible control u verifies the mass conservation of the controllable demand allocation:

c∈Ω -1 ({o,s}) γ c (k) = 1 ∀k ∈ [[0, T ]], ∀(o, s) ∈ B × S (11)
We will denote with U be the set of all physically feasible controls.

3 Forward system

Junction model

The junction model defines the dynamics of the flow between neighboring cells. We require that is satisfies the following properties.

Requirement 1 (Multicommodity first-in first-out (FIFO) condition). For any outgoing link i, the distribution of its flow across the different commodities must be proportional to the ratio of vehicles of each commodity on the link. If ρ i (k) = 0 we must have:

f out i,c (k) = f out i (k) ρ i,c (k) ρ i (k) (12)
Requirement 2 (Consistency with split ratios). Let f ij,c (k) be the flow of commodity c from cell i to j at time step k. The outflow must be consistent with the split ratios.

f ij,c (k) = f out i,c (k) • β ij,c (k) (13)
Requirement 3 (Maximum flow constraint). The outflow cannot exceed the demand and the inflow cannot exceed the supply:

0 ≤ f in i (k) ≤ σ i (k) ∀k ∈ [[0, T ]] (14) 0 ≤ f out i (k) ≤ δ i (k) ∀k ∈ [[0, T ]] (15)
We wish to define a multi-commodity junction flow solver that assigns flows across the network in a manner that is consistent with the above requirements. The multi-commodity junction solver we consider is based on the source destination model (SDM) in [START_REF] Garavello | Traffic flow on networks[END_REF].

Definition 16 (Priority vector). In the case of a junction where there is more than one incoming cell and the aggregate demand of these cells is greater than the aggregate supply of the outgoing cells, the available supply needs to be distributed among the competing demands according to some priority vector as follows. The priority vector P j for cell j defines the allocation of its supply over the incoming cells i ∈ Γ -1 (j). The priority for a given incoming cell i is given by P ij such that i∈Γ -1 (j) P ij = 1.

Definition 17 (Aggregate split ratio). The aggregate split ratio β ij (k) over all commodities for a given path through a junction is defined as follows:

β ij (k) = c∈C ρ i,c (k) ρ i (k) β ij,c (k) = 1 ρ i (k) c∈C ρ i,c (k) β ij,c (k) (16)
The aggregate split ratio is only defined for positive aggregate densities, i.e. ρ i (k) > 0.

Diverge solver (1 × m)

We consider a diverging junction z with one incoming link i and m outgoing links. There are |C| commodities that flow through the junction each with their own time-varying split ratio β ij,c (k). If ρ i (k) = 0, then δ i (k) = 0 and f out i (k) is zero. We only consider the case of ρ i (k) = 0. Given the split ratios and densities of the cells at a junction, we wish to maximize the flow across the junction subject to the maximum flow constraints.

max f out i (k) (17) subject to 0 ≤ f in j (k) ≤ σ j (k) ∀j ∈ J out z 0 ≤ f out i (k) ≤ δ i (k)
We replace f in i (k) using the following relation:

f in j (k) = c∈C β ij,c (k) • f out j,c (k) [mass conservation] = f out i (k) c∈C ρ i,c (k) ρ i (k) • β ij,c (k) [by the FIFO constraint] = f out i (k) • β ij (k)
[by definition of the aggregate split ratios] ( 18)

Which gives us a trivial maximization problem that implies the following equality.

f out i (k) = min σ j (k) β ij (k) , ∀j ∈ J out z | β ij (k) > 0 , δ i (k) (19)
The total outflow f out i (k) for each incoming link i is then divided among the commodities according to the FIFO law:

f out i,c (k) = ρi,c(k) ρi(k) f out i (k).
The commodity flows are split among the outgoing links according to the split ratios constraints:

f in j,c (k) = β ij,c (k) f out i,c (k).
Existence and uniqueness of solution. A non-zero solution exists if none of the constraints of the optimization/feasibility problem imposes a zero flow. In other words, as long as the demand is non-zero and none of the outgoing links with positive demand (β ij (k) > 0) have non-zero supply, a non-zero solution exists. Since the solution to the maximum junction flow is given by equation ( 19) and the outflows are uniquely determined by the split ratios, the solution is unique.

Merge solver (n × 1)

We consider a merging junction z with n incoming links and one outgoing link j. A priority vector P j (s.t.

P ij = 1) prescribes the priorities at which the outgoing link accepts flows from the n incoming links when the junction is supply constrained.

If the problem is demand constrained (i.e.

i∈J in z δ i (k) < σ j (k))
, then the solution is given by:

f out i (k) = δ i (k) ∀i ∈ J in z (20)
Otherwise the problem is supply constrained and the solution to the junction problem is given by solving the following quadratic optimization problem that finds the flow-maximizing solution with the smallest violation of the priority vector, where the violation is measured using the L 2 distance:

min t,{f out i (k)} i∈J in z i∈J in z f out i (k) -t • p ij 2 (21) subject to i∈J in z f out i (k) = σ j (k) 0 ≤ f out i (k) ≤ δ i (k) ∀i ∈ J in z
The total outflow f out i (k) for each incoming link i is then divided among the commodities according to the FIFO law:

f out i,c (k) = ρi,c(k) ρi(k) f out i (k).
See figure 2a for a graphical illustration of the solution to a 2 × 1 junction.

The priorities are satisfied exactly when the intersection of the maximum flow isoline and the priority constraint is feasible. When this point is outside the feasible set, the flow-maximizing feasible point that is closest to the priority constraint (in euclidean distance) is chosen 1 .

Existence and uniqueness of the solution. In the demand-constrained case, existence and uniqueness are trivial. In the supplied constrained case, the solution is the feasible point that lies on the boundary of the feasible supply set (a segment) and minimizes the euclidean distance to the priority vector (a line), where the feasible set is given by the supply constraint (an n dimensional hyperplane:

i∈J in z f out i (k) ≤ σ j (k))
and demand constraints (an n-dimensional hyperrectangle:

f out i (k) ≤ δ i (k), ∀i ∈ J in z ).
A solution exists when the feasible set is non-empty, which is the case if the supply/demand constraints are greater than zero. This proves the existence of a solution in all non-degenerate (zero supply or demand) cases. The boundary of the supply constraint hyperplane intersects each coordinate axis at x i (k) = σ j (k) and can not be parallel to the priority constraint P , which is a line that goes through the origin. Therefore, since the solution must lie on a segment that is not parallel to the priority constraint line P , the solution that minimizes the distance to the P must be unique. This concludes the proof. See figure 2a for an illustration of the (2 × 1) case.

Merge-diverge solver (2 × m)

We consider a junction with 2 incoming links and m outgoing links 2 .

Assumption 5. The priority vectors P j for each outgoing link j are identical. This implies that the inflow priorities are allocated with respect to the total flow that enters the junction and that the priority does not depend on which outgoing link the vehicles will enter.

The priority vector P j prescribes the ratios at which the m outgoing links allocate their available supply to the 2 incoming links. It satisfies P i = P i1 = P i2 and i∈J in z P i = 1.

Let J in z and J out z be the sets of incoming and outgoing links at the junction. If the problem is demandconstrained (i.e.

i∈J in z β ij (k) δ i (k) ≤ σ j (k), ∀j ∈ J out z )
, then the solution is given by:

f out i (k) = δ i (k) ∀i ∈ J in z ( 22
)
1 The solution violates the priority rule only in the case where the demand for one or more of the incoming links is less than what its flow-maximizing allocation is, based on the priority vector. In other words, the priority rule is only violated when an incoming link doesn't have enough flow to satisfy its priority-based allocation. It is reasonable in the physical sense to maximize flow and only violate the priority when it is a lack of demand that causes the violation. The model is not denying any vehicles with priority the ability to pass through the junction. This is an important property to note, because it avoids having to solve a multi-objective optimization problem to come up with a physically meaningful set of flows through the junction.

2 We limit our analysis to merge-diverge junctions of no more than two incoming links because our model does not prescribe a unique solution when the number of incoming links is greater then two. Thus, our model can only be used with 1 × m, n × 1 and 2 × m junctions.

f out 2 (k) f out 1 (k) δ 2 (k) δ 1 (k) C : f out 1 (k) + f out 2 (k) = σ 1 (k) f out 1 (k) = P1 P2 f out 2 (k) (a) 2 × 1 junction f out 2 (k) f out 1 (k) δ 2 (k) δ 1 (k) C 1 C 2 f out 1 (k) = P1 P2 f out 2 (k) (b) 2 × 2 junction
Figure 2: An illustration of the solutions to merging junctions. The dashed lines denote the demand constraints imposed by the density on the incoming links. The solid lines (C) denote the supply constraints imposed by the density in the outgoing links. The solid lines going through the origin denote the merge priority vector. If the priority vector intersects any supply constraint inside the feasible demand set, the solution will be the feasible intersection point. If the intersection is outside the feasible demand set, then the solution will be the nearest feasible point where the supply and the demand constraints intersect (marked with a dot).

Otherwise, the flows through the junction are given by the following optimization problem.

min t,{f out i (k)} i∈J in z i∈J in z f out i (k) -t • P i 2 (23) subject to i∈J in z β ij (k) f out i (k) ≤ σ j (k) ∀j ∈ J out z max j   i∈J in z β ij (k) f out i (k) -σ j (k)   = 0 ∀j ∈ J out z f out i (k) ≤ δ i (k) ∀i ∈ J in z
The total outflow f out i (k) for each incoming link i is then divided among the commodities according to the FIFO law:

f out i,c (k) = ρi,c(k) ρi(k) f out i (k).
The commodity flows are split among the outgoing links according to the split ratio constraints: 2b for a graphical illustration of the solution to a 2 × 2 junction.

f in j,c (k) = i:(i,j)∈A β ij,c (k) f out i,c (k). See figure
Existence and uniqueness of solution. In the demand constrained case, existence and uniqueness are trivial. In the supplied constrained case, the solution is the feasible point (with respect to the supply and demand constraints) that lies on the boundary of the feasible supply set (a union of segments) and minimizes the euclidean distance to the priority vector (a line). A solution exists when the feasible set is non-empty, which is the case if the supply/demand constraints are greater than zero. This proves the existence of a solution in all non-degenerate (zero supply or demand) cases. The feasible supply set is the intersection of m two dimensional hyperplanes, which is a convex set. Therefore, the boundary of the feasible supply set (a union of segments) is also convex. Furthermore, the boundary of the feasible supply set intersects each coordinate axis at x i (k) = min j σj (k) βij (k) and therefore can not be parallel to the priority constraint P , which is a line that goes through the origin. Therefore, since the solution must lie on a convex union of segments and none of these segments is parallel to the priority constraint line P , the solution that minimizes the distance to P must be unique. This concludes the proof.

Boundary conditions

The boundary conditions at each source link of the network dictate the flows that enter the network. Each boundary condition is given as a flow rate at the boundary.

Definition 18 (Boundary demand). The number of vehicles of commodity c leaving cell i ∈ B at time step k is the boundary demand of commodity d i,c (k). Let c n be the commodity corresponding to non-controllable flow. The non-zero terms are defined as:

d i,c (k) = ∆t • D i,c (k) ∀i ∈ B, c = c n (Non-controllable demand) d i,c (k) = ∆t • D Ω(c) (k) • γ c (k) ∀i ∈ B, ∀c ∈ CC (Controllable demand)
Since the inflow to the network is limited by the maximum flow capacity and density of the immediate downstream link, all of the demand at a given time step might not make it into the network. A source buffer is used to accumulate the flow that cannot enter the network to guarantee the conservation of boundary flows. In the single commodity case, this model is sufficient. However, in the multi commodity case, we also need to make sure that the flow through the boundary respects the multicommodity FIFO condition given in requirement 1.

As stated in Definition 18, d i,c (k) is the boundary demand per commodity on cell i at time step k. The FIFO condition dictates that the vehicles entering the boundary buffer at time k should enter link i at the ratio di,c(k) c∈C di,c(k) for each commodity c. The simplest solution is to have a single buffer l at the boundary, as in the single commodity case, and keep track of how many vehicles of each commodity are at the buffer. The flow into the boundary cell will be as follows:

f in i,c (k) = l i,c (k) l i (k) f in i (k)
This equality satisfies the FIFO condition assuming that the vehicles in the buffer are uniformly distributed. However, in reality the buffer can accumulate vehicles arriving at the boundary at different time steps with different commodity ratios di,c(k) di(k) . Thus, this model can violate the FIFO property across multiple time steps if some vehicles cannot leave the buffer in one time step. See Appendix C for a more detailed discussion on the FIFO condition and a multi buffer model that satisfies the FIFO condition.

As the length of a buffer can be seen as the density of a cell of length 1, we use the same notation ρ i,c (k) for a buffer i. The speed of this cell is then v i = 1 ∆t because of requirements 4 and 6.

System dynamics

For a given control u, we can determine the evolution of the network using the following equations that prescribe the system dynamics. Let x(u) give the state of the network under these dynamics subject to the control u.

The system of equations governing the evolution of the network (implicit definition of x) are written formally in the form H(x, u) = 0, thus x is an implicit function of u. The discretized system dynamics can be described using six types of constraints, given by H h = 0, h ∈ {1, . . . , 6}, listed below. These six constraints have different individual instantiations depending on the specific setting such as link type or junction type. The explicit formulation is given below.

Mass conservation

H 1 k,i,c : ρi,c (k) = ρi,c (k -1) + ∆t Li f in i,c (k -1) -f out i,c (k -1) ∀i ∈ A \ (B ∪ S), ∀k ∈ [[1, T f ]], ∀c ∈ C (H1a) H 1 k,i,c : ρi,c (k) = ρi,c (k -1) + ∆t Li • f in i,c (k -1) ∀i ∈ S, ∀k ∈ [[1, T f ]], ∀c ∈ C (H1b)
with initial conditions

H 1 0,i,c : ρi,c (0) = ρ 0 i,c ∀i ∈ A \ S, ∀c ∈ C (I1a) H 1 0,i,c : ρi,c (0) = 0 ∀i ∈ S, ∀c ∈ C (I1b) Boundary conditions 3 H 1 k,i,c : ρi,c (k) = ρi,c (k -1) + ∆t Li D Ω(c) (k) • γc(k) -f out i,c (k -1) ∀i ∈ B, ∀k ∈ [[1, T f ]], ∀c ∈ CC (H1c) H 1 k,i,c : ρi,c (k) = ρi,c (k -1) + ∆t Li Di,c (k) -f out i,c (k -1) ∀i ∈ B, ∀k ∈ [[1, T f ]], c = cn (H1d)
with initial conditions

H 1 0,i,c : ρi,c (0) = ρ 0 i,c + ∆t Li • D Ω(c) (0) • γc(k) ∀i ∈ B, ∀c ∈ CC (I2a) H 1 0,i,c : ρi,c (0) = ρ 0 i,c + ∆t Li • Di,c (0) ∀i ∈ B, c = cn (I2b) Flow propagation Recall that ρi (k) = C c=1 ρi,c (k) is the total density of cell i at time step k. H 2 k,i : δi (k) = min (Fi, viρi (k)) ∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] (H2a)
H 2 k,i : δi (k) = min Fi, ρi (k) Li ∆t ∀i ∈ B, ∀k ∈ [[0, T f ]] (H2b) H 3 k,i : σi (k) = min Fi, wi ρ jam i -ρi (k) ∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] (H3a) H 3 k,i : σi (k) = Fi ∀i ∈ S, ∀k ∈ [[0, T f ]] (H3b)

Junction solution

The derivation of the explicit solutions to the 1 × 2, 2 × 1 and 2 × 2 junctions listed below are given in appendix D. The general solutions for the junction model are given in section 3.1.

To simplify the notation, we use the following shorthand:

• We drop the time index k

• We abbreviate i 1 = 1 and i 2 = 2
• We use the following notation i 1 = i 2 and i 2 = i 1

When ρ i = 0, there are no vehicles in the incoming link i and the outflow of this link is zero for all commodities. The following equations only apply for ρ i = 0. Aggregate split ratios:

H 4 k,i,j,z : βij = 1 ρi c∈C ρi,cβij,c ∀z ∈ J , ∀i ∈ J in z , ∀j ∈ J out z , ∀k ∈ [[0, T f ]] (H4)
Flow out of incoming links by commodity:

H 5 k,i,c : f out i,c = ρi,c ρi min σj βij , ∀j ∈ J out z | βij > 0 , δi ∀z ∈ J1×n, ∀i ∈ J in z (H5a) H 5 k,i,c : f out i,c = ρi,c ρi      δi if Pi (min (δ1 + δ2, σ1) -δi) > δiPi min (δ1 + δ2, σ1) -δi if Pi (min (δ1 + δ2, σ1) -δi) > δiPi Pi min (δ1 + δ2, σ1) otherwise ∀z ∈ J2×1, ∀i ∈ J in z (H5b) H 5 k,i,c : f out 1 =            δ1 if P 1 P 2 > δ 1 min δ 2 , σ 1 -β 11 δ 1 β 21 , σ 2 -β 12 δ 1 β 22 min δ1, σ 1 -β 21 δ 2 β 11 , σ 2 -β 22 δ 2 β 12 if P 1 P 2 < min δ 1 , σ 1 -β 21 δ 2 β 11 , σ 2 -β 22 δ 2 β 12 δ 2 min P 1 σ 1 P 1 β 11 +P 2 β 21 , P 1 σ 2 P 1 β 12 +P 2 β 22 otherwise ∀i ∈ J in 2×2 (27a) f out 2 is obtained by symmtery. ( 27b
)
f out i,c = ρ i,k ρi f out i ∀i ∈ {1, 2}, ∀c ∈ C (H5c)
Flow into outgoing links by commodity:

H 6 k,i,c : f in i,c = x∈J in z βxi,cf out x,c ∀z ∈ J , ∀i ∈ J out z , ∀c ∈ C (H6)
As shown in figure 3, the dynamics equations have a topological ordering that allows for the following efficient forward simulation algorithm 4 . Time 0

Time T f = T -1 Time k ••• (ρ 1,c (0)) ••• (ρ N,c (0)) δ 1 (0) ••• δ N (0) σ 1 (0) ••• σ N (0) (β 1j (0)) ••• (β Nj (0)) f out 1,c (0) ••• f out N,c (0)
f in 1,c (0) ••• f in N,c (0)
(ρ 1,c (T f )) ••• (ρ N,c (T f )) δ 1 (T f ) ••• δ N (T f ) σ 1 (T f ) ••• σ N (T f ) (β 1j (T f )) ••• (β Nj (T f )) f out 1,c (T f ) ••• f out N,c (T f ) f in 1,c (T f ) ••• f in N,c (T f )
Figure 3: Dependency diagram of the variables in the system.

4 For any given time step k, each of the three internal loops are trivially parallelizable problems.

4 Adjoint based optimization

Problem formulation

For a given control u, we can determine the evolution of the network using the equations for the system dynamics. Let x(u) be the state of the network under these dynamics subject to the control u.

The total travel-time J(x(u)) is defined as:

J = T -1 k=0 i∈A\S ρ i (k) • L i (28)
The system optimal dynamic traffic assignment with partial compliance (SO-DTA-PC) is a physically acceptable (see Definition 14) division of the compliant agents among the different commodities that minimizes the total travel-time (including the travel-time of the non-compliant commodities). The solution is obtained by solving the following optimization problem.

min u∈U J(x(u))
subject to system dynamics control constraints

where the system dynamics are given in Section 3.3 and the control constraints are the following.

γ c (k) ≥ 0 ∀c ∈ CC, k ∈ [[0, T f ]] c∈Ω -1 {(o,s)} γ c (k) = 1 ∀k ∈ [[0, T f ]]
Note that this is a non-convex optimization problem that might contain multiple local minima. Therefore, gradient methods will not guarantee global optimality. However, descent algorithms can still be used to obtain locally optimal solutions and can be improved by using multiple starting points [START_REF] Boese | A new adaptive multi-start technique for combinatorial global optimizations[END_REF][START_REF] Marti | Multi-start methods[END_REF]. Furthermore, non-convex optimization techniques such as subgradient and interior point methods [START_REF] Giles | Adjoint equations in CFD : duality , boundary conditions and solution behaviour[END_REF][START_REF] Wachter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF] require the gradient of the system. We use the discrete adjoint method, which will be explained in the next section, to efficiently solve for the gradient of the system. The control constraints can be satisfied either using a projected gradient descent or a barrier function. In our implementation, we use the projected gradient descent approach.

Overview of the adjoint method

We consider the following general optimization problem:

min u∈U J(x, u) subject to H(x, u) = 0 (29)
where x ∈ X denotes the state variables and u ∈ U denotes the control variables.

The adjoint method [START_REF] Duffy | An Introduction to Gradient Computation by the Discrete Adjoint Method[END_REF]) is a technique to compute the gradient ∇ u J(x, u) = dJ du of the objective function without fully computing ∇ u x = dx du . The gradient is then used to perform a gradient descent. We suppose that for any control u, ∂H ∂x (x, u) is not singular. Under equality constraints H(x, u) = 0, the Lagrangian

L(x, u, λ) = J(x, u) + λ T H(x, u) (30)
coincides with the objective function for any feasible point (x(u), u). The problem is then equivalent to computing the gradient of the Lagrangian:

∇ u L(x, u, λ) = ∂J ∂u + ∂J ∂x dx du + λ T ∂H ∂u + ∂H ∂x dx du = ∂J ∂u + λ T ∂H ∂u + ∂J ∂x + λ T ∂H ∂x dx du (31)
In particular, if λ satisfies the adjoint equation:

∂J ∂x + λ T ∂H ∂x = 0 (32)
then the gradient is,

∇ u L(x, u) = ∂J ∂u + λ T ∂H ∂u (33)
The solution for λ exists and is unique if ∂H ∂x is not singular, which is the case in our forward system, as explained in the following section.

Applying the adjoint method

To be able to use the adjoint method to compute the gradient, the partial derivatives of the forward system with respect to the state variables ∂H ∂x must not be singular. We can rewrite our system of equations in the form H(x, u) = 0 and verify this condition trivially.

All the diagonal terms of ∂H ∂x are non zero (since equal to 1 or -1 depending on the way we rewrite H v ). As seen in the dependency chain shown in Figure 3, the non zero derivative terms of H v depend only on variables that have a smaller index in x. This means that ∂H ∂x is lower triangular with no zero terms on the diagonal and is thus non singular. Therefore, we can apply the adjoint method to compute the gradient of this system.

Reduced state space. The forward system dynamics that were described in section 3.3 had a large number of state variables. However, the only required state variables of the system are the partial densities ρ i,c (k). All the others variables were introduced to make the forward system easier to understand. We will now drop most of these auxiliary variables to simplify the computation of the adjoint system. We only use ρ i,c (k), f out i,c (k) and f in i,c (k) to describe the system and replace the other variables by their expressions as a function of the three state variables that we retain. We define:

x =     (ρ i,c (k)) c∈C i∈A f out i,c (k) c∈C i∈A f in i,c (k) c∈C i∈A     k∈[[0,T ]] H =         H 1 k,i,c c∈C i∈A H 5 k,i,c c∈C i∈A H 6 k,i,c c∈C i∈A         k∈[[0,T ]]
Computational complexity. Let n be the dimension of the state vector x ∈ R n , m be the dimension of the control vector is u ∈ R m and N c = |C| be the total number of commodities. From the above definition of the state vector, we can see that n = |A| • T • N c . The dimension of H is also n as defined above.

Direct computation of the gradient ∇ u J(x, u) takes O(n 2 m) time.

∇ u J(x, u) = ∂J ∂x • ∇ u x + ∂J ∂u (34)
Computing ∇ u J requires solving the system H(x, u) = 0 ⇒ ∂H ∂x dx du + ∂H ∂u = 0, which is equivalent to solving m different n × n linear systems and takes O(n 2 m) time. The final step of multiplying ∂J ∂x dx du and adding ∂J ∂u takes O(nm) time, but is dominated by the time to compute dx du . The discrete adjoint method reduces this complexity to O(n 2 + nm) by solving for λ in the adjoint system. Computing the adjoint variables λ T ∈ R n using equation (32) only takes O(n 2 ) time since it only requires solving one n × n linear system. Multiplying λ T ∂H ∂u and adding ∂J ∂u takes O(nm) time, so the total computation time is O(n 2 + nm).

The structure of our system allows for further reduction of the complexity to O(n + m|C|). As shown in section 4.3, ∂H ∂x is a lower triangular matrix and therefore we can compute the solution to equation (32) using backwards substitution. We will exploit the fact that the matrix ∂H ∂x is extremely sparse. The maximum row cardinality is four because the forward system does not contain any constraints with more than four variables. Therefore, equation ( 32) can be solved in O(n) time. If the maximum in-degree of the network is d in , the maximum column cardinality is 2 + |C|(1 + d in ), as will be clear in the next section from equation (38). Assuming that d in is a small constant, the multiplication step in equation ( 33) takes O(m|C|) time. This leads to a total computation time of O(n + m|C|).

Adjoint equations

The adjoint equations are given by the system:

∂J ∂x + λ T ∂H ∂x = 0 (35) ⇒ ∂J ∂x + x ′ ∈x λ ′ x ∂H x ′ ∂x = 0 ( 36
)
where x is the state vector and H ′ x is the forward system equation corresponding to the variable x ′ ∈ x. To write the adjoint system equation corresponding to x ′ , we have to look at all the forward system equations where x ′ appears and consider all the non-null ∂H x ′ ∂x terms. In particular we write the equations such that

∂H x ′ ∂x ′ = -1.
Note that this can be done because the Godunov scheme provides an explicit expression for the forward system constraints.

Computing ∂J ∂x ∂J ∂ρ i,c (k) = L i ∀c ∈ C, ∀i ∈ A \ S, ∀k ∈ [[0, T ]] 0 otherwise (37) Computing λ T ∂H ∂x ∂H ∂ρ i,c (k) : x ′ ∈x λ x ′ ∂H x ′ ∂ρ i,c (k) = λ ρi,c(k) ∂H ρi,c(k) ∂ρ i,c (k) + λ ρi,c(k+1) ∂H ρi,c(k+1) ∂ρ i,c (k) + c ′ ∈CC   λ f out i,c ′ (k) ∂H f out i,c ′ (k) ∂ρ i,c (k) + λ f out x,c ′ (k) x:(x,i)∈A ∂H f out x,c ′ (k) ∂ρ i,c (k)   (38) ∂H ∂f out i,c (k)
:

x ′ ∈x λ x ′ ∂H x ′ ∂f out i,c (k) = λ ρi,c(k+1) ∂H ρi,c(k+1) ∂f out i,c (k) + λ f out i,c (k) ∂H f out i,c (k) ∂f out i,c (k) + j:j∈(i,j) λ f in j,c (k) ∂H f in j,c (k) ∂f out i,c (k) (39) ∂H ∂f in i,c (k)
:

x ′ ∈x λ x ′ ∂H x ′ ∂f in i,c (k) = λ ρi,c(k+1) ∂H ρi,c(k+1) ∂f in i,c (k) + λ f in i,c (k) ∂H f in i,c (k) ∂f in i,c (k) (40)
The next section shows how to compute all the individual partial derivatives that are required in the above equations. Once they are computed, we can simply plug them into the above equations and solve the system via backwards substitution since ∂H ∂x is lower triangular.

Partial derivatives

Computing the gradient of the system via the adjoint method requires computing the partial derivatives ∂J ∂u , ∂J ∂x , ∂H ∂u and ∂H ∂x . The first three of these can be computed trivially.

• Partial derivatives of the cost function with respect to the control ( ∂J ∂u ) from equation ( 28)

∂J ∂γ c (k) = 0 (41)
• Partial derivatives of the cost function with respect to the state variables ( ∂J ∂x ) from equation ( 28)

∂J ∂ρ i,c (k) = L i if c ∈ C, i ∈ A \ S, k ∈ [[0, T ]] 0 otherwise (42)
• Partial derivatives of the constraints with respect to the state variables ( ∂H ∂u ) from equation ( H1c)

∂ρ i,c (k) ∂γ c (k) = ∆t Li • D Ω(c) (k) if c ∈ CC, i ∈ B, k ∈ [[0, T ]] 0 otherwise (43)
Finally, we need to compute the partial derivatives of the constraints with respect to the state variables ( ∂H ∂x ), which is much more algebraicly involved and therefore described in appendix E.

Numerical Results

To illustrate the effectiveness of our framework for computing the system optimal dynamic flow allocation with partial control, we have implemented the algorithm and tested it on both synthetic and practical traffic rerouting scenarios using experimental field data. Our implementation uses the discrete adjoint method to compute the gradient and uses a log barrier function with a projection step to keep the solution in the physically feasible control set as explained in section 4.1. We use the Rprop [START_REF] Riedmiller | Rprop-a fast adaptive learning algorithm[END_REF] algorithm as our gradient descent technique. All the experiments were run on a 1.8 GHz Intel Core i5 dual-core processor with 8GB of RAM. The performance cost C of each scenario is measured using the total travel time of all the vehicles passing through the network.

C = T -1 k=0 i∈A\S ρ i (k) • L i • ∆t
We present numerical results for two network scenarios. 1. A network adapted from the synthetic network used in [START_REF] Ziliaskopoulos | A linear programming model for the single destination system optimum dynamic traffic assignment problem[END_REF]. The network is illustrated in figure 4.

2. A subsection of Interstate 210 with a parallel arterial route, as depicted in figure 6.

Synthetic network

The synthetic network is a simple 10 cell network adapted from the example used in [START_REF] Ziliaskopoulos | A linear programming model for the single destination system optimum dynamic traffic assignment problem[END_REF] 5 . It contains three paths over which vehicles can be routed. The time discretization is set to one time unit and the length of each cell is also normalized to one unit. Therefore, the total capacity of each cell in terms of the number of vehicles N is equal to the jam density ρ jam . Each cell in figure 4 is annotated with its cell capacity N , while the edge weights in the network prescribe the max flow F between the cells. The free flow speed v of each cell is also normalized to one and the congestion speed w is equal to the free flow speed. The demand at the origin at the first three times steps is respectively 8, 16 and 8 vehicles. The network is simulated for 10 time steps, which gives enough time for all the entering flow to exit the network. First we use the discrete adjoint optimization framework to compute the system optimal flow allocation for the network assuming that all of the flow is compliant. Table 1(a) shows the optimal route allocation for the origin demands at each time step with non-zero demand. The total travel time cost (C) with the optimal flow allocation is 178 time units. The solution converges to within 0.5% of the optimal solution in three iterations. Capacity reduction due to incident. We now consider the case where the capacity between cell 3 and cell 4 is temporarily reduced due to some incident. The corresponding capacities for link (3, 4) are given in table 2.

Time step 1 2 3 4 5 6 7 8 9 10 F (3,4) 6 6 0 0 3 3 6 6 6 6

Table 2: Capacity reduction due to incident

If the vehicles continue to be routed using the previous path allocation, the total travel cost (C) will now be 244 time units. The total cost increases by 37% because a large percentage of vehicles are routed along the path that is temporarily closed and then subjected to a reduced capacity. If we recompute the system optimal flow allocation, the total cost decreases to 211 time units and the corresponding flow allocation is given in table 1(b). This solution is actually a local minimum in the system due to the FIFO condition for vehicles departing cell 3. Whenever there is some non-zero flow for route 1 when the capacity is zero, the flow of vehicles that take route 2 is also restricted to zero. This causes a non-convexity that results in a discontinuity of the gradient at the point where the flow of vehicles on route 1 is zero. Gradient-based methods are not well suited to deal with such conditions because the information obtained from the gradient only provides local information. The global optimal solution occurs with the flow allocation given in table 1(c) and results in a total travel time cost of 207. As mentioned in Section 4.1, the effect of local minima can in general be mitigated using multi-start strategies [START_REF] Boese | A new adaptive multi-start technique for combinatorial global optimizations[END_REF][START_REF] Marti | Multi-start methods[END_REF] and the efficient gradient computation obtained via the adjoint method can be combined with non-convex optimization techniques such as interior point methods [START_REF] Wachter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF]. Partial control. In many situations, it might not be possible to reroute all the vehicles in the system. Therefore, we also analyze the behavior of the system when only some fraction of the vehicles are rerouted. Figure 5 shows how the total travel time changes with the percentage of vehicles that can be rerouted. In this example, we see that the system optimal (local minimum) can be achieved by controlling only 60% of the vehicles.

Interstate 210 network

The experimental analysis was conducted on a 8 mile corridor of Interstate 210 in Arcadia, California with a parallel arterial route, as illustrated in figure 6. The network has 24 cells corresponding to satisfying the CFL condition for a time step of 30 seconds. The physical properties of the network such as the capacity were obtained using the Scenario Editor software developed as part of the Connected Corridors project, a collaboration between the University of California Berkeley and California Partners for Advanced Transportation Technology (PATH). Calibrated fundamental diagram parameters, split ratios, and boundary data were also obtained from other parallel research efforts at Connected Corridors. The data used for calibrating these parameters was obtained from the Freeway Performance Measurement System (PeMS) [START_REF] Chen | Freeway performance measurement system: mining loop detector data[END_REF]. We consider a prototypical one hour time horizon during the morning commute. The density profile of the freeway under the calibrated parameters and estimated boundary flows is shown in figure 7(a).

Capacity drop during morning commute. We analyze the behavior of the freeway corridor in the event of a capacity drop caused by some incident. We assume that the capacity drop occurs at the fifth freeway road segment 10 minutes into the simulation and that it lasts for 20 minutes, as illustrated in figure 7 (b). The freeway capacity at segment five will be assumed to be reduced by half during this period, corresponding to a closure of two lanes (out of four) at the location of the incident. Figure 7 shows the density profile corresponding to; (a) normal operation with capacity drop, (b) a capacity drop due to a two lane closure during the incident with no traffic diversion, (c) the same capacity drop with traffic being diverted to the parallel arterial, and (d) the change in the density profile due to the traffic diversion. As the figure shows, rerouting the excess flow to the parallel arterial eliminates the bottleneck during the incident and improves the throughput of the freeway corridor. In this example, the parallel arterial is assumed to prioritize vehicles being routed from the freeway and the full arterial capacity is used for this purpose. However, in certain situations, municipalities may want to allocate some capacity of the parallel arterial for local traffic. In this case, the optimizer can be limited to only use a certain fraction of the capacity of the parallel arterial. Figure 8 shows the density evolution when the arterial capacity allocated for rerouting freeway traffic is limited to 40% and 50% in comparison to full arterial utilization. The arterial capacity allocation can be controlled via the traffic signal controls along the arterial. Adjoint method vs finite differences. To demonstrate the increased efficiency of computing the solution via the discrete adjoint method, we also implemented the gradient computation using a simple finite differences method6 by perturbing each variable and measuring the response of the system. This approach, which approximates the gradient at a given point, has a runtime complexity of O(nm) where n is the dimension of the state vector and m is the dimension of the control vector. Recall that the size of the state vector is n = |A| • T • |C| and that the size of the control vector is m = T • |C|, where |C| is the total number of commodities (feasible paths) in the problem. Therefore, O(nm) = O(|A| • T 2 • |C| 2 ) and the finite difference method has a computation time that is quadratic in the number of time steps. In comparison, the adjoint method has a time complexity of O(n

+ m|C|) = O(|A| • T • |C| + T • |C| 2 )
, which is linear in the number of time steps. The complexity of both methods is quadratic in the total number of feasible paths, but this is assumed to be a small number in practical routing problems, since vehicles that travel between a fix origin-destination pair will only typically have a small number of reasonable paths. Figure 9 shows the time taken for one gradient computation as a function of the network number of time steps in the problem for the I-210 network. The simulations are run by changing the time discretization of the problem to control the total number of time steps. Reducing the time discretization also increases the number of cells due to the maximum cell length imposed by the CFL condition. The results show that the finite differences approach quickly becomes computationally intractable as the number of time steps in the problem increases and highlights the value of the discrete adjoint method for solving large problems in a tractable manner.

Conclusion

This article has presented a model and optimization framework for solving the System Optimal Dynamic Traffic Assignment problem with Partial Control (SO-DTA-PC) for general networks with horizontal queuing dynamics. The model only requires full origin-destination (OD) information for the fraction of the agents that are controllable, with aggregate split ratios being sufficient for the non-controllable (selfish) agents. We have utilized a Godunov discretization of the Lighthill-Williams-Richards (LWR) partial differential equation with a triangular flux function and a corresponding multi-commodity junction solver for the dynamics of the system, but the proposed framework can also be used with other dynamics models. Furthermore, we show that the sparsity pattern of the forward system allows us to compute the gradient of the system with linear computational complexity and memory using the discrete adjoint method. Finally, we apply this framework to find the optimal vehicles rerouting strategy in response to a capacity loss in the network, and show the congestions reductions that can be achieved. Numerical results are presented for a test network and Interstate 210 in Southern California The cells (including buffers and sinks) B

The buffer cells S

The sink cells OD

The set of origin-destination (OD) pairs Γ -1 (i) The predecessors of cell i Γ(i)

The successors of cell i C

The commodities CC The controllable commodities Inputs

ρ i,c (
f in i (k) Total flow into cell i, time step k f out i (k) Total flow out of cell i, time step k f in i,c (k) Flow of commodity c into cell i, time step k f out i,c (k) Flow of commodity c out of cell i, time step k ρ i (k) Density on cell i, time step k ρ i,c (k) Density contribution of commodity c on cell i, time step k σ i (k) Supply on cell i, time step k δ i (k) Demand on cell i, time step k d i (k) Boundary demand on cell i, time step k γ c (k)
Demand allocation for commodity c, time step k 7 Time from t = 0 and the final time step T = T -1.

B CFL conditions

Requirement 4 (CFL condition 1). ∀i ∈ A \ S, v i ≤ Li ∆t For numerical stability, the vehicles in a given cell should only be able to travel forward at most one cell in a single time step. Requirement 4 ensures that this condition is satisfied by imposing an upper bound on the velocity.

Requirement 5 (CFL condition 2). ∀i ∈ A \ B, w i ≤ Li ∆t Each cell can not have a density greater than ρ jam . Requirement 5 ensures that this condition is satisfied. While the requirement that v i ≤ Li ∆t comes from the positive density, this one comes from the fact that the density has to be smaller than ρ jam i . Indeed, in the case of a cell with no outflow (because of an extreme congestion of the next cell), the inflow can be limited by the supply. If there is enough demand we have at the next time the density which has to be smaller than ρ jam i :

ρ i (k) + ∆t L i w i (ρ jam i -ρ i (k)) < ρ jam i ρ i (k) 1 -w i ∆t L i < ρ jam i 1 -w i ∆t L i And because 0 ≤ ρ i (k) ≤ ρ jam i we also have 1 -w i ∆t
Li ≥ 0 which is the requirement.

Requirement 6 (Finite time density discharge). ∀i ∈ A, v i ≥ Li ∆t This condition guarantees that the density of a given cell discharges in a finite amount of time when there is no incoming flow. Requirement 6 ensures that this condition is satisfied. If v i < Li ∆t , we can have exponential decrease of the density in some cells when they should be emptied in only 1 steps. Taking the case of a cell without inflow, we have ρ

i (k + 1) = ρ i (k) -∆t Li ρ i (k) v i which gives ρ i (k + t) = ρ i (k) (1 -∆tvi Li ) t
. This is not an acceptable physical solution and thus should be excluded.

Remark 2. Satisfying requirements 4 and 6 implies imposing v i = Li ∆t . Given that the velocity v is an exogenous parameter, this imposes a strict requirement on the space discretization of the road segments. However, any given road segment might not be strictly into cells of exact length v i • ∆t, and in most cases a cell of length L ∈ (0, v i ∆t) will remain. There are multiple solutions to this issue:

• Approximate the length of each road segment to be a multiple of v i • ∆t. The relative rounding error decreases as the road gets longer and the descretization ∆t gets smaller.

• Change the dynamics of the last cell in each road segment to have a special case that allows for vehicles to be fully discharged when the supply of the downstream link allows it. This makes the dynamics equations and optimization problem more complicated.

• Accept this model limitation and have a small amount of density stuck in the network. This is not so bad in practice, since the number of vehicles stuck in a link decreases exponentially with time.

C FIFO condition at the origin

It is important to note that the FIFO condition is violated in its strict sense even within the network. The flow propagation model assumes that all the flow within a cell is uniformly distributed according to the individual commodity ratios regardless of when the vehicles arrived at the cell.

Example 1. Consider the following simple example. There are two commodities a, b in cell i with 10 vehicles of each commodity at time k. At time k + 1, 10 vehicles exit the cell (5 of a and 5 of b by the FIFO rule) and 10 new vehicles (3 of a and 7 of b) enter the cell. The new ratio of vehicles at i is 8 a to 12 b. At time k + 2, once again 10 vehicles exit the network. According to the cell level FIFO rule, the 10 vehicles will consist of 4 a's and 6 b's. However, the first 10 cars of those currently in cell i came at the ratio of 1:1 and truly satisfying the FIFO rule would require the 10 exiting vehicles to consist of 5 a's and 5 b's.

Remark 3 (Cell-level multicommodity FIFO condition). The strict multicommodity FIFO condition is not satisfied in most traffic flow models, but this is considered to be acceptable by the traffic modeling community.

The multicommodity FIFO condition is therefore in practice at best limited to the cell level.

While this argument is satisfactory at the interior of the network, due to the bounded number of vehicles in a given cell (due to the jam density), the FIFO violation can be significant at the boundaries of the network.

The number of vehicles in a source buffer at the boundary of the network can be arbitrarily large, and thus, there is no bound on how badly the multi-commodity FIFO constraint can be violated.

A simple extension that limits the mixing of vehicles entering at different time steps is to have a series of source buffers at the boundary, where the multicommodity FIFO constraint is enforced when vehicles move between the buffers. In this model, the commodity ratios are maintained separately for each buffer. Any vehicles that enter the network at a given time step are added to the last active (non-empty buffer) buffer. This restricts the violation of the FIFO condition across multiple steps to the capacity of a single buffer.

The capacity of each buffer is chosen such that the buffer can satisfy the maximum supply of the first cell of the network boundary that the buffer serves. This prevents artificial delays at the source. The buffer capacity is set to ∆t • F b where b is the boundary cell.

The limitation of this model is that we need to maintain l max ∆t• F b buffers per source, where l max is the maximum queue length at the boundary and ∆t • F b is the capacity of each buffer. Let l b i,c (k) be the number of vehicles of commodity c in the b th buffer for source node i at time step k and l b i (k) its sum over all the commodities.

The buffers are updated as follows:

• First, given f in i (k), move flow out of the initial buffer

f in i,c (k) = l 1 i,c (k) l 1 i (k) f in i (k) ∀c ∈ C (44) l 1 i,c (k + 1) = l 1 i,c (k) -f in i,c (k) ∀c ∈ C (45)
• Let B be the number of buffers in use. Iterate through the buffers and push flow upstream using the following algorithm. The only demand that is exposed to the system dynamics and the optimization problem is the demand that is captured in the first buffer.

D Explicit solutions for 1 × 2, 2 × 1 and 2 × 2 junctions

In this section we present the explicit solutions used for the dynamics of the optimization problem and in the software implementation. We restrict our attention to 1 × 2, 2 × 1, and 2 × 2 junctions, since most road networks can be modeled using only these junction types.

(1 × 2) diverge junction

The explicit solution was already derived in section 3.1.1 for the general (1 × n) case.

(2 × 1) merge Let i 1 be one incoming link, i 2 be the other incoming link and j be the outgoing link. As shown in Figure 10a, we distinguish 3 cases based on where the priority vector intersects the demand constraints. The cases are

P1 P2 > P1 P2 max , P1 P2 min ≤ P1 P2 ≤ P1 P2 max and P1 P2 < P1 P2 min
, where

P1 P2 max = δi 1 f in j -δi 1 and P1 P2 min = f in j -δi 2 δi 2 .
In call cases, the flow into the outgoing link is the minimum of the total supply and demand values.

f in j (k) = min δ i1 (k) + δ i2 (k) , σ j (k) (46)
The flow out of incoming links depends on the priority vector and is given by the solution to the optimization problem in equation ( 21).

f out i1 (k) =      δ i1 (k) if P i1 f in j (k) -δ i1 (k) > δ i1 (k) P i2 f in j (k) -δ i2 (k) if P i2 f in j (k) -δ i2 (k) > δ i2 (k) P i1 P i1 f in j (k) otherwise f out i2 (k) = f in j (k) -f out i1 (k) (47)
The flow out by commodity can then be computed as follows. (2 × 2) merge and diverge

f out i,c (k) = ρ i,c (k) ρ i (k) f out i (k) ∀i ∈ {i 1 , i 2 } (48) f out 2 (k) f out 1 (k) δ2 (k) δ1 (k) f out 1 (k) = σ1 (k) -f out 2 (k) P 1 P 2 max P 1 P 2 min (a) 2 × 1 junction f out 2 (k) f out 1 (k) δ2 (k) δ1 (k) C1 C2 P 1 P 2 max P 1 P 2 min (b) 2 × 2 junction
Let i 1 , i 2 be the incoming links and j 1 , j 2 be the outgoing links. To simplify the notation, we use the following shorthand:

• drop the time index k

• δ 1 = δ i1 , δ 2 = δ i2 • σ 1 = σ j1 , σ 1 = σ j2 • P 1 = P i1 , P 2 = P i2
The aggregate split ratios are computed as follows.

β ij (k) = 1 ρ i (k) C c=1 ρ i,c (k) β ij,c (k) ∀(i, j) ∈ {1, 2} × {1, 2} (49)
As shown in Figure 10b, we once again distinguish 3 cases based on where the priority vector intersects the demand constraints. The cases are ) can be obtained by plugging substituting δ 1 (and resp. δ 2 ) for f out 1 (and resp. f out 1 ) in the equations C j for the supply constraints.

C j : β 1j f out 1 + β 2j f out 2 = σ j (50)
Thus, we obtain:

P 1 P 2 max = δ 1 min δ 2 , σ1-β11δ1 β21 , σ2-β12δ1 β22 (51)
P 1 P 2 min = min δ 1 , σ1-β21δ2 β11 , σ2-β22δ2 β12 δ 2 (52)
We can then compute the flow out of incoming links using the optimization problem in equation ( 23). Since the problem is symmetric in i 1 and i 2 , we just solve it for i 1 = 1.

• If P1 P2 > P1 P2 max
, the solution is the intersection of f out 1 = δ 1 and the most constraining supply constraint (C 1 in figure 10b), and results in the trivial solution of f out

1 = δ 1 . • If P1 P2 < P1 P2 min
, the solution is the intersection of f out 2 = δ 2 and the most constraining supply constraint (C 2 in figure 10b), and therefore we obtain:

f out 1 = min δ 1 , σ1-β21δ2 β11 , σ2-β22δ2 β12 • If P1 P2 min ≤ P1 P2 ≤ P1 P2 max
, the solution lies at the point where the priority vector intersects the most constraining supply constraint ((C 2 in figure 10b)), and therefore we obtain:

f out 1 = min P1σ1 P1β11+P2β21 , P1σ2 P1β12+P2β22
This gives us the explicit following explicit solution for f out 1 at a 2 × 2 junction.

f out 1 =            δ 1 if P1 P2 > δ1 min δ2, σ 1 -β 11 δ 1 β 21 , σ 2 -β 12 δ 1 β 22 min δ 1 , σ1-β21δ2 β11 , σ2-β22δ2 β12 if P1 P2 < min δ1, σ 1 -β 21 δ 2 β 11 , σ 2 -β 22 δ 2 β 12 δ2 min P1σ1 P1β11+P2β21 , P1σ2 P1β12+P2β22 otherwise (53) f out 2 is obtained by symmtery.
Finally, we can now compute the flow out of incoming links by commodity.

f out i,c = ρ i,k ρ i f out i ∀i ∈ {1, 2}, ∀c ∈ C (54)
E Computing the partial derivatives of the constraints with respect to the state variables ( ∂ ∂x )

We first iterate through the three classes of variables density, outflow and inflow. All unlisted derivatives evaluate to zero.

Commodity density ρ i,c (k) from equations (H1a, H1b, H1c, H1d))

∂ρ i,c (k) ∂ρ i,c (k -1) = 1, ∀c ∈ C, ∀i ∈ A \ B, ∀k ∈ [[1, T ]] (55) ∂ρ i,c (k) ∂f in i,c (k -1) = ∆t L i , ∀c ∈ C, ∀i ∈ A \ B, ∀k ∈ [[1, T ]] ( 56
) ∂ρ i,c (k) ∂f out i,c (k -1) = - ∆t L i , ∀c ∈ C, ∀i ∈ A \ S, ∀k ∈ [[1, T ]] (57) Flow in f in i,c from equation (H6) ∂f in j,c ∂f out i,c = β ij,c ∀i ∈ J in z , ∀j ∈ J out z , ∀z ∈ J ∀k ∈ [[1, T -1]] (58)
Computing the partial derivatives of the flow out f out i,c is requires a much more involved process. We begin by computing the following intermediate partial derivatives:

• Computing ∂ ∂ρi,c(k) ρi,c(k)δi(k) ρi(k)
:

From equations (H2a, H2b),

δ i (k) = min (F i , v i ρ i (k)) ∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] min F i , ρi(k)Li ∆t ∀i ∈ B, ∀k ∈ [[0, T f ]] (59)
which gives the following equations:

ρ i,c (k) ρ i (k) δ i (k) =    min ρi,c(k) ρi(k) F i , ρ i,c (k) v i ∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] min ρi,c(k) ρi(k) F i , ρi,c(k)Li ∆t ∀i ∈ B, ∀k ∈ [[0, T f ]] (60)
Using equation ( 1),

∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] ∂ ∂ρ i,c (k) ρ i,c (k) ρ i (k) δ i (k) = (ρi(k)-ρi,c(k)) ρi(k) 2 F i if F i < v i ρ i (k) v i otherwise (61) ∀i ∈ B, ∀k ∈ [[0, T f ]] ∂ ∂ρ i,c (k) ρ i,c (k) ρ i (k) δ i (k) = (ρi(k)-ρi,c(k)) ρi(k) 2 F i if F i < Li ∆t ρ i (k) Li ∆t otherwise (62) Remark 4. If ρ i (k) = 0, then F i > Li ∆t ρ i (k) or F i > v i ρ i (k)
, so the derivatives are well defined.

• Computing

∂ ∂ρ i,c ′ (k) ρi,c(k)δi(k) ρi(k)
: (Note the two different commodities c and c ′ .)

∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] ∂ ∂ρ i,c ′ (k) ρ i,c (k) ρ i (k) δ i (k) = -ρi,c(k) ρi(k) 2 F i if F i < v i ρ i (k) 0 otherwise (63) ∀i ∈ B, ∀k ∈ [[0, T f ]] ∂ ∂ρ i,c ′ (k) ρ i,c (k) ρ i (k) δ i (k) = -ρi,c(k) ρi(k) 2 F i if F i < Li ∆t ρ i (k) 0 otherwise (64) Remark 5. If ρ i (k) = 0, then F i > Li ∆t ρ i (k) or F i > v i ρ i (k)
, so the derivatives are well defined.

• Computing ∂ ∂ρi,c(k) ρi,c(k)σj (k) ρi(k)
and

∂ ∂ρ i,c ′ (k) ρi,c(k)σj (k) ρi(k)
when i = j:

Note that σ j (k) does not contain ρ i (k) terms.

∂ ∂ρ i,c (k) ρ i,c (k) ρ i (k) σ j (k) = (ρ i (k) -ρ i,c (k)) ρ i (k) 2 σ j (k) (65) ∂ ∂ρ i,c ′ (k) ρ i,c (k) ρ i (k) σ j (k) = -ρ i,c (k) ρ i (k) 2 σ j (k) (66)
Remark 6. This partial derivative is only needed in cases where the junction is strictly supply constrained and ρ i (k) > 0, so we can ignore the fact that the derivative is undefined at ρ i (k) = 0.

• Computing

ρ i,c (k) ρ i (k) σ j (k) = 0 (70)
Remark 7. This partial derivative is only needed in cases where the junction is strictly supply constrained and ρ i (k) > 0, so we can ignore the fact that the derivative is undefined at ρ i (k) = 0.

• Computing

β ij (k) = 1 ρ i (k) c ′ ∈C ρ i,c ′ (k) β ij,c ′ (k) ∀k ∈ [[0, T f ]] (71) Let κ ij (k) = c ′ ∈C ρ i,c ′ (k) β ij,c ′ (k) ρ i,c (k) ρ i (k) σ j (k) β ij (k) = ρ i,c (k) σ j (k) κ ij (k) (72) ∂ ∂ρ i,c (k) ρ i,c (k) ρ i (k) σ j (k) β ij (k) = κ ij (k)σ j (k) -ρ i,c (k) σ j (k) β ij,c (k) κ ij (k) 2 (73)
Remark 8. This partial derivative is only needed in cases where the junction is strictly supply constrained and ρ i (k) > 0, so we can ignore the fact that the derivative is undefined at ρ i (k) = 0.

• Computing

∂ ∂ρ i,c ′ (k) ρi,c(k)σj (k) ρi(k)βij (k) : ∂ ∂ρ i,c ′ (k) ρ i,c (k) ρ i (k) σ j (k) β ij (k) = -ρ i,c (k) σ j (k) β ij,c ′ (k) κ ij (k) 2 (74)
Remark 9. The ρ i (k) = 0 condition is just as in the previous case.

• Computing

∂ ∂ρ j,c ′ (k) ρi,c(k)σj (k) ρi(k)βij (k) : ∂ ∂ρ j,c ′ (k) ρ i,c (k) ρ i (k) σ j (k) β ij (k) = 1 β ij (k) ∂ ∂ρ i,c ′ (k) ρ i,c (k) σ j (k) ρ i (k) (75)
which can then be simplified using equations (69, 70).

Now we can proceed to computing the partial derivatives of f out i,c . Definition 19 (Demand-constrained junction). A junction is demand-constrained if the flow through the junction is limited by the incoming flow of cell i. We denote this condition by DC(i).

Definition 20 (Supply-constrained junction). A junction is supply-constrained if the flow through the junction is limited by the outgoing flow into some outgoing cell j. We denote this condition by SC(j).

Solution for 1 × 2 junctions

The solutions to all the partial derivatives that appear in the expressions below have already been solved explicitly.

From equation (H5a), if SC(j) (78)

f out i,c = ρ i,c ρ i min σ j (k) β ij (k) , ∀j ∈ J out z | β ij (k) > 0 , δ i (k) ∀z ∈ J 1×n , ∀i ∈ J in
∂f out i,c ∂ρ j,c ′ (k) = 0 if DC(i) ∂ ∂ρ j,c ′ (k) ρi,c(k) ρi(k) σi(k) βij (k) if SC(j) (79)
Remark 10. It is important to note that these derivatives are undefined if the junction is both supply and demand-constrained. However, this can only occur if the density of the cell is exactly equal to the value at which the demand and supply constraints meet. This is extremely unlikely in practice with floating point numerical operations. In the rare event that it does occur, we assume that the junction is supply-constrained.

Solution for 2 × 1 junctions

The solutions to all the partial derivatives that appear in the expressions below have already been solved explicitly.

From equation (H5b)

f out i,c = ρ i,c ρ i     
δ i if P i min δ i + δ i , σ jδ i > δ i P i min δ i + δ i , σ jδ i if P i min δ i + δ i , σ jδ i > δ i P i P i min δ i + δ i , σ j otherwise ∀z ∈ J 2×1 , ∀i ∈ J in z (80) The solution for the 2 × 2 junctions can be obtained using a similar set of computations, but is omitted here for readability and due to length constraints. This concludes the computation of all the partial derivatives required for computing the gradient of the system using the discrete adjoint method.

Remark 11. If we do not have closed form solutions for the junctions, it may not be possible to compute the explicit partial derivatives of the outgoing flow with respect to the partial densities of the incoming and outgoing links of the junction. However, for any junction that cannot be solved explicitly, it is still possible to compute ∂f out i,c (k) k) for all i, i ′ ∈ J in z ∪ J out z and c, c ′ ∈ C with a finite differences method using |J in z | • |J out z | • |C| local simulations of just the junction dynamics (not the entire system). We can then continue to use the adjoint method while numerically differentiating these junctions that do not admit an explicit solution. This local finite differences method can still allow a very efficient computation of the gradient for heterogeneous networks with some junctions or sub-networks that contain complex dynamics.

∂ρ i ′ ,c ′ (

Figure 4 :

 4 Figure4: The synthetic network adapted from[START_REF] Ziliaskopoulos | A linear programming model for the single destination system optimum dynamic traffic assignment problem[END_REF]. There are 10 cells marked C = 1, • • • , 10 and the jam density ρ jam of each cell is denoted by the maximum number of vehicles (N ), since the length of each cell is normalized to 1. The edge weights represent the max flow F between the neighboring cells.

Figure 5 :

 5 Figure 5: The change in total travel time vs percentage of vehicles that can be rerouted. All performance measures are with respect to the local minimum found by the optimizer.

Figure 6 :

 6 Figure 6: The Interstate 210 sub-network

Figure 7 :

 7 Figure 7: The density evolution along the 14 freeway road links with; (a) no incident, (b) a two lane capacity drop from minutes 10-30 at link 5, and (c) flow being rerouted to the parallel arterial due to the capacity drop, and (d) the density difference between the incident profiles with and without rerouting.

  Figure 8: A comparison of the density evolution for different rerouting capacities on the parallel arterial route; (a) only 40% of the arterial capacity can be utilized for rerouting freeway vehicles, (b) only 50% of the arterial capacity can be utilized for rerouting freeway vehicles, and (c) The entire arterial capacity is utilized for rerouting freeway vehicles (i.e. the parallel arterial temporarily closed for other traffic).

Figure 9 :

 9 Figure 9: The base 10 logarithm (log 10 ) of the total computation time for solving the I-210 network vs the number of time steps in the problem. The total time horizon is fixed, so a larger number of time steps implies a smaller time discretization. This also results in a larger number of cells (smaller in length) due to the CFL condition.

  . The system is implemented in the Connected Corridors system at UC Berkeley, a partnership between the California Department of Transportation (CalTrans) and California Partners for advanced transportation technology (PATH) for Integrated Corridor Management (ICM).

  Max flow capacity of cell i P ij Merge priority factor from cell i to cell j Sets J in z Incoming links to junction z J out z Outgoing links to junction z A

  0) Initial density of commodity c on cell i β ij,c (k) Split ratio for commodity c on cell i to cell j, time step k D (o,s) (k) Demand rate of controllable agents going from o ∈ B to s ∈ S, time step k D i,c (k) Demand rate of non-controllable agents of commodity c ∈ N C on cell i, time step k Variables

Algorithm 1

 1 Update buffers for b = 1 to B -1 do ∆l b = min Ll b i (k + 1) ∆l b ∀c ∈ C ⊲ Per-commodity flow leaving buffer b + 1 end for ∆d = di (k) ⊲ Total demand at time step k + 1 b = B ⊲ Allocate the demand to last buffer and create new buffers if needed while ∆d > 0 do ∆l b = min Ll b i (k + 1) , ∆d ⊲ Maximum flow that can enter buffer b l b i,c (k + 1) = l b i,c (k + 1) + d i,c (k+1) d i (k+1) ∆l b ∀c ∈ C ⊲ Per-commodity flow entering buffer b b = b + 1 end while

Figure 10 :

 10 Figure 10: Graph displaying the constraints for different junctions.

  min F i , w j ρ jam j ρ j (k) ∀j ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] w j ρ jam j ρ j (k) ∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] ρi,c(k) ρi(k) F j ∀i ∈ S, ∀k ∈ [[0, T f ]] (68) ∀i ∈ A \ (B ∪ S), ∀k ∈ [[0, T f ]] ∂ ∂ρ j,c ′ (k) ρ i,c (k) ρ i (k) σ j (k) =    0 if F j < w j ρ jam j ρ j (k) -ρi,c(k) ρi(k) w j otherwise (69) ∀i ∈ S, ∀k ∈ [[0, T f ]] ∂ ∂ρ j,c ′ (k)

case 1 :

 1 Pi (min (δi + δi, σj) δi) > δiPi Pi (min (δi + δi, σj)δi) > δiPi c ′ (k) ρi,c(k)δi(k) ρi(k) if δ i + δ i < σ j ∂ ∂ρ i,c ′ (k) ρi,c(k)σj (k) ρi(k) ∂ρ j,c ′ (k) = 0 if δ i + δ i < σ j ∂ ∂ρ j,c ′ (k) ρi,c(k)σj (k)

  is called the demand allocation for commodity c at time step k. The number of vehicles with origin o and destination s that are allocated to commodity c at time step k is D (o,s)

Table 1 :

 1 Optimal allocation of demand across routes

	Time step	1	2	3	Time step	1	2	3	Time step	1	2	3
	Route 1	0.75	0.5	0.5	Route 1	0.25	0.417	.417	Route 1	0	0.25	0.25
	Route 2	0	0	0	Route 2	0	0	0	Route 2	0.25	0	0
	Route 3	0.25	0.5	0.5	Route 3	0.75	0.583	0.583	Route 3	0.75	0.75	0.75
	(a) normal operation		(b) incident local minimum	(c) incident global minimum

For notational simplicity, the equations given here are for the single buffer model.

We have modified the network to increase the capacity of links 8 and 9 such that they can accommodate flow-from both route 2 and route 3, and changed some of the other parameters to satisfy the CFL conditions in section 2.2.

See[START_REF] Morton | Numerical solution of partial differential equations: an introduction[END_REF] for a detailed analysis of finite difference methods.

Acknowldegements

The authors have been supported by the California Department of Transportation under the Connected Corridors program, NSF CAREER grant CNS-0845076 under the project 'Lagrangian Sensing in Large Scale Cyber-Physical Infrastructure Systems', the European Research Council under the European Union's Seventh Framework Program (FP/2007-2013) / ERC Grant Agreement 257661, the INRIA associated team 'Optimal REroute Strategies for Traffic managEment' and the France-Berkeley Fund under the project 'Optimal Traffic Flow Management with GPS Enabled Smartphones'.


Rodrigo Cofre
email: rodrigo.cofre_torres@inria.fr

Bruno Cessac

Can we hear the shape of a Maximum entropy potential from Spike trains ?

come

Can We Hear the Shape of a Maximum Entropy

Potential From Spike Trains?

Maxent is a method for finding the probability distribution that is consistent with known constraints expressed in terms of averages of patterns, but is otherwise as unbiased as possible.

(i) It assumes stationarity in the data; (ii) The choice of constraints is ad-hoc; (iii) The interpretation of the coefficients shaping the Gibbs distribution (Lagrange multipliers) is controversial.

Our results allow to establish a link between synaptic connectivity, stimulus and effective connectivity. Our method works as well for finite size rasters, where transition probabilities can be estimated from empirical averages. Our method opens up new possibilities which allow a better understanding of the role of different neural In order to relate spike trains to the neural tissue properties and stimulus that generate these spiking activity, dynamical neural networks modeling is neccesary.

  Estimating Maximum Entropy distributions from Periodic Orbits in Spike Trains, R. Cofre, B. Cessac, Research Report, open access through HAL-INRIA "Hearing the Maximum Entropy Potential of a spike-generating Markov process", R. Cofre, B. Cessac, submitted, open access through HAL-INRIA NeuroMathComp Team Inria Sophia-Antipolis *rodrigo.cofre_torres@inria.fr

Acknowledgements

We consider a spike-generating stationary Markov process whose transition probabilities are known. We show that there is a canonical potential whose Gibbs distribution, obtained from the Maximum Entropy Principle (MaxEnt), is the equilibrium distribution of this process. We provide a method to compute explicitly and exactly this potential as a linear combination of spatio-temporal interactions. The method is based on the Hammersley Clifford decomposition and on periodic orbits sampling. An explicit correspondence between the parameters of MaxEnt model and the parameters of Markovian models like the Generalized-Linear Model can be established. We also provide numerical results.

This work was supported by the French ministry of Research and University of Nice (EDSTIC), INRIA, ERC-NERVI number 227747, KEOPS ANR-CONICYT and European Union Project # FP7-269921 (BrainScales), Renvision # 600847 and Mathemacs # FP7-ICT-2011.


Rubén Herzog

Joaquin Araya

Michael Pizarro

Bruno Cessac

Cesar Ravello

Maria-Jose Escobar

Adrian G Palacios

Higher Probability of Spatio-Temporal Pairwise Correlation in NM Methods Modeling Spatio-Temporal Pairwise Correlations Fitting the Coefficients (λ) of the Gibbs Potential to Empirical Data. Population Selection and Modeling Pairwise Spatio-Temporal Correlations

Three retina patch were extracted from 2 different O.degus eyes and then recorded using MEA-256 (MultiChannel Systems) during i) spontaneous activity (SA) (5-10 min), ii) WN (30 min, 60 fps, 50μm) and iii) NM (30 or 60 s sequence, 30 trials, 60 fps). Spike sorting were carried out off line using OfflineSorter (Plexon inc.), detecting spikes on each channel with a threshold of -5 S.D. of the signal mean. Then, semi-automatic sorting were performed using ValleySeek or T-Dist E-M algorithm, fitting each class manually. We detect ~ 200 units per patch, whose STA was calculated, obtaining the temporal and 2D spatial filter. RGC were classified according their temporal filter, choosing a group of 14-16 units sharing the same time-to-peak and a monophasic OFF temporal profile.

The selected populations were then analyzed using ENAS software, which fits a Gibbs distribution with spatio-temporal constraints (4) to the empirical data using the maximum-entropy principle. It uses as constrains the empirical spike rate, the pairwise correlations in a 20ms time window (spatial pairwise correlation) and the pairwise correlations with a delay of one time window (temporal pairwise correlation).The free terms fitted are the Gibbs coefficients, which are multipliers (λ) for every observable event (ω) of the Gibbs potential: D is the successive temporal windows in which an observable, ω, occurs. In this case, we use D = 2. i and j are the neuron indexes. Then, we computed a symmetric variant of the Kullback-Leibler divergence (5) and took only the real part to make comparisons between conditions and between observables types.

Background

Goals Results

Acknowledgment

Population retinal responses to natural stimuli is still an open research field, nevertheless retinal ganglion cells (RGC) spatiotemporal correlations may be related to population coding (1). Here, we use a diurnal rodent retina (O. degus), which has the advantage of present a 3:1 proportion of rods and cones (2), respectively, to study the RGC population responses to habitat-based natural stimuli. In order to do this, we have developed a mobile robot that is capable to record movies in the natural habitat of this rodent, simulating both his movements and the eye-ground distance, which allows us to stimulate and record an in vitro retina patch using MEA (multi electrode array) with a sequence of images taken from the animal natural habitat. The analysis of spike statistics has been done using the Enas software (http://enas.gforge.inria.fr/v3/) allowing to characterize spatio-temporal pairwise correlations, in contrast to e.g. Ising models, that analyze only instantaneous spatial correlations.

>

Classify and select a RGC population according to their temporal filter profile.

>

Study the spatio-temporal pairwise correlations of the selected RGC population modeling their activity during WN and NM using Gibbs potential.

> Look for statistical differences/equivalences between the model for NI and WN.

Using our robot we obtained a movie of the natural habitat composed of trees, grass, mountains, blue sky, optic flow in several directions and different simulations of the head movement. Natural movies (NI) presents high spatial and temporal correlation (Figure 1c, 1d)(3) whereas white noise (WN) presents no spatial or temporal correlation (by construction). We found that Spike rates coefficients distributes along the identity, showing that the is no big changes in the firing rate. Spatial and Temporal correlation coefficients are distributed mostly under the identity, which means that there is higher probability to find spatio-temporal correlated events in NM condition than in WN. So, we suggest that the statistical properties of the natural stimuli are increasing correlated events in the studied population and this is in agreement with K-Pairwise models for both conditions (6).

DegusPrime: DegusPrime: habitat recording robot habitat recording robot

The inset panel shows the closeness between the distribution of observables types coefficients. We use a symmetric variant of Kullback-Leibler divergence to compute the closeness between two distributions.

Conclusions

Thanks to Prof. Patricio Orio for providing the computational cluster. Thanks to Selim Kraria for the assistance and support while running ENAS on cluster. Thanks to Monica Otero and Carlos Sepulveda for the help on STA analysis.

  Figure 1. a) White noise stimuli sequence and the response of a RGC population to 120s of the ~1800s complete sequence. b) One frame of the ~1900 frames of the natural movie and the responses of the same RGC population during 120s, when the movie was repeated 4 times. Spatio-temporal correlated firing patterns are visible. c) Joint-pixel intensity distribution of all pixels and the neighbor ones in one frame of the NM. d) Joint-pixel intensity distribution of one pixel and the same one along successive frames of the NM.

>

  Financial support: FONDECYT 1110292, 1140403, ANR-47 CONICYT and Millennium Institute ICM-P09-022-F. CONICYT Fellowship to RH. Plots of the Gibbs Coefficients. Blue diagonal line is identity. The two panels show WN and NM vs IN, where we found that the Spike rates coefficients are bigger for the two visual stimulation conditions. It also shows clustering above and below the identity Since IC is an scotopic/mesopic condition, the inputs for the RGC are different, so it becomes difficult to compare with the photopic conditions. However, those scatter plots are useful to easily detect changes in the probability of the events between 2 conditions. c) Closeness Between Gibbs Coefficients Distributions and its Event Types. The two photopic conditions (WN and NM) has closer distributions, in contrast with IC, which is farthest from the other two. Even though All Coefficients of IC distribution deviates from the photopic ones, they are closer in Spike rates coefficients. In contrary, correlated coefficients are farther. There is big closeness between WN and NM, being the Spatial correlation coefficients which are farthest (inset panel in bStatistical properties of visual stimuli modulates the pairwise correlations of a population of RGC which shares common characteristics, such as their spatial and temporal receptive field filters. > Spatial correlated spiking patterns represents the mayor deviation between WN and NM conditions. Thus, we purpose that population coding for this monophasic OFF RGC population is mostly (but not only) based on spatial correlation when stimulated with Natural Movies.

Tkačik, G., et

	2006). Nature Reviews. Neuroscience, 7(5), 358-66. doi:10.1038/nrn1888.
	2. Jacobs, G. H., et al. (2003). Journal of Comparative Physiology. A, 189(5), 347-61. doi:10.1007/s00359-003-
	0408-0.
	3. Simoncelli E. P., Olshausen	B. A.	(2001). Annu. Rev. Neurosci. 2001. 24:1193-216.
	doi:10.1146/annurev.neuro.24.1.1193.
	4. Nasser, H., & Cessac, B. (2014). Entropy, 16(4), 2244-2277. doi:10.3390/e16042244.
	5. Johnson D. H. , Sinanovic S.(2000) IEEE Transactions on Information Theory. Doi:10.1.1.26.2327.
	6.

1. Averbeck, B. B. et al. (al. (2014). PLoS Computational Biology, 10(1), e1003408. doi:10.1371/journal.pcbi.1003408