Optimal control of urban air pollution related to traffic flow in road networks

Air pollution is one of the most important environmental problems nowadays. In large metropolitan areas, the main source of pollution is vehicular traffic. Consequently, the search for traffic measures that help to improve pollution levels has become a hot topic today. In this article, combining a 1D model to simulate the traffic flow over a road network with a 2D model for pollutant dispersion, we present a tool to search for traffic operations that are optimal in terms of pollution. The utility of this tool is illustrated by formulating the problem of the expansion of a road network as a problem of optimal control of partial differential equations. We propose a complete algorithm to solve the problem, and present some numerical results obtained in a realistic situation posed in the Guadalajara Metropolitan Area (GMA), Mexico.


1.
Introduction. Air pollution is one of the most troubling environmental problems in the world today, where the untimely death of more than three million people per year has already been noted for this reason [17]. In big cities, one of the main causes of this pollution is the vehicular traffic. Governments have become aware of the problem, and we are currently witnessing the incessant search for measures to address the situation: for instance, promoting the use of public transport and electric vehicles, restricting the number of cars that can enter the center of the cities (entry restricted, for example, by even/odd number plate), applying the prohibition of entry to those cars that most pollute, etc.
In large metropolises, atmospheric pollution is fully monitored, and all of them have action protocols that are applied when pollution levels exceed certain thresholds considered to be hazardous. Even so, the effectiveness of the measures that are being taken is in question and, to help analyze and design actuation techniques, numerous numerical simulation and optimal control works have been developed. The use of models based on partial differential equations is common, both in the study of traffic flow in an urban network [8,9,15,13,12], and in the analysis of atmospheric pollution [1,10,23,25]. However, the studies that combine both aspects are not so common, and those who do it (usually only interested in simulating the contamination caused by the traffic, as can be seen, for example, in [23] or [4]), they usually start from a previously known vehicular flow, which makes difficult the search for actuations on the road that are optimal in terms of contamination.
Recently, combining a 1D model for vehicular flow with a 2D model for the dispersion of pollutants, the authors have proposed ( [3]) a methodology that allows to simulate how affects the atmospheric contamination any action that could be taken on the model of traffic. In order to try to illustrate how this coupled model can become a useful tool in the search for optimal actuations, in this work we address the problem of the expansion of a road network: for a given network, with congested traffic problems, studies are being carried out on where to build and how to manage a new avenue, that will help to reduce (minimize) all possible existing pollution. This work can be considered within the framework of the optimal control of partial differential equations, and, with that idea in mind, we start (Section 2) by reformulating the model proposed in [3], in order to give a state system that explicitly includes all the design variables (controls) of the problem. From this system, in Section 3 we formulate the problem in a well-posed mathematical way and, using classical adjoint state techniques [19], we derive an alternative expression of the objective function, which will be very useful in later computations. We propose a complete algorithm to solve the problem (Section 4) and, finally, in Section 5, we present some numerical results obtained when applying the proposed methodology in a realistic situation posed in the Guadalajara Metropolitan Area (GMA), the second largest metropolitan area in Mexico, with a total population of more than four million inhabitants, and a car fleet of more than two million vehicles. 2. Mathematical model for air pollution due to vehicular traffic. We consider an urban domain Ω ⊂ R 2 including a road network composed of N R unidirectional avenues (segments) that meet at a number N J of junctions, such that the endpoints of each segment are either on the boundary of Ω or touch one of the junctions. Each segment A i ⊂ Ω, i = 1, . . . , N R , is modelled by an interval [a i , b i ], and we denote by: a parametrization of the segment A i preserving the sense of motion on the avenue. For simulating the pollution due to traffic flow, we are going to use a mathematical model combining a 1D traffic flow model with a 2D air pollution model (see [3]). Next, we briefly present this combined model.
2.1. The traffic flow model. In order to model the traffic flow in a road network we consider the conservation law formulation proposed by Lighthill and Whitham [18], and by Richards [24], the so-called LRW model. Then, we denote by ρ i (s, t) ∈  t) [number of cars/h] the corresponding flux. Moreover, we denote by N J the number of junctions and, for each junction j ∈ {1, . . . , N J }, we denote by n j and m j the number of incoming and outgoing avenues, respectively. So, we assume that the incoming roads are A νj (1) , . . . , A νj (nj ) , and that the outgoing ones are A νj (nj +1) , . . . , A νj (nj +mj ) .
This model is based on assuming that the velocity v i depends exclusively on the density ρ i . In this case, it is assumed that there exist a Lipschitz continuous and concave function f i : [0, ρ max i ] → R (the so-called static relation), satisfying the following properties: In this situation, the conservation law of the number of vehicles is described by the following system: for i = 1, . . . , N R , j = 1, . . . , N J , k = 1, . . . , n j , l = 1, . . . , m j , and y, z ∈ {1, . . . , N R } denoting, respectively, the indices of incoming and outgoing avenues: β j kl f νj (nj +l) (ρ νj (nj +l) (a νj (nj +l) , .)) = α j lk f νj (k) (ρ νj (k) (b νj (k) , .)) in (0, T ), (2e) where: • Function ρ 0 i represents the initial vehicles density at road A i . • Function ρ in s (respectively, ρ out t ) corresponds to the inlet (respectively, the outlet) density condition for the incoming avenue A y (respectively, the outgoing avenue A t ).
• Parameters α j lk represent the preference of drivers arriving to a junction, that is, α j lk gives the percentage of drivers that, arriving from the incoming avenue A νj (k) , are going to take the outgoing avenue A νj (nj +l) . Consequently, the following constraints should be verified: 0 ≤ α j lk ≤ 1, and, for k = 1, . . . , n j , mj l=1 α j lk = 1.
• Parameters β j kl represent the ingoing capacity in outgoing avenues, that is, β j kl gives the percentage of vehicles that, coming from the A νj (k) , can enter the outgoing avenue A νj (nj +l) . In a similar way to previous case, these parameters should verify: 0 ≤ β j kl ≤ 1, and, for l = 1, . . . , m j , nj k=1 β j kl = 1.
Remark 1. Boundary condition (2d) may be not required in any of the outgoing avenues, assuming free traffic through that exit. In this case, it is assumed that avenues leaving the domain are infinite. Consequently, those avenues are modelled by the infinite segment [a z , ∞) (that is, equations (2a) and (2b) should be posed, respectively, in (a z , ∞) × (0, T ) and in [a z , ∞)), although we are only interested in the behaviour of the model along the finite interval (a z , b z ).
Remark 2. Condition (2e) is related to the conservation of cars in junctions. In fact, taking into account (3) and (4), the sum in l (from 1 to m j ) and in k (from 1 to n j ) gives the classical Rankine-Hugoniot relation (see [8,9,15]): which states the continuity of the total flux across a junction, and thus the conservation of the mass of the vehicles through the junction itself.
Existence of solution for system (2a)-(2d) is guaranteed in different situations. For example (see [3] and the references therein), if ρ 0 i ∈ L ∞ (a i , b i ), ρ in k ∈ L ∞ (0, T ), and for each junction j we assume one of these situations: 1. For case n j > m j = 1 (usually known as on-ramp), the following additional conditions must be satisfied: if not all cars can go through the junction, then is the capacity of the outgoing avenue. 2. For case n j ≤ m j (the number of outgoing roads is greater that the number of incoming roads), matrices A j = (α j lk ) must verify a technical hypothesis (see, for example, [3,8]), and the following additional conditions must be satisfied: (a) f νj (nj +l) (ρ νj (nj +l) (a νj (nj +l) , .)) = nj k=1 α j lk f νj (k) (ρ νj (k) (b νj (k) , .)), , .)) is maximum, subject to previous property (a), then, the existence of, at least, one solution ρ = (ρ 1 , . . . , ρ N R ) of the traffic flow model can be proved: Proof. The existence of, at least, one solution was proved initially by Coclite, Garavello and Piccoli [8], assuming junctions with at most two incoming roads and two outgoing roads. However, it was later extended to the general case by Garavello and Piccoli [9].
2.2. The air pollution model. In this Subsection we are interested in simulating air pollution due to vehicular traffic. Then, we assume that vehicle emissions depend on the traffic density given by the solution of model (2a)-(2e). Implicitly, this means considering avenues as curves in Ω (i.e., sets of zero measure) and, consequently, emissions on the avenues should be included in the model through measures on these sets (see [22]).
So, the concentration φ(x, t) [kg/km 2 ] of an air pollution indicator (for instance, in our particular case, carbon monoxide CO), is given by the following initial/boundary value problem [10,25]: where is the CO extinction rate corresponding to the (first order) reaction term, φ 0 is a known function giving the initial CO concentration, n denotes the unit outward normal vector to the boundary ∂Ω = stands for the source of pollution due to vehicular traffic on the avenue A i , and it is the Radon measure (that is, an element of the dual of the space of continuous functions) given by: where G Ai (., t) represents the vehicular emissions on the avenue A i . If we assume that emissions are proportional both to traffic density and to its flux, for the parametrization given by (1) , with ρ i the solution of model (2a)-(2e), and γ i and η i weight parameters representing contamination rates. In this case, and assuming the functions to be smooth enough, we have that: As a consequence of Theorem 2.1, we have the bounded solutions ρ i ∈ C(0, T ; Moreover, functions f i are continuous and, consequently, measures given by (6) are well defined. In this way, the second member of the partial differential equation (5a) is a well-defined Radon measure. So, under standard regularity hypotheses on the problem data (mainly, v ∈ (L ∞ (Ω×(0, T ))) 2 , µ, κ ∈ L ∞ (Ω × (0, T )), φ 0 ∈ L 2 (Ω)), the existence of a unique weak solution for the air pollution model can be demonstrated: The problem (5a)-(5d) has a weak solution φ, which is the unique transposition solution of (5a)-(5d). Moreover, φ ∈ L r (0, T ; W 1,p (Ω)), ∀r, p ∈ [1, 2) such that 2 r + 2 p > 3. Proof. The existence of a unique transposition solution for the problem can be obtained by using the technical reasonings introduced in Casas [5], also employed in Casas et al. [6,7]. The full details can be found in a recent paper by the authors [20], where it is also shown that this transposition solution is a weak solution.
3. Construction and management of a new avenue: The expansion of a road network as an optimal control problem. The coupled model (2a)-(2e), (5a)-(5d) is very useful in order to study what consequences will have, in terms of atmospheric pollution, any action taken on vehicular traffic. For instance, if we want to observe the incidence on pollution, of the actuation relating to restricting the entry into the network of cars with even/odd plate number, we will only have to solve that coupled model, replacing the inlet boundary conditions (2c) for the traffic flow model with the new inlet conditions: and then to analyze the values obtained for pollution concentration φ (cf. [11]), making a comparison with the values previously obtained for the original boundary conditions (2c). The model (2a)-(2e), (5a)-(5d) can also be very useful for designing operations on the road network, which are optimal in terms of pollution control. Consider, for example, a situation in which the network of roads is intended to be expanded, building a new exit avenue that decongests the current traffic in a given area. To fix ideas, we take as domain Ω ⊂ R 2 the one shown in Figure 1, designed to study atmospheric pollution in the Guadalajara Metropolitan Area (GMA), Mexico (see further details in [1]). In this domain we consider a network that initially consists of 12 unidirectional avenues which intersect in 6 nodes (see Figure 1). It is assumed that the exit avenues A old 11 and A old 12 are suffering from congested traffic, and is under study the construction of a new exit road (A 13 , which joins the zone Z N EW with a given point q f ix ∈ ∂Ω, where that artery will connect with an existing highcapacity pathway). In addition, in order to decongest the avenues A old 11 and A old 12 , it is intended to construct, on each of them, a node that divides them into two, so that A old 11 = A 11 ∪ A 14 and A old 12 = A 12 ∪ A 15 . From these nodes (nodes 8 and 7, respectively) there will be two new avenues, A 16 and A 17 , (one from each node) that will intersect at another node (node 9), located at a point p ∈ Z N EW , from which the new exit road is started (see Figure 1). In this way, the new expanded network will have N R = 17 unidirectional avenues that intersect in N J = 9 nodes.
Starting from the already existing network and the known point q f ix ∈ ∂Ω, if we denote by [a old 11 , b old 11 ] and [a old 12 , b old 12 ] the intervals representing the old avenues A old 11 andA old 12 , the new expanded network is uniquely determined by the point p = (p 1 , p 2 ) ∈ Z N EW where the node 9 will be located, and by the parameters s 11 ∈ [a old 11 , b old 11 ] and s 12 ∈ [a old 12 , b old 12 ], which determine the points in which the old avenues A old 11 and A old 12 will be split (that is, the points at which nodes 8 and 7 will be installed).
On the other hand, in order to determine the traffic that will circulate through the new avenues, it will also be necessary to determine the values of the parameters α 7 , α 8 ∈ [0, 1], that determine the preferences of the drivers arriving at nodes 7 and 8, and the parameter β 9 ∈ [0, 1], which indicates how the traffic in node 9 is regulated. In the node j = 9 there are two incoming and one outgoing roads (i.e., n 9 = 2 and m 9 = 1). So, in this case, Thinking about expanding the network in order to reduce pollution in the domain Ω we define: • The vector of decision variables u (the control) that collects the variables (both of construction and management type) that determine the operation of the new network, which is given by: u = (p 1 , p 2 , s 11 , s 12 , α 7 , α 8 , β 9 ) ∈ R 7 .
Then, Green's formula can be used, and we obtain the following equalities giving the expression (10):

Numerical solution.
To solve the problem (8) very different numerical methods can be used. In this paper, since no information on the objective function gradient is available, and since the dimension of the control variable u is relatively low, a classical direct search algorithm has been used, the simplex algorithm of Nelder and Mead [21]. In order to use the Nelder-Mead method, the constrained problem (8) needs to be previously transformed into an unconstrained one by using, for instance, a penalty function. To do it, we introduce the function Ψ : U ad −→ R 14 collecting all the constraints corresponding to the characterization of U ad , in such a way that u ∈ U ad ⇔ Ψ(u) ≤ 0. For instance, if the admissible region Z N EW is given by the rectangle Z N EW = [p min where > 0 determines the relative contribution of J and the penalty term. In this way, Φ is an exact penalty function (cf. Han [14]) in the sense that, for sufficiently small , the solutions of our constrained problem (8) are equivalent to the minimizers of function Φ.
In order to minimize Φ, since the computation of Ψ(u) is completely straightforward, the only difficulty in applying the method relies in the computation of J(u). So, to evaluate J, we must first solve the adjoint system (9a)-(9d). For this, in this work we have used a method which combines characteristics for the time discretization with a Lagrange P 1 finite element method for the space discretization (see [10] for more details). Once the adjoint state g has been obtained, whenever we want to evaluate J(u), all we have to do is solving the traffic flow system (2a)-(2e) in the road network determined by u, and applying the expression (10) derived in Theorem 3.1. Therefore, we continue this Section by describing the method used to solve the system (2a)-(2e).

4.1.
Numerical solution of the traffic flow model on a single avenue: System (2a)-(2d). In this first case, we consider a scenario with only one single avenue A ⊂ Ω, defined from the interval I = [a, b] ⊂ R. So, we consider the system (2a)-(2d) without subscripts, and we are going to approximate its solution ρ, using a classical finite volume approximation of the conservation law.
To do this, we divide the spatial domain I = [a, b] into M cells I k = [s k− 1 2 , s k+ 1 2 ] of length ∆s > 0, and denote by s k = (s k− 1 2 + s k+ 1 2 )/2 the midpoint of each cell. The time interval [0, T ] is also divided into N ∈ N subintervals of length ∆t = T /N , and define t n = n∆t. Then, we consider the control volume V ≡ I k × [t n , t n+1 ] and, integrating equation (2a) in this domain, we obtain that with: •ρ n k = 1 ∆s where: • supp : [0, ρ max ] −→ R is the supply function given by • dem : [0, ρ max ] −→ R is the demand function given by Remark 4. Expression (12) for k = M requires the knowledge of dataρ n M +1 , for n = 0, . . . , N − 1. If at the outgoing endpoint there exist available data providing traffic density -that is, if we have boundary condition (2d) -then we can takē ρ n M +1 = ρ out (t n ). If we do not have those data, it is assumed that there is free traffic downstream and, consequently, we can take supp(ρ n M +1 ) = C. The supply-demand method is mainly based on assuming that the flow at the intersection of two cells is the minimum between the flow that upstream cell can supply and the flow that downstream cell demands.
The definitions of supply and demand functions are based on the assumption that the upstream cell can accommodate the maximum flow (C) for free traffic (ρ ≤ ρ C ) and must restrict to actual traffic if this is congested (ρ C ≤ ρ ≤ ρ max ), while the downstream cell demands only the actual flow in this cell if there is free traffic and the maximum flow if traffic is congested. A detailed explanation of the method can be found in the recent monograph of Treiber and Kesting [26].
Remark 5. It is worthwhile remarking here that supp(ρ) and dem(ρ) are monotone functions (non-increasing and non-decreasing, respectively). Thus, the flux f n k+ 1 2 = min{dem(ρ n k ), supp(ρ n k+1 )}, given by (12), is a non-decreasing function with respect to the first argumentρ n k , and non-increasing with respect to the second oneρ n k+1 . This fact makes the numerical method given by (11)- (12) to become a monotone first order scheme and, consequently, if a proper CFL condition is satisfied, we have that (see [8, (12) when defining the flows in the first order numerical scheme (11) results of its ability to be implemented on road networks with unidirectional avenues. So, each avenue A i in the network is divided into a number M i of cells, and the traffic densities at the midpoints of these cells are calculated as described in previous Section. The only difference is when that avenue begins or ends at a junction. If the avenue begins at a junction, then the incoming boundary condition at that point is not known, and the density in the first cell of that avenue (ρ n+1 i,1 , n = 0, . . . , N − 1) needs to be calculated. This computation is done using the expression (11) for k = 1, but in this case it is necessary to specify how to calculate the incoming upstream boundary flux of the first cell (f n i, 1 2 , n = 0, . . . , N − 1) -note that in (12) this flux is not considered.
Similarly, if the avenue ends at a junction, we do not know the outgoing boundary condition at that point, and we cannot assume that downstream there is free traffic, as in the case of an infinite avenue (see Remark 4). Then, in this case it is necessary to redefine the outgoing flux downstream boundary of the last cell (f n i,Mi+ 1 2 , n = 0, . . . , N − 1). In below paragraphs we can find how to define these numerical flows from condition (2e). So, let us consider the junction j ∈ {1, . . . , N J }, where, in order to simplify the notation, we will assume that the n j incoming avenues are A 1 , . . . , A nj and the m j outgoing ones are A nj +1 , . . . , A nj +mj (that is, we suppose that ν(l) = l, for l = 1, . . . , n j + m j ). Then, • For each avenue ending at the intersection (for each k = 1, . . . , n j ) the flow of the last of its cells is calculated as in (12), but taking into account that there are now m j cells that can accommodate vehicles, and that what each of them can support (for each l = n j + 1, . . . , n j + m j ) will be affected by the parameter β j kl given by the percentage of vehicles that, arriving on the avenue A k , can take the avenue A nj +l . In this way: • For each avenue beginning at the intersection (for each l = 1, . . . , m j ) the flow of the first of its cells is also calculated as in (12), but bearing in mind that there are now n j cells that demand to enter, and that, in addition, the demand for each of them is affected by the parameter α j lk which gives the percentage of vehicles that, arriving by the avenue A k , wants to take the avenue A nj +l . In this way: 4.3. Numerical optimization. Finally, once we know how to compute the values of penalty function Φ(u) by the method proposed in above Subsections, we have only to explain how to find the minimizers of this function Φ. In order to compute the minimum values of Φ, as already commented at the beginning of this Section, we propose to use the Nelder-Mead algorithm [21], which allows us to find an unconstrained minimizer, starting from a set of admissible points as a first iterate. This algorithm is a derivative-free method, which merely compares values of function Φ: the values of the objective function, being taken from a set of sample points (simplex), are used to continue the sampling. A detailed description of the simplex algorithm of Nelder and Mead to obtain the minimizers of Φ can be found, for instance, in [2]. Its simplicity and efficiency have made the algorithm probably the most widely cited of the direct search methods.
Although the Nelder-Mead algorithm is a heuristic search method that is not guaranteed to converge in the general case, it presents good convergence properties in low dimension cases. Moreover, to prevent stagnation at a non-optimal point, we use a classical modification: when stagnation is detected, we modify the simplex by an oriented restart, replacing it by a new smaller one. Interested readers can find further remarks on this method, for example, in [2].
5. Numerical results. The method described in above Section has been used to analyze the expansion of the existing road network in the domain Ω ⊂ R 2 (see Figure 1), designed to study the atmospheric pollution of the GMA. A simulation of T = 24 h has been carried out, and, to facilitate the interpretation of the results, it has been assumed that initially there are neither cars in the network (ρ 0 i = 0) nor pollution in the domain Ω (φ 0 = 0). Related to the traffic model, it has been assumed that all avenues of the network (including those to be built) have the same characteristics, so that the theoretical flow is the same for all of them. In order to define this theoretical flow, we have used the well-known Triangular Fundamental Diagram (TFD), which consists of assuming that the function f (ρ) provided by that flow (the so-called static relation) is the triangle shaped, piecewise linear function (see Figure 2), given by: For the boundary conditions, we have assumed an equal input of vehicles by the two entry roads (ρ in 1 = ρ in 2 ), corresponding to a standard weekday, as given in Figure 3, and for the exit avenues, congested traffic in A 14 and A 15 (ρ out 14 = ρ out 15 = 100 u/km), and free (but known) traffic in A 10 and in the new avenue A 13 (ρ out 10 = ρ out 13 = 10 u/km). Finally, for the already existing network, in the management of the intersections (parameters α j lk and β j kl ) it has been assumed that A old 11 and A old 12 were the main exit avenues. . Boundary condition for the LWR model on A 1 and A 2 (functions ρ in 1 (t) = ρ in 2 (t)), corresponding to a weekday with typical peak and valley hours.
In this situation, for the already existing network, the CO concentration isolines, after T = 24 hours, can be seen in Figure 4(a). In terms of contamination, the best network expansion is illustrated in Figure 4(b). The vector of corresponding design variables (solution of the problem (8)), reached after 98 iterations of the Nelder-Mead algorithm, is: u opt = (−10.6318, 11.1209, 6.6155, 3.0508, 0.5012, 0.5253, 0.7025).
In this comparison of the levels of contamination, for the final time of simulation, of the unexpanded network (Figure 4(a)) with those CO levels expected for the optimized network defined by the vector u opt (Figure 4(b)), it can be observed how the traffic distribution determined by the parameters α 7 = 0.5012, α 8 = 0.5253 and β 9 = 0.7025 leads, after T = 24 hours, to a much better distribution of pollution in the western area of the domain. In addition, for this expansion of the network we obtain, by comparing J(u opt ) with the value of J corresponding to the original network, a global decrease in pollution of 7.98%. To illustrate in a clearer way the behaviour and characteristics of the new road network, we have also included a comparison of several average values obtained for the optimized network defined by u opt with those that were for the original unexpanded network. So, in Figure 5 the mean values of the flux (Figure 5(a)) and the velocity (Figure 5(b)), averaged for the whole road network, are compared for the time interval of T = 24 hours, showing how the expansion given by u opt increases traffic fluidity, especially in the rush hour, avoiding that the mean velocity drops to zero at any period of time, and remaining always under the maximum allowed velocity V 0 = 50 km/h. Finally, in Figure 6(a) we can see the evolution along the time interval of the mean car emmisions to atmosphere (averaged on the road network), and in Figure  6(b) the mean CO concentration (in this case, averaged on the whole domain Ω).
In an overview of these figures we observe how the network expansion significantly reduces emmisions to the atmosphere and improves contamination levels as time increases. In addition, if we look closely at Figure 6(a) we can see that, in the first times, the mean emissions are exactly the same in the initial network as in  the extended network. This is because the network is initially empty and cars take a few moments to reach nodes 8 and 7. When cars begin to reach these nodes, and before peak time arrives, there is free traffic in the two output avenues (A 11 and A 12 ), and as those routes are shorter than any of the variants included in the extended network, density and flow (and thus also emissions) are slightly lower in the original network than in the enlarged one. However, when rush-hour traffic arrives, avenues A 11 and A 12 begin to undergo congested traffic and it is from this moment when the emissions are lower in the extended network than in the original one (since the extended network suffers fewer congested traffic problems). 6. Conclusions. In this paper, the mathematical model previously analyzed by the authors in [3] is reformulated to simulate the vehicular traffic on a road network, and the contamination that this traffic generates in a two-dimensional domain, in a way suitable for optimization purposes. This model shows to be very useful for dealing with optimal control problems that seek operations on traffic in order to improve pollution levels.
This fact is made visible, formulating the problem of the expansion of a metropolitan area road network as a problem of optimal control of partial differential equations. It is shown that, from a mathematical viewpoint, the optimal control problem is well posed, the computation of the objective function by means of the adjoint state is rigorously justified and, based on this result, a complete algorithm is proposed to numerically solve the proposed problem. Finally, the applicability of the method is illustrated using it in a realistic scenario, posed in the Guadalajara Metropolitan Area, Mexico.
The complexity of this realistic problem means that other theoretical aspects, such as demonstrating existence/uniqueness of optimal solutions or deriving optimality conditions that allow solving the problem with gradient-type methods, are still open questions in which we intend to proceed in the future. Moreover, collaboration with traffic experts to apply the model in real situations, as well as the optimal management of the road network from an environmental point of view, are other lines of work that are intended to be addressed in the near future.