The Traffic Reaction Model: A kinetic compartmental approach to road traffic modeling

In this work, a family of finite volume discretization schemes for LWR-type first order traffic flow models (with possible on-and off-ramps) is proposed: the Traffic Reaction Model (TRM). These schemes yield systems of ODEs that are formally equivalent to the kinetic systems used to model chemical reaction networks. An in-depth numerical analysis of the TRM is performed. On the one hand, the analytical properties of the scheme (nonnegative, conservative, capacity-preserving, monotone) and its relation to more traditional schemes for traffic flow models (Godunov, CTM) are presented. Finally, the link between the TRM and kinetic systems is exploited to offer a novel compartmental interpretation of traffic models. In particular, kinetic theory is used to derive dynamical properties (namely persistence and Lyapunov stability) of the TRM for a specific road configuration. Two extensions of the proposed model, to networks and changing driving conditions, are also described.


Introduction
Macroscopic traffic flow models play an important role in the analysis, simulation, and control of traffic flows on networks [9,26,27,28,38,41,45,52].In this line of research, first order traffic flow models, that describe the spatiotemporal evolution of the vehicular density through conservation laws, predominate in the literature [21,51].These models, e.g.[14,21,24,34], express the conservation of vehicles on the road by linking the density and the flux of vehicles in space and time.In particular, Lighthill and Whitham [34], and Richards [42] suggested to further assume that the flux of vehicles can be expressed as function of the density called flux function (or fundamental diagram).The resulting traffic model, called Lighthill-Whitham-Richards (LWR) model, has been extensively studied and takes the form of a nonlinear hyperbolic Partial Differential Equation (hPDE) satisfied by the density of vehicles [3,52].In general, no closed-form solution is available to solve such hPDEs and therefore appropriate numerical approximations must be introduced.Finite Volume Methods (FVMs) have been particularly praised for deriving numerical solutions of conservation laws, especially for their ability to accurately describe the singularities specific to these hPDEs (e.g.shock waves) [3,32].The Godunov scheme [22] is an instance of such methods widely used in traffic flow modeling for its advantageous theoretical and convergence properties [32].For instance, Daganzo proposes to equip the traditional Godunov scheme [22] with static saturation functions (capacitated sending-receiving flows between cells), thus yielding the so-called Cell Transmission Model (CTM) [14].This model is convenient in terms of describing capacity changes and has been the basis of many traffic network control solutions e.g.[10].The choice of numerical scheme is instrumental in deriving proper approximations of the solution of the considered hPDEs.First, the numerical approximations at hand need to preserve several important properties of the original hPDE like the nonnegativity of the solution or the fact that it should be bounded by the road capacity at all times.Besides, when it comes to FVMs, their stability [17] can be assessed by analyzing the systems of ordinary differential equations (ODEs) resulting from the space discretization of the PDE.These systems of ODEs may then have a connection to an already existing modeling framework thus facilitating their analysis.On the other hand, positive (or nonnegative) systems having the property that the state variables are always nonnegative have a distinguished role in systems and control theory due to the wide field of possible applications e.g., in thermodynamics, chemistry, biology or economy [18].Moreover, the positivity of the states efficiently supports dynamical analysis and control design [23].Within nonnegative systems, kinetic models (also called chemical reaction networks or CRNs) and the closely related compartmental models are good descriptors of complex nonlinear dynamics, and their mathematical and corresponding graph structure are intuitive and useful to state strong results on qualitative dynamics [16,25].The theory of kinetic systems is a rapidly developing area of increasing interest where the most significant contributions are related to the existence and uniqueness of equilibria and robust stability even when several model parameters (possibly including time delays) are not precisely known or unknown [20].Therefore, it is of general interest to transform originally non-chemical models to kinetic form, and interpret or discover their properties from this point of view [43].Building on these ideas, we propose and study a family of finite volume approximations of the LWR model, called Traffic Reaction Model (TRM), with direct links to the positive kinetic systems used to model chemical reaction networks [20,49].The benefits of this analogy are multiple.On the one hand, the resulting model provides a new interpretation of the discrete dynamics of the LWR model as a compartmental chemical system in which the concentration of "particles" of free and occupied space are exchanged between adjacent road cells.This analogy also allows to straightforwardly extend the LWR model to better represent traffic conditions that can occur in real life (e.g.merging, changing road conditions, traffic lights).Invariance and symmetry properties have extensively been researched in kinematic wave theory e.g.[15,31,35,36,37].In the latter two papers, duality plays a key role to bridge the gap between coordinates in which traffic flow models are presented.A key consequence is that variational theory enables us to connect homogeneous first order traffic flow models with proper initial and boundary value problems.In [30], the authors explore the symmetry of kinematic wave models under coordinates transformation and dual variable changes.When dualized, it is shown that travel time, delay, and flow remain invariant.In our paper, flow invariance is in the epicenter of TRM (through the re-definition of the flux by the dual state variable).The dual state variable is the free space concentration or dual density, which ensures the invariance of the nonlinear flux function via a non-unique parametrization and decomposition.However, unlike in [30], we do not try to explicitly transform the conservation laws and define new coordinates to ease delay or capacity calculations (for triangular fundamental diagrams).Instead, we use a nonlinear flux function, and decompose and parameterize it with the dual state variable (free space), where no coordinates transformation is required.Via the dual state variable and with the help of the flux decomposition, the numerical segmentation of the hyperbolic PDE (TRM) preserves several fundamental properties (e.g.non-negativity, capacity, etc).Furthermore, TRM can handle (under some mild conditions) inhomogeneity in PDEs and can accommodate a wide set of nonlinear flux functions.Finally, the system of nonlinear ODEs obtained by the TRM scheme can be interpreted as particular instances of kinetic dynamical systems which can be formally represented as chemical reactions.This special description opens the way to leverage the extensive literature on reaction networks to derive and analyze the dynamical properties of the system (such as robustness or stability) [5,19,44,46], and even propose novel optimal control strategies and traffic management techniques [8,44] towards [11,12,13].The main contributions of this work are as follows.Firstly, we propose an in-depth analysis of the TRM as a numerical scheme.We prove that it defines is consistent, monotone, nonnegative, capacity-preserving, and conservative numerical schemes.We also draw parallels between the TRM and existing methods, showing in particular that the TRM is in fact equivalent to a generic formulation of Daganzo's CTM, and that the Godunov scheme can be seen as particular instance of TRM.Consequently, both methods can also be treated as particular kinetic systems and the results from kinetic system theory can also be extended to them.Secondly, we provide the arguments allowing to bridge the gap between numerical schemes for traffic models and kinetic systems for chemical reaction networks.We provide detailed physical interpretations of how traffic flow is apprehended when seen as a kinetic system.We then leverage the link between the TRM and kinetic systems to provide an analysis of some dynamical properties of the TRM in a particular case.Persistence is introduced as new concept for traffic flow models based on this parallel.Also the structure and stability of the equilibria of the system are studied using Lyapunov theory for kinetic systems.Finally, we present two direct extensions of the TRM that directly follow from the kinetic interpretation of traffic proposed in this work.The first one accounts for changing road conditions.Such extensions can be used to accurately model real traffic data set in state estimation and short-term prediction applications, as advocated in [39,40].The second extension allows us to describe traffic on road networks.In this aspect, the TRM hints to model the link and node behaviour identically using reaction rate governed mixing compartments, giving rise to a smooth interface between link and node modeling framework.Furthermore, the network extension of TRM fits in the genealogy of first order node models presented in [50].The network extended TRM fulfills the requirements for node models described in [50], where most of the properties are preserved thanks to the applied non-negative discretization scheme.Interestingly, many of the requirements are of continuous nature (e.g.conservation, capacity, non-negativity, or invariance) and no external constraints need to be enforced to guarantee them.It is important to note that the Network TRM presented in this paper covers only non-signalized intersections.The outline of the paper is as follows.In Section 2, we introduce the TRM as a finite volume scheme for the LWR model and describe its main properties.We then expose in Section 3 the coupling that exists between kinetic reaction theory and the TRM.In Section 4, we provide an analysis of the TRM through its interpretation as particular case of CTM and a numerical experiment testing the convergence and accuracy of the TRM.In Section 3.3, we derive some dynamical properties (namely persistence and Lyapunov stability) of the TRM when considered on a circular road.Finally, we present in Section 5 the two extensions of the TRM dicussed above.

Nonnegative discretization of hyperbolic PDEs
The Lighthill-Whitham-Richards (LWR) traffic model [34] is the first-order macroscopic traffic flow model defined by the following partial differential equation (PDE) on the domain This PDE models the space-time evolution of a conserved quantity ρ(x, t) ∈ Ω = [0, ρ max ] corresponding to to the traffic density on a road, where ρ max denotes the maximal vehicle density.The flux of vehicles is represented f • ρ, for some continuous and at least once differentiable f : R → R called flux function.In particular, we denote by ρ crit the critical density value at which the flux attains its maximal value f max (i.e.f max = f (ρ crit )).Finally, the functions r, s model the on-and off-ramps, respectively.Remark 2.1.We consider physically-relevant (weak) solutions ρ of PDE (1) also called entropy solutions.We refer the reader to Appendix A for a result justifying the existence and uniqueness of an entropy solution, and the regularity of this solution under non-zero sink and source terms.
To mirror the vehicle density ρ, we introduce a function ν modeling the density of free space (i.e.unoccupied by a vehicle) along the road, and given by ν = ρ max − ρ.
We set ν max = ρ max .Note then that the free space density ν can be seen as a dual or adjoint variable to the density ρ. 2 In this work, we consider the case where the flux function f has the form where g : Ω × Ω → R + satisfies the following assumptions3 : (A1) g is Lipschitz continuous w.r.t.both arguments, with associated Lipschitz constants K 1 > 0 and K 2 > 0, (A2) g is non-decreasing in each arguments, (A3) g(ρ, 0) = g(0, ν) = 0 for all ρ, ν ∈ Ω (which ensures that no vehicles is removed (resp.added) if the road is empty (resp.at capacity).Example 1.Any flux function f that can be written as where g 1 and g 2 are non-decreasing Lipschitz-continuous functions such that g 1 (0) = g 2 (0) = 0 satisfy Assumptions (A1)-(A3).In particular, if we take g 1 (ρ) = ρ, then the function ρ → g 2 (ρ max − ρ) can be interpreted as the speeddensity relationship of the fundamental diagram.Within this framework, taking g 2 (ν) ∝ ν yields the Greenshields  (quadratic) flux function (cf. Figure 1).We can also retrieve trapezoidal fundamental diagrams by considering g 2 as where v max denotes the free flow speed, and ρ 1 = ρ max − ν 1 , ρ 2 = ρ max − ν 2 denote the critical densities of the fundamental diagram which satisfy 0 < ρ 1 ≤ ρ 2 < ρ max .Figure 1 gives an example of a trapezoidal flux function obtained using this decomposition.
Regarding the source and sink terms, we assume that they can be written for any where • the functions ρ on and ν off are the traffic density of the on-ramp and the free space density of the off-ramp, respectively, and are assumed to be taking values in Ω; • the functions g on : Ω × Ω → R + and g off : Ω × Ω → R + are the traffic flows of the on-and off-ramp, respectively, and are assumed to satisfy Assumptions (A1)-(A3); • the spatial position of the on-and off-ramp is described by indicator functions 1 on and 1 off defined as for some x on ≤ x on and x off ≤ x off .We use the Finite Volume Method (FVM) to spatially discretize (i.e.semi-discretize) the conservation law model (1) with flux (2), source and sink (3) functions.To do so, we start by slicing the road into cells of length ∆x > 0, and centered at the points x i = i∆x, for i ∈ Z.We denote by I the set of road cell indices (hence I = Z for an infinite road) and denote by x i−1/2 and x i+1/2 the left and right boundary of the i-th cell, respectively, and by ] the i-th cell.Then, the FVM scheme can be written as where ρ i is an approximation of the average traffic density over the i-th cell, the so-called numerical flux F is given by F (u, v) = g(u, ρ max − v), and the numerical source R i and sink S i terms are given (for t ≥ 0, i ∈ I) by We call Traffic Reaction Model (TRM) the resulting model which describes traffic along a (discretized) road.Note that this model is control-oriented since changing g and/or g on , g off corresponds to speed limit control and/or ramp metering.We also take as initial condition: As such, the TRM defines a particular instance of finite volume scheme for scalar conservation laws [4,33].Finally, note that in practice, finite roads should be considered.To be able to use the TRM to model traffic on a finite road, boundary conditions describing the incoming and outgoing traffic on the road have to be defined.Assuming that the road is now discretized into P cells with indices i ∈ I = {1, . . ., P }, the density of vehicles on these cells is modeled using the system of ODEs defined by (4) while assuming that the quantities ρ 0 and ρ P +1 are some specific time-dependent functions taking values in Ω.For instance, these functions can be set so that the quantities F (ρ 0 , ρ 1 ) and F (ρ P , ρ P +1 ) match some known incoming and outgoing fluxes in the road.Another possibility is to set ρ 0 = ρ P and ρ P +1 = ρ 1 , which corresponds to periodic boundary conditions, or equivalently a circular road.We now present a few numerical properties of the TRM, which hold whether the road is finite or not (cf.Appendix F.1 for proof).Theorem 2.2.The numerical flux F is consistent with the flux function f (i.e. it satisfies for any u ∈ Ω, F (u, u) = f (u)) and is monotone.As for the TRM defined in (4), it preserves nonnegativity and capacity, and it is conservative finite volume scheme.
It is important to note that, for a given f the choice of g is not unique.Hence, the TRM defines a family of numerical discretization schemes by means of the decomposition of the flux f (ρ) in ( 2), which all share the properties described in Theorem 2.2.
Example 2. Consider the so-called Greenshield flux function defined by where v is the average vehicle speed, and ω = v max /ρ max with the free flow speed v max > 0. Different decompositions of f as in ( 2) can be proposed, among which: are the supply and demand functions, respectively, and f max denotes the maximal value of the flux.Note in particular that the decomposition g Gdnv (ρ, ν) in ( 6) is the Godunov flux, meaning that the Godunov scheme can be seen as a particular instance of TRM.In the remainder of this text, we refer to the decomposition g MAK in (5) as a mass action kinetic (MAK) decomposition.We further discretize the TRM by considering a time step ∆t > 0 and taking t k = k∆t, k ∈ N 0 .We then define the fully-discrete TRM to be sequence of values {ρ k i : i ∈ I, k ∈ N 0 } defined by the recurrence relations: ∈ Ω, where the term q k i is given by In particular, if the road is finite, the quantities ρ k j , j ∈ {0, P + 1} are defined from the functions ρ j as ρ k j = (1/∆t) The next result exposes the conditions under which fully-discrete TRM converges to the solution of the PDE.It relies on the following assumption on the source r and sink s terms.Assumption 2.3.The source and sink term q = r − s of PDE (1) is bounded and such that for any ρ ∈ R, (x, t) → q(x, t, ρ) is Lipschitz-continuous (with a constant independent of u) and for any (x, t) ∈ R × R + , ρ → q(x, t, ρ) is locally Lipschitz-continuous (with a constant independent of (x, t) and bounded K 1 and K 2 as defined in (A1)).Theorem 2.4.Let us assume that the source and sink term q = r − s of PDE (1) satisfy Assumption 2.3 and that the following Courant-Friedrichs-Levy (CFL) condition is fulfilled: where K 1 and K 2 are defined in (A1).Then, the fully discrete TRM converges (in L 1 loc (R × R + )) towards the entropy solution of PDE (1) as ∆x, ∆t → 0, with ∆t/∆x kept fixed and satisfying the CFL condition.
Proof.This result is a direct application of Theorem 1 of [4], which holds since the discretized scheme is monotone (cf.Theorem 2.2) and following from the assumptions on the regularity of the source and sink term q.
Remark 2.5.In the absence of source and sink terms (and with an infinite road), but under the CFL condition (8), the discretized numerical scheme in ( 7) is L ∞ -stable, meaning that for any [17]).Consequently, the discrete numerical scheme preserves nonnegativity and capacity (as long as the initial condition is in Ω).Furthermore, since the scheme is monotone and L ∞ -stable, it is Total Variation Diminishing (TVD) [32,Theorems 15.4 & 15.5].As a consequence, the TRM is in particular shock-capturing.It implicitly incorporates correct jump conditions (without oscillations) near discontinuities when ∆x and ∆t tend to 0. Besides, the convergence of the (fully-discrete) TRM towards the PDE solution is of order at most 1 [33].We retrieved this behavior in numerical experiments exposed in Appendix B.

Kinetic and compartmental description
In this section, we will show that the discretized traffic flow model (4) (on a finite road) is formally kinetic with appropriate assumptions on g, and it can be interpreted in a compartmental context.

Kinetic models
Here we characterize kinetic system models based on the notations used in [20], where more details can be found on this model class.We assume that a kinetic model contains M species denoted by X 1 , . . ., X M , and the species vector is X = [X 1 . . .X M ] T .The basic building blocks of kinetic systems are elementary reaction steps of the form where for j = 1, . . ., R, C j = y T j X and C ′ j = y ′ j T X are the complexes and y j , y ′ j ∈ N M 0 are integer coefficients.The transformation shown in (9) means that during an elementary reaction step between the reactant complex C j and product complex C ′ j , [y j ] i items (molecules) of species X i are consumed, and [y ′ j ] i items of X i are produced for i = 1, . . ., M .The reaction ( 9) is an input (resp.output) reaction of species the state vector corresponding to X for any t ≥ 0 (in a chemical context, χ is the vector of concentrations of the species in X ).Then the ODEs describing the evolution of χ in the kinetic system containing the reactions ( 9) are given by where is the rate function corresponding to reaction step i, and determines the velocity of the transformation.For the rate functions, we assume the following for i = 1, . . ., R: (A4) K i (•, t) is locally Lipschitz continuous, and K i (χ, •) is piecewise locally Lipschitz continuous with a finite number of discontinuities, (A5) for j = 1, . . ., m, K i (•, t) depends on χ j if and only if y j > 0, (A6) K i ≥ 0, and for j = 1, . . ., m, K i (χ, t) = 0 whenever χ j = 0 and y j > 0, (A7) There exist continuous nonnegative and strictly monotone (w.r.t. a partial order on R for all t ≥ 0 These properties ensure the local existence and uniqueness of the solutions as well as the invariance of the nonnegative orthant for the dynamics in (10).From now on, a reaction from complex C i to complex C ′ i with rate function K i will be denoted as We will also suppress the t argument in the rate function if it does not explicitly depend on time.A set of nonlinear ODEs given as χ = f (χ) is called kinetic if it can be written in the form (10) with appropriate rate functions K i .We remark that the representation (10) of a kinetic ODE is generally non-unique even if the rate functions are a priori fixed [1].

Kinetic representation of the traffic flow model
Using the above notions, we can give a kinetic interpretation for the TRM, where the species and the reaction steps have a clear physical meaning.For this purpose, we consider the TRM (4) on a finite finite road (with P cells).Let us model each road cell i ∈ I = {1, . . ., P } as a compartment containing two homogeneously distributed species: units of occupied space N i and units of free space S i .The flow of vehicles between the consecutive cells (compartments) i − 1 and i (for i ∈ {2, . . ., P }) is then represented by chemical reactions converting occupied space in one cell into occupied space on the next one, as follows is seen as a sequence of compartments (Bottom row) containing "molecules" of free space S and occupied space N (represented as circles).The flow of vehicles along the road is then modeled by chemical reactions between the compartments (written in red).
where K 1−1,i is the corresponding reaction rate.This model is represented in Figure 2. Note that, in this representation, units of space can flow from cell i − 1 to cell i provided that a unit of occupied space is present in cell i − 1 and that a unit of free space is present in cell i.By recalling condition (A6) in Subsection 3.1, it is easy to see that the model in Eq. ( 11)naturally enforces that there can be no flow from cell i − 1 to cell i if cell i − 1 is empty (i.e.there is no unit of occupied space N i−1 ) or if cell i is full (i.e.there is no unit of free space S i ).
Regarding the on and off-ramps in a given cell i ∈ I, they are respectively described by the following reaction steps where K on,i and K off,i are the reaction rates.These transformations express that during an elementary step of on-ramp (resp.off-ramp) reaction, a unit of free space (resp.occupied space) is consumed, and turned into a unit of occupied space.As for the boundary compartments j ∈ {1, P }, they are characterized by similar reaction steps, namely where the reaction rates K in and K out are set to reflect the boundary conditions.Let n = (n i : i ∈ I) and s = (s i : i ∈ I) denote the concentrations of occupied and free space in the compartment i ∈ I. Within our analogy, these quantities can respectively be seen as the equivalent of the density of vehicles and the density of free space in the cell.Then the kinetic differential system (10) corresponding to the kinetic model with reaction steps ( 11)-( 12) is the system of ODEs defined by where we adopt the notation K 0,1 (n 0 , n 1 ) ≡ K in (n, t) and K P,P +1 (n P , n P +1 ) ≡ K off (n, t).Note that ṅi + ṡi = 0. Hence, c i = n i + s i is a first integral (conserved quantity) of the dynamics which can be interpreted as the maximal vehicle density in cell i.It is also clear that the conservation of c i guarantees the boundedness of the solutions of ( 14)- (15).Substituting then s i = c i − n i into (14) (and assuming that 0 This last equation is equivalent to the TRM equation ( 4) when taking n i = ρ i to be the vehicle density of the i-th cell, c i = ρ max to be the maximal cell capacity, and defining the reaction rates as and for i ∈ I, K on,i (ρ, t) = R i (ρ, t)/∆x, and K off,i (ρ, t) = S i (ρ, t)/∆x.A notable case is when we consider TRM with the decomposition g = g MAK introduced in Example 2.Then, the resulting reaction rates are exactly those obtained when considering the kinetic system under the so-called mass action kinetics [6] assumption.
More generally, the kinetic interpretation of the TRM provides a new insight into the definition of fundamental diagrams, as modeled by the function ρ → f (ρ) linking the flux to the density of vehicles.Indeed, the decomposition g used to define f , which involves both the density of vehicles ρ and its dual the free space density ν, has a clear physical and kinetic interpretation: it models how free space is turned into occupied space by the flow of vehicles.
It is clearly visible from Eq. ( 16) that there is a simple proportional relationship between the kinetic reaction rates and the speed-density relationship of the fundamental diagrams.This gives a transparent and physically meaningful mapping between traffic flows and kinetic models.Therefore, two main conclusions can be drawn.First, the TRM offers a complementary microscopic point-of-view on traffic where the particles being modeled are not vehicles but units of free and occupied space, and which has a clear and straightforward link to the quantities modeled by the macroscopic traffic model.Second, the TRM allows to propose and to interpret fundamental diagram relations based on a microscopic model for the transformations of free space and occupied space.

Application: Dynamical analysis of the ring topology
We conclude thus section with an example of how the kinetic interpretation of the TRM can be leveraged to deduce some properties of the model using kinetic system theory.In particular, we consider the special case when the TRM with P compartments has a ring topology (i.e.periodic boundary conditions) without on and off-ramps, and study the persistence of the dynamics and then the stability of the equilibria.

Persistence of the dynamics using kinetic theory
Following the kinetic interpretation made in Section 3, the TRM can be analyzed using the theory of kinetic systems.Indeed, the TRM is equivalent to the kinetic system defined by the reactions ( 11) and (13).Note in particular that the boundary reactions ( 13) can be merged into a single reaction linking the first and last cells, and given by Since the reaction rates K i,i+1 are given by ( 16) and do not depend explicitly on time, and assuming that they are real analytic functions, we can use the results of [2] on the persistence of kinetic systems.A nonnegative kinetic system is called persistent if no trajectory that starts in the interior of the orthant R P + has an ω-limit point on the boundary of R P + .A non-empty set of species Σ ⊂ X is called a siphon if each input reaction associated to Σ is also an output reaction associated to Σ.A siphon is minimal if it does not contain (strictly) any other siphons.Naturally, the union of siphons is a siphon.The sufficient conditions for the persistence of the dynamics are the following [2]: (1) there exists a positive linear conserved quantity (i.e., a first integral) for the dynamics, and (2) each siphon contains a subset of species, the state variables of which define a nonnegative linear first integral for the dynamics.It is straightforward to see that Condition (1) is fulfilled, since I 0 = P i=1 s i +n i is a conserved quantity for the system.For Condition (2), it can be checked that the only minimal siphons in the ring topology are: Σ N = {N 1 , . . ., N P }, Σ S = {S 1 , . . ., S P }, and Σ i = {N i , S i } for i = 1, . . ., P .Since I N = P i=1 n i , I S = P i=1 s i , and I i = n i + s i for i = 1, . . ., P are also linear first integrals containing the state variables of Σ N , Σ S , and Σ i respectively, Condition (2) is fulfilled.Therefore, the dynamics of the ring topology are persistent for any numerical flux F for which the corresponding reaction rates K i,i+1 satisfy the mild conditions described above.Note that in order to prove the persistence of the dynamics for a wide class of flux functions just by exploiting the ring structure, we needed both sets of state variables representing the density of vehicles as well as the density of free space.

Structure of equilibrium and stability
We first analyze the possible equilibrium points of the system.The proof of the next result can be found in Appendix F.2. Proposition 3.1.For the system of ODEs (4) resulting from the mass action kinetic TRM with a ring topology, the only possible equilibrium point (ρ * 1 , . . ., ρ * P ) is the one satisfying for any Hence, as one could have expected, the only possible equilibrium point for the system is the one where the vehicles distribute themselves uniformly across the (circular) road.Clearly there exist infinitely many equilibria for the system, but there is one unique positive equilibrium within each equivalence class H c = {ρ ∈ Ω | P i=1 ρ i = c} defined by the conservation of vehicles.The next result elaborates on the stability of this equilibrium.Proposition 3.2.The equilibrium ρ is stable on Ω.
Proof.We use the entropy-like Lyapunov function candidate well-known from the theory of reaction networks [19]: Then, which gives then we use the inequality e a (b − a) ≤ e b − e a to give an upper bound such that Note that the entropy-like Lyapunov function V does not depend on the model parameters, and V < 0 outside the equilibria .We also remark that the theory of kinetic and compartmental systems can still be used when the spatial discretization is non-uniform along the ring.In such a case, the analytical computation of the equilibrium is not straightforward.However, we know that the dynamics is persistent, and for any c > 0, there exists a unique strictly positive equlibrium in H c which is stable with the Lyapunov function defined in Eq. ( 17), and asymptotically stable within H c .
It is important to add that the special case studied here can be greatly generalized using recent results.Persistence and stability can be similarly proved not only for the MAK decomposition used here, but for very general reaction rates and thus for a really wide set of fundamental diagrams for arbitrary strongly connected networks [47].Moreover, we can handle explicitly time-varying transition rates (possibly resulting from changing conditions and/or external control) and determine a whole family of logarithmic Lyapunov functions in that case [54].

Analysis of TRM: Comparison with the CTM
In this section, we formally show that the TRM is equivalent to the Cell Transmission Model (CTM) [14] under some specific parametrization of the input capacity (i.e.cell supply function).Further analysis of the TRM is proposed in Appendices C and D where we show how the TRM can be interpreted as a REA algorithm and we analyze how it approximates shock waves (modified equations analysis).
In the CTM, the road is discretized into homogeneous cells of size ∆x (indexed in the ascending order in the direction of circulation of the road).The number of vehicles inside a cell i is represented by a time-dependent function η i .Consider some time step ∆t.The evolution of the number of vehicles η i from a time t to a time t + ∆t is described by the recurrence relation where the quantity y i (t) represents the number of vehicles entering the cell i (from the cell i−1) between t and t+ ∆t.These flows of vehicles are then assumed to be given, for any time t, by the relation where N i (t) denotes the maximal number of vehicles the cell i can hold at time t, and Q i (t) is the so-called input capacity of the cell i, which describes the maximal number of vehicles that could flow into the cell.Note in particular that for any i and any t, η i (t) ∈ [0, N i (t)].
Consider now the case where the maximal number of vehicles a cell can contain is constant, i.e. there exists some constant N such that for any i and any t, N i (t) = N .Let us then take the input capacities Q i defined by Then, assuming that the CFL condition ( 8) is satisfied, we have y i (t) = Q i (t) for any i and any t.Indeed, using the properties of g, we have on the one hand, and similarly on the other hand, Using then the fact that the CFL condition (8) implies that K 1 ∆t/∆x ≤ 1 and K 2 ∆t/∆x ≤ 1, we retrieve that Note then that the density of vehicles ρ i (t) inside the i-th cell is obtained by dividing the number of vehicles η i (t) by the cell size ∆x.Similarly, the maximal density of vehicles is ρ max = N/∆x.Hence, by dividing ( 18) by ∆x, we obtain the relation: where for any i, Hence, (21) coincide exactly with the recurrence relation defining the TRM (7).We then conclude that the CTM, when defined with the input capacities (20), fits into the TRM framework proposed in this work.This particular choice of input capacities simplifies greatly the expression (19) defining the flow of vehicles: the minimum disappears and we are left with the numerical fluxes of the TRM.In conclusion, the TRM elegantly incorporates in a compact (analytic) form the restrictions on the flux of vehicles allowed to travel from cell to the other.The resulting smoothness opens the door for the application of results on polynomial systems of ODEs for the analysis of the numerical properties of the TRM.

Possible extensions of the TRM
In this section we present two straightforward extensions of the TRM that allow to model traffic in more complex situations that those covered by the LWR model (1) on a unidirectional road.

Accounting for changing driving conditions
The TRM can be adapted to model roads with non-homogeneous driving conditions, may they be due to changes in the maximal vehicle density admissible along the road (due for instance to varying number of lanes) or more generally to external factors affecting the optimal flow of vehicles.Indeed, consider once again a road discretized into cells i ∈ I, and let ρ max denote the capacity of the i-th cell.We call Extended TRM the system of ODEs defined by where for each i ∈ I, the numerical flux F i is given by , the function g is defined in the same way as in Section 2, and the function t → C i (t) takes values in (0, 1] can be seen as local capacity drop factor which scales down the "ideal" flow of vehicles F i (ρ i−1 , ρ i ) between cells i − 1 and i, i.e. the flow of vehicles one should expect between these cells given their density states.By identification with the usual TRM, this factor can essentially be seen as a time-varying normalized free flow speed v max at the interface between cells i−1 and i, or within the kinetic compartmental interpretation, as a time-varying normalized reaction rate coefficient between compartments i−1 and i.On the other hand, allowing the maximal capacity ρ (i) max to change across compartments allows for instance to account for changes in the number of lanes or for non-uniform discretizations of the road.Following the same arguments as the ones used in Theorem 2.2, it is clear that the Extended TRM will yield solutions that preserve non-negativity and capacity.Besides, the extended TRM has been shown to accurately describe realworld traffic datasets in applications related to traffic state estimation [39] and short-term prediction [40].This is manly due to the local capacity drop factors C i which allow to locally change the conditions at which traffic flows on the roads, hence allowing the TRM (which inherits from the seemingly simplistic LWR model) to recreate complex traffic patterns observed in the data.To illustrate this last property of the Extended TRM, we present in Figure 3 the results of a simulation study where traffic along a unidirectional road with a traffic light is modeled.To do so, we start by discretizing the road into compartments and consider that at the position x = 0 where the traffic light is located, the corresponding capacity drop coefficient C(t) oscillates between two values: 1 when the traffic light is green, and 0 when it is red (cf.Figures 3a  and 3b).Following the kinetic interpretation of the TRM, this choice will effectively stop the flow of vehicles at x = 0 when the traffic light is red.The mass action kinetic decomposition is chosen, which in particular yields a polynomial system of ODEs which can be solved at arbitrary time steps with high accuracy using standard libraries.
The Extended TRM is run on a road of 5km discretized into compartments of 5m.The free flow speed of vehicles is set to 30km/h and the simulation is run on a time horizon of 10min, with light switches happening approximately every 2 minutes.The results presented in Figures 3c to 3e seem to show that the model manages to recreate both the wave propagations and the oscillatory nature the density variations due to the traffic light changes.

Accounting for road networks
The TRM can also be directly adapted to model traffic on networks.Indeed, following the kinetic and compartmental interpretation given in Section 3, traffic on a network can be assimilated to a compartmental chemical reaction network for which each compartment is now allowed to have more than two other neighboring compartments.Starting from a road network, modeled as a directed graph, each directed edge is discretized into compartments, and each intersection is replaced by an additional compartment.An example of such discretization is given in Figure 4, and show in particular how the TRM can be extended to handle complex intersections.The Network TRM then takes the form of a system of ODEs defined as where for each compartment i ∈ I, N in (i) denotes the set of compartments j that point to i (i.e.such that vehicles can move from j to i) and N out (i) denotes the set of compartments k that i points to (i.e.such that vehicles can move from i to k).Once again, following the same arguments as the ones used in Theorem 2.2, it is clear that the Network TRM will also yield solutions that preserve non-negativity and capacity.This qualifies the Network TRM to act as a generic link and node model for traffic networks according to the classification of [50].Following the categories in [50], the Network TRM extension is an unsignalized node model where supply and capacity constraints are handled in a balanced and continuous way.This latter property roots in fact in the balanced interaction of chemical species changing the concentrations in the compartments.In terms of supply constraint interaction rules [50], the Network TRM generates them implicitly: free space concentrations are transformed into occupied ones dynamically, smoothly.
In terms of the solution algorithm, the set of nonlinear ODEs guarantees the existence of a unique solution to the density propagation accross all link and node compartments (via Lipschitz continuity [29]).Node and link models are described in a unified way with TRM.Links and nodes rely on the reaction rate governed exchange of material (vehicle) in-and outflows.Merging and diverging is handled by in-and outgoing flows (see, eq. ( 22)).Besides, the Network TRM can readily be further extended in order to account for varying road conditions using the same approach as the one outlined for the Extended TRM.In practice, the link between the TRM and chemical reaction networks can be leveraged to deduce persistence and Lyapunov stability results that generalize even to time-varying networks, the properties of which are outlined in Section 3.3 [47,48,53].

Conclusions
Traffic Reaction Model (TRM), a family of Finite Volume Methods for segmenting certain hyperbolic Partial Differential Equations arising in traffic flow modelling, has been presented in the paper.First, the proposed numerical approximation scheme is consistent, monotone, nonnegative, capacitated, and conservative, even under non-zero sink and source terms.Second, the resulting system of nonlinear ODE models enables us to formally view traffic flow models as chemical reaction networks.In the kinetic model formulation, two dual state variables appear representing traffic density and the density of free space, respectively.These dual densities enable a proper factorization of the nonlinearity in the original PDE, as done by the function g(), and also allow us to define multiple meaningful numerical fluxes.
We have compared the TRM to other known popular discretization schemes, and clarified their relations.Here we emphasize that TRM is equivalent to the CTM for a particular choice of input capacities.Numerical results show that the convergence properties of the proposed model are comparable to other popular schemes.The main benefits of the proposed modeling approach are the possibility to represent the system model in the form of smooth ODEs which is often required to apply advanced nonlinear control techniques, and to use the tools of compartmental systems and reaction network theory in the dynamical analysis of traffic models.The latter one was illustrated through the example of the ring topology, where the robust persistence and stability of the dynamics were shown using a Petri net representation (containing also both sets of state variables), and a parameter-free logarithmic Lyapunov function, respectively.Finally, we have presented the possible extensions of the model firstly in the form of a time-varying nonlinear system, with basically the same physical interpretation, to describe capacity changes due to e.g., traffic lights, and secondly, as a straightforward network formulation which can be used to model an arbitrary road (directed graph) structure, thus allowing to better model real traffic data [39,40].Further work will be focused on structured stability analysis, state estimation, and robust control design.An interesting research direction can be to further analyze the dual variable triggered parametrization of f (ρ) within the framework of variational theory in order to see how the Hamilton-Jacobi representation (see, [31]) is changed.

B.1 Shock wave
Consider the Riemann-problem with the initial condition ρ(x, 0) = 0.1ρ max if x < L/2 0.8ρ max otherwise, which induces a solution consisting of a shock wave.Figure 5 shows the discretization errors of various monotone discretization schemes (both semi and fully discrete): the mass action kinetic (MAK) TRM ( 5), the Godunov (Gdnv) scheme (6), and the modified Lax-Friedrichs (LxF) scheme 4 .The errors decrease approximately as O(P −1 ) for all (semi and fully discrete) schemes, and the semi-discrete formulation of a given scheme systematically over-performs its fully discrete formulation.This is expected as the fully-discrete formulation add additional error due to the time discretization.If we compare the semi or fully discrete schemes with one another, the Godunov scheme systematically overperforms the other two, which display similar errors.

B.2 Rarefaction wave
Consider the Riemann-problem with the initial condition ρ(x, 0) = 0.8ρ max if x < L/2 0.1ρ max otherwise, which induces this time a solution consisting of a rarefaction wave.The discretization errors associated with this choice of initial condition are given in Figure 6.Once again, the semi-discrete formulations overperform compare to the fully-discrete ones, and the Godunov scheme overperforms compared to the other two schemes, that display similar errors.But this time, the errors at a slightly lower order of O(P −3/4 ).This difference could be explained by the fact that (spatially) piecewise constant reconstructions are used here to approximate a continuous rarefaction fan.

Appendix C Interpretation of the TRM as a REA algorithm
Let us consider the fully discrete TRM with no source terms.Then, if ρ k i denotes the density of vehicles and ν k i = ρ max − ρ k i denotes the density of empty space in the i-th cell at time t k (for any i ∈ Z, k ∈ N), we have Let us then define the quantities α l (i) and α r (i) by where by convention we take α l (i) = 0 (resp.α r (i Note in particular that since g is non-decreasing with respect to each of its argument, we have α l (i) ≥ 0 and α r (i) ≤ 0. We can then write the TRM recurrence relation as where Hence, the update of the vehicle density within the i-th cell can be written as the average between two quantities: one depending only on the vehicle density of the cells i − 1 and i, and one depending only on the density of empty space of the cells i and i + 1.It turns out that these quantities can be linked to the solutions of some specific advection problems.Indeed, consider the i-th cell, centered at x − i, its left border x i−1/2 , and right border x i+1/2 (cf. Figure 7).Let us then introduce the functions η l and η r respectively centered around the left and right border of the cell, and given by .
(i − 1)-th cell i-th cell (i + 1)-th cell Left half Right half Consider the advection problem defined by Then, the average over the left half of the i-th cell of the solution at time t k+1 = t k + ∆t of this advection problem is given by the quantity Similarly, the quantity I − α r (i), ν k i , ν k i+1 corresponds exactly to the cell average over the right half of the i-th cell of the solution at time t k+1 of the advection problem with coefficient α r (i) ≤ 0 and (piecewise constant) initial condition η r .Hence, the TRM can be seen as a particular instance of Reconstruct-Evolve-Average (REA) algorithm [33].Indeed, assuming the density values at time t k are known, the values at time t k+1 can be seen as resulting from the following three-step process: Step R Reconstruct a density function ρ(t k , •) which is constant over each cell from the current cell estimates ρ k i and define the free space as Step E Evolve, at each inter-cell boundary x i+1/2 , the system of (decoupled) advection equations given by ρ t + α l (i + 1)ρ x = 0 ν t + α r (i)ν x = 0 from the initial conditions (ρ(t k , •), ν(t k , •)), and for a time ∆t.
Step A Compute the cell estimate ρ k+1 i of the i-th cell by averaging the value of the cell average of ρ over the left half of this cell, and the value of cell average of (ρ max − ν) over its right half.REA algorithms include the Godunov scheme [33], which can be described as an implementation of the update rule (24) where the advection problems are replaced by coupled system of conservation laws This adjoint system can actually be decoupled and then yields that ρ satisfies our initial PDE (without no sink and source terms) and the relation ν = ρ max − ρ.Hence, we actually proved that solving the resolution of the adjoint system as implemented in the Godunov scheme is equivalent to that of advection problems with coefficients as defined in ( 25) and Other choices of decompositions g(ρ, ν) result in different advection coefficients and consequently may trigger different numerical approximation errors.A systematic and comparative analysis of different g(ρ, ν) functions can be obtained by evolving the time-space bounds to their respective advection coefficients and compare them.
Let us consider the fully discrete TRM defined from the mass action kinetic decomposition, and assume that there exists some constant δ = (∆t/∆x)v max ∈ (0, 1/2).Then, the recurrence relation defining this scheme can be written in the following form (for i ∈ Z, k ∈ N 0 ) Replacing then, in the left-hand side of ( 27), the estimates ρ n j by the corresponding evaluations ρ(t n , x j ) of the PDE solution ρ, yields the so-called local truncation error (LTE) L ρ (t, x), given for any t ≥ 0, x ∈ R by: The LTE is the difference between the true solution, one time step ahead, and the approximation of this quantity using the recurrence relation of the scheme.Hence, it provides an indication of how well the recurrence relation defining a numerical scheme approximates the variations (in time) of the true solution.Indeed, if L ρ (t, x) = 0 for any t ≥ 0, x ∈ R, and assuming that ρ 0 i = ρ(0, x i ) for any i ∈ Z, then the approximate solution TRM recurrence (27) will coincide with the true solution at any time t k and location x i .Building on this idea, the modified equation approach [32] proposes to find smooth functions ρ, for which estimates of the magnitude of LTE (and hence of the ability of the scheme to approximate ρ) are available.This done by expanding the LTE L ρ in (28) using a Taylor expansion of ρ, and truncating the result at a given order n.Taking then ρ to be the solution of the PDE defined by the remaining terms of the expansion yields that the corresponding LTE L ρ must be of order O(∆x n ).For instance, when truncating at an order 1, we get that ρ = ρ is the solution of the conservation law (26) (with no source terms), and therefore that the LTE of the TRM for ρ satisfies L ρ (t, x) = O(∆x) (which is in accordance with the fact that monotone schemes for conservation laws converge at most as fast as O(∆x) [33]).Considering a truncation of order 2 yields the following PDE to define ρ, for which we know that the corresponding LTE L ρ satisfies L ρ(t, x) = O(∆x 2 ).It is then natural to consider, for simple cases, the solutions of the so-called modified equation (of order 2) (29), in order to get a better understanding of the numerical properties of the TRM scheme, since the scheme yields better approximations of the solutions of this PDE (L ρ(t, x) = O(∆x 2 )) than of the solutions of the conservation law (L ρ (t, x) = O(∆x)).
One such case is the study of the approximation of shocks by the numerical scheme.When the initial condition of the conservation law contains a discontinuity such that the left value ρ l is less than the right value ρ r , then the solution consist of a shock traveling with a speed v s = v max (1 − (ρ l + ρ r )/ρ max ).A similar behavior is observed for some solutions ρ of (29), as stated in the next proposition (proof in Appendix F.3). Proposition D.1.Let u denote the sigmoid function, which is defined by u where x 0 ∈ R is an arbitrary constant and Then, ρ is the only traveling wave solution of (29) satisfying lim Note that the traveling wave solution described in Proposition D.1 travels at the same speed as the shocks of PDE (26).
The wave in question consists in a smooth approximation of a shock from ρ l to ρ r using a sigmoid function.As for the arbitrary constant x 0 , it sets the position, at time t = 0, of the inflection point of the traveling wave.The parameter σ, on the other hand, sets how sharp the transition between ρ l and ρ r is: indeed, for any ϵ ∈ (0, 1), the errors |ρ l − ρ| and |ρ r − ρ| are smaller than ϵ(ρ r − ρ l ) at a distance (from the inflection point) greater than σ| log(ϵ/(1 − ϵ))|.For instance, errors smaller than 0.01(ρ r − ρ l ) are achieved 5σ away from the inflection point.Consequently, the quantity σ can be used to get an idea of the extent with which the TRM smoothes out the discontinuities.In particular, since σ is proportional to ∆x, the expression of σ given in Proposition D.1 can be used to directly estimate the number of cells needed to transition from the left state of a shock to the right state.Also, since σ is inversely proportional to the shock size ρ r − ρ l , we expect this transition to be longer when approximating small shocks than when approximating bigger shocks.Note that similar conclusions can be drawn for the Lax-Friedrichs (LxF) scheme.Indeed, using the same approach, one can show that the modified equation of order 2 of the LxF scheme takes the form which differs from the TRM only by the factor in front of ρxx .Hence, the same reasoning as the one used in Proposition D.1 can be applied, and we can deduce that the traveling wave solution of (31) can once again be written as in (30), but with a parameter σ = σ LxF now equal to In particular, since δ ∈ (0, 1/2) we have that σ TRM < σ LxF , and therefore, we can expect the mass action kinetic TRM to smooth less the shocks than the LxF scheme.Finally, we present in Figure 8 a comparison between approximations of a solution of PDE (26) using the mass action kinetic TRM and the LxF scheme.We also plot the solutions of the corresponding modified equations of order 2. As expected, the mass action kinetic TRM and LxF approximate solutions are close to their respective solutions to the modified equation.Since we consider ρ l = 0.2 and ρ r = ρ m , we expect the mass action kinetic TRM solution to be close to the left and right state if we are at a distance of more than 5σ TRM ≈ 3.5∆x from the inflection point.As for the LxF scheme, we expect this distance to be 5σ LxF ≈ 9∆x.The figure confirms these estimates.(26) with no source and sink term, and initial condition equal to ρ l = 0.1 if x < 0 and ρ r = 0.9 if x > 0 (ρ max = 1, v max = 1).The solutions are computed at time t = 1, and with parameters ∆x = 10 −3 and δ = (∆x/∆t)v max = 0.4.The LxF solution is denoted by LxF, the solution of modified equation of LxF is denoted by Meq LxF, the mass action kinetic TRM solution is denoted by MAK TRM and the solution of the modified equation of the mass action kinetic TRM is denoted by Meq MAK TRM.

Appendix E Solutions of advection problems
Let w l , w r ∈ R, and define w 0 by w 0 (x) = w l if x < 0 w r if x > 0 For α ∈ R, consider the following advection problem: The solution of ( 32) is given by w(t, x) = w 0 (x − αt), t ∈ R + , x ∈ R,  Appendix F Additional proofs F.1 Proof of Theorem 2.2 Proof.The numerical flow F is consistent since by definition and following (2), we have F (u, u) = g(u, ρ max − u) = f (u) for any u ∈ Ω.It is monotone due to the continuity assumption on f , consistency and the non-decreasing properties of g.

Figure 1 :
Figure 1: Comparison between quadratic and trapezoidal fundamental diagrams.The maximal density ρ max and the free flow speed v max are taken equal to 1.The black dotted lines are placed at the critical density values ρ 1 = 0.25 and ρ 2 = 0.6.

Figure 2 :
Figure 2: Representation of the Traffic Reaction Model interpretation of traffic flow.The discretized road (Top row)is seen as a sequence of compartments (Bottom row) containing "molecules" of free space S and occupied space N (represented as circles).The flow of vehicles along the road is then modeled by chemical reactions between the compartments (written in red).
Representation of the road and its discretization (b) Time evolution of the capacity drop factor C(t) (c) Space-time plot of the TRM density simulation (d) Time evolution of the density at x = ±52.5 (e) Time evolution of the density at x = ±1002.5

Figure 3 :
Figure 3: Extended TRM simulation of the normalized density on a road with a traffic light.The red and green lines indicate when the capacity drop factor is set to 0 (red lines) or to 1 (green lines).

Figure 5 :
Figure 5: L 1 -norm and L ∞ -norm errors of the shock wave discretization for various numbers of segments P .

Figure 6 :
Figure 6: L 1 -norm and L ∞ -norm errors of the rarefaction wave discretization for various numbers of segments P .

Figure 7 :
Figure 7: Left and right halves of the cells.

Figure 8 :
Figure8: Comparison of the Lax-Friedrichs (LxF) scheme and the mass action kinetic TRM.The true solution (in black) corresponds to PDE(26) with no source and sink term, and initial condition equal to ρ l = 0.1 if x < 0 and ρ r = 0.9 if x > 0 (ρ max = 1, v max = 1).The solutions are computed at time t = 1, and with parameters ∆x = 10 −3 and δ = (∆x/∆t)v max = 0.4.The LxF solution is denoted by LxF, the solution of modified equation of LxF is denoted by Meq LxF, the mass action kinetic TRM solution is denoted by MAK TRM and the solution of the modified equation of the mass action kinetic TRM is denoted by Meq MAK TRM.

2 IFigure 9 :( 33 ) 2 0−∆x/ 2 w
Figure 9: Cell averages at the discontinuity.The discontinuity is represented in red for the initial time and in blue after ∆t.