Stability analysis of a hyperbolic stochastic Galerkin formulation for the Aw-Rascle-Zhang model with relaxation

We investigate the propagation of uncertainties in the Aw-Rascle-Zhang model, which belongs to a class of second order traffic flow models described by a system of nonlinear hyperbolic equations. The stochastic quantities are expanded in terms of wavelet-based series expansions. Then, they are projected to obtain a deterministic system for the coefficients in the truncated series. Stochastic Galerkin formulations are presented in conservative form and for smooth solutions also in the corresponding non-conservative form. This allows to obtain stabilization results, when the system is relaxed to a first-order model. Computational tests illustrate the theoretical results.


Introduction
Nowadays traffic models have become an indispensable tool in the urban and extraurban management of vehicular traffic. Understanding and developing an optimal transport network, with efficient movement of traffic and minimal traffic congestions, will have a great socioeconomical impact on the society, in particular in pandemics situations.
Besides guaranteeing optimal transport in the presence of pandemic situations, there is a second major aspect, where our work on traffic flow modelling may contribute. It is clear that in a pandemic situation the spreading of possible infections correlates with the number of contacts as e.g. modelled in SIR dynamics [30,39]. Traffic flow provides valuable information on possible contacts and on possible points of high population density in urban and extraurban areas. The prediction of the flow into and from those areas can help to calibrate the transmission coefficients in typical SIR models for disease propagation. Here, however, deterministic predictions are of little to no use in an a priori assessment of possible critical points of high traffic density. Therefore, it is mandatory to expand the current theory on macroscopic deterministic traffic flow models towards realistic but uncertain models. The current paper precisely tackles this point.
A vast amount of literature about vehicular traffic modeling has flourished in the last decades. Nevertheless, there are still several limitations for obtaining trustful traffic forecasts. This is possibly due to the fact that the evolution of traffic is described by highly nonlinear dynamics that is also exposed to the presence of various sources and types of uncertainties [45,26,25,51]. For example, the uncertainty may stem from real data affected by errors in the measurements or the reaction time of drivers. A pandemic scenario adds additional uncertainties, but needs reliable estimates. In particular, in view of the discussion of possible measures to reduce traffic and accumulation in certain areas, the reliable and quantifiable prediction is of high importance. The approach presented in this paper allows to quantify the complete statistics of the uncertain solution and hence it also allows to compute e.g. rare events. Quantifying the propagation of uncertainty in nonlinear models is therefore of interest and the purpose of this paper.
Uncertainty quantification in the sense used here is concerned with the propagation of input uncertainty through traffic models. Several approaches are presented in the literature and can be classified in non-intrusive and intrusive methods. The main idea underlying the former approach is to solve the model for fixed number of samples using deterministic numerical algorithms. Then, the statistics of the quantities of interest are determined by numerical quadrature. Typical examples are Monte-Carlo and stochastic collocation methods [37].
In contrast, we consider the intrusive stochastic Galerkin method. Here, stochastic processes are represented as piecewise orthogonal functions, for instance Legendre polynomials or multiwavelets. These representations are known as generalized polynomial chaos (gPC) expansions [4,16,17,54,56]. Expansions of the stochastic input are substituted into the governing equations and a Galerkin projection is used to obtain deterministic evolution equations for the coefficients of the series expansions.
Results for nonlinear hyperbolic systems are only partial, since desired properties like hyperbolicity are not necessarily transferred to the intrusive formulation [8,27]. A problem is posed by the fact that the deterministic Jacobian of the projected system differs from the random Jacobian of the original system. We refer the interested reader to [15] for examples of the Euler as well as shallow water equations. Furthermore, it is remarked in [42,Sec. 5] that simulations for Euler equations may break down for high Mach numbers unless auxiliary variables and wavelet-based expansions are used.
Still, stochastic Galerkin methods applied to hyperbolic equations is an active field of research. Those can be successfully applied to scalar conservation laws, since the resulting Jacobian is symmetric. In the scalar case, well-balanced schemes have been developed [29] and a maximum-principle can be ensured [32].
Furthermore, entropy-entropy flux pairs and hence hyperbolicity can be transferred to a stochastic Galerkin formulation by introducing auxiliary variables [8], which require expensive variable transforms. Although there are many attempts to make the transform more efficient and stable [33,34], the computational cost remain a drawback of this approach. To this end, an expansion in Roe variables has been proposed [42]. Since it exploits quadratic relationships, the necessary transforms are numerically cheap and stable. These auxiliary variables enable also a hyperbolic stochastic Galerkin formulation for isothermal Euler equations for arbitrary gPC expansions. Moreover, it has been observed that the shallow water equations allow for a hyperbolic stochastic Galerkin formulation which neither requires auxiliary variables nor any transform [6].
Additional results are available for certain wavelet-based gPC expansions, including the Wiener-Haar basis and piecewise linear multiwavelets [38,42]. These wavelet expansions are motivated by a robust expansion for solutions that depend on the stochastic input in a nonsmooth way and are used for stochastic multiresolution as well as adaptivity in the stochastic space [1,31,52].
In this paper, we consider hyperbolic systems used in vehicular traffic modeling, namely second order macroscopic models [2,57]. The main feature is that they take into account the non-equilibria states, assuming that accelerations are not instantaneous. They are able to recover typical traffic phenomena as generating capacity drop, hysteresis, relaxation, platoon diffusion, or spontaneous congestions like stop-and-go waves [18,11,24].
The first results in this direction were proposed by Payne and Whitham [53] taking into account that the speed of each car does not change instantaneously. However, their model has the drawback that the driver's decision is influenced by the road conditions behind. A second order model is due to Aw, Rascle [2] and Zhang [57]. By taking into account the differences between traffic and fluid flows, they designed models to simulate the anisotropic traffic behaviour.
The inhomogeneous Aw-Rascle-Zhang (ARZ) model includes a relaxation term that allows drivers to achieve the equilibrium speed [19]. In the small relaxation limit the ARZ model approaches to the Lighthill-Whitham-Richards (LWR) model [36,46], which can be obtained by means of a Chapman-Enskog-type expansion. Here, the stability and well-posedness of solutions to the hyperbolic ARZ model is governed by the study of the sign of the diffusion coefficient, which requires the so-called sub-characteristic condition [5,28]. The diffusion term vanishes in the zero-relaxation limit and the LWR model is recovered [47,23,24]. This paper analyzes stochastic Galerkin formulations for the Aw-Rascle-Zhang model in conservative and non-conservative form. The non-conservative form allows to state eigenvalues and hence ensures hyperbolicity. Furthermore, the stability of the system is investigated if it is relaxed to a first-order model. As basic tool we follow the approach in [47,23,24] and study definiteness properties of the corresponding diffusion coefficient by using a Chapman-Enskog-type expansion.
Section 2 introduces the deterministic Aw-Rascle-Zhang model in conservative and nonconservative form. Section 3 presents stochastic Galerkin formulations. For a special class of wavelet-based gPC expansions an auxiliary variable that does not cause any computationally expensive transforms is introduced to ensure hyperbolicity. Section 4 is devoted to a stability analysis of the inhomogeneous ARZ model. The theoretical results are derived only for classical smooth solutions with deterministic relaxation. Riemann problems to weak solutions with uncertainties in the relaxation parameter are illustrated numerically in Section 5.

Second order traffic flow models with relaxation
Typical macroscopic traffic flow models describe the density ρ = ρ(t, x) and the mean velocity v = v(t, x) of vehicles at a location x ∈ R and time t > 0. The natural assumption that the total mass is conserved leads to impose that the density ρ satisfies the continuity equation In first-order models the velocity v = v(ρ) is given as a function of the density alone, e.g. the LWR model [36,46]. Second-order models describe the velocity by an additional differential equation. In particular, we consider the inhomogeneous Aw-Rascle-Zhang model [2,19] with relaxation   Here, h(ρ) : R + → R + is called hesitation or traffic pressure [10]. It is a smooth, strictly increasing function of the density. The relaxation term with parameter τ > 0 on the right hand side makes the drivers tend to a given equilibrium velocity V eq (ρ). This is important, since the homogeneous ARZ model without relaxation has no mechanism to move drivers when initially are at rest. By introducing the variable z = ρ v + h(ρ) , the system (2) can be written in conservative form as Here, the velocity v(ρ, z) is a driver dependent property. The conservative formulation (3) is abbreviated as The eigenvalues of the Jacobian are λ 1 (ρ, z) = v(ρ, z) − ρh (ρ) and λ 2 (ρ, z) = v(ρ, z). Hence, the ARZ model is strictly hyperbolic under the assumption ρ > 0. The (local) equilibrium velocity V eq (ρ) satisfies the scalar conservation law ∂ t ρ + ∂ x f eq (ρ) = 0 for f eq (ρ) = ρV eq (ρ) and f eq (ρ) = V eq (ρ) + ρV eq (ρ).
Stability requires that the full system propagates information faster than the local equilibrium, i.e. the sub-characteristic condition leads to a dissipative advection-diffusion equation. For the deterministic ARZ model [57,23], this reads as (DI) In the sequel, we will extend these results to the stochastic case.

Stochastic Galerkin formulation
We extend the hyperbolic balance law (3) to account for uncertainties that arise from random initial conditions. The hesitation function and the equilibrium velocity, however, remain given deterministic functions. Uncertainties are summarized in a random variable ξ, defined on a probability space Ω, F(Ω), P , and propagated by the random system For fixed time and space coordinates we expand the solution in terms of the generalized polynomial chaos (gPC) expansion The piecewise polynomial functions φ k (ξ) form an orthonormal basis with respect to the weighted inner product If the random solution u(t, x, ξ) is known, the gPC modes can be determined by the orthonormal projection u(t, x, ·), φ k (·) . Under mild conditions on the probability measure the truncated expansion (gPC) converges in the sense 12,9]. A challenge occurs, since only the gPC modesû(0, x) corresponding to the initial data are known. To determine them for t > 0, we derive a differential equation, called stochastic Galerkin formulation, that describes their propagation in time and space.

A semi-intrusive approach as introductory example
A naive approach would be to substitute the truncated expansion (gPC) into the random system (6) and then use a Galerkin ansatz to project it onto the space spanned by the basis functions. The resulting system, without relaxation term, reads as and Jacobian Dûf consisting of block matricesf α,β û(t, Here, the Jacobian Dûf (û) consists of the projected entries of the deterministic Jacobian (4). The Jacobian (8), however, has not necessarily real eigenvalues and a full set of eigenvectors.
In the case of the Aw-Rascle-Zhang model, the flux function (7) and its Jacobian (8) are not even directly specifiable, since the deterministic expressions (3) and (4) envolve the terms z 2 /ρ, z /ρ and the possibly nonpolynomial hesitation function h(ρ). Computing numerically the integrals in equation (7) and (8) would lead to an expensive, non-hyperbolic, semi-intrusive scheme.

Intrusive formulation for general gPC expansions
Instead, we follow the approaches [27,6] to handle the terms z 2 /ρ and z /ρ. We introduce the Riemann invariant w := z /ρ in the original ARZ model [2]. While the semi-intrusive approach in Section 3.1 computes the gPC modes w by the orthonormal projection w, φ k , we project the product ρw and determine the modes by the pseudo-spectral Galerkin product Similarly to [7,37,50], we express it bŷ The matrix P(ρ) is strictly positive definite and hence invertible provided that the gPC expansion G K [ρ] > 0 is strictly positive [49,55,15]. The strict positive definiteness of the matrix P(ρ) is assumed throughout this paper. This assumption excludes vacuum states. We have for the inverse terms the pseudo-spectral gPC approximations w = P −1 (ρ)ẑ andẑ * w, i.e.
This yields for general gPC bases a stochastic Galerkin formulation for the homogeneous ARZ model, without relaxation, as whereĥ(ρ) ∈ R K+1 denotes a given gPC formulation of a hesitation function. For example, the linear hesitation function h(ρ) = ρ has the gPC modesĥ(ρ) =ρ. By using the following calculation rules, see e.g. [42,27,14], we obtain the Jacobian of the gPC formulation (10) as where 1 := diag{1, . . . , 1} denotes the identity matrix. The matrices M k and hence the linear operator P : R K+1 → R (K+1)×(K+1) , defined in equation (9), are exactly computable in an offline stage. Therefore, the stochastic Galerkin formulation (10) is intrusive and no numerical quadrature is needed during a simulation. Furthermore, the eigenvalues can be exactly computed. However, eigenvalues are not proven real which motivates the following subsection.

Hyperbolic and intrusive formulation for wavelet-based gPC expansions
Under additional assumptions on the bases functions, hyperbolicity can be guaranteed. We consider basis functions φ k that satisfy the following properties: (A1) The precomputed matrices M and M k commute for all , k = 0, . . . , K.

Stochastic Galerkin formulation for the inhomogeneous ARZ model
The hyperbolic formulation, presented in Subsection 3.3, is directly extendable to a stochastic Galerkin formulation for the inhomogeneous ARZ model. To this end, we assume an arbitrary, but consistent gPC expansion V eq (ρ) of the random equilibrium speed V eq ρ(ξ) , satisfying Then, we introduce a stochastic Galerkin formulation of the relaxation term in the conservative formulation (3) by This auxiliary variable also allows to obtain a stochastic Galerkin formulation for the nonconservative formulation (2). Altogether we have the hyperbolic stochastic Galerkin formulations for the inhomogeneous ARZ model in a conservative form non-conservative form We show in Theorem 3.1 that these two formulations are equivalent for smooth solutions, as it holds in the deterministic case [2]. However, if there is a jump in the solution, the nonconservative form contains the product of the discontinuous matrix-valued function P(v) with the distributional derivative of the termv +ĥ(ρ), which may contain a Dirac mass at the point of the jump. In general, such a product is not well-defined [3, Sec. 1]. Theorem 3.1 ensures that the system is strongly hyperbolic, which means that eigenvalues of the Jacobian Dûf (û), i.e. the characteristic speeds of the hyperbolic system are real. Moreover, the Jacobian Dûf (û) admits a complete set of eigenvectors which implies that classical solutions are well-posed [21].
where D(v) denote the eigenvalues of the matrix P(v). Furthermore, the stochastic Galerkin formulations (N ) and (C) are strongly hyperbolic in the sense that the characteristic speeds are real and the Jacobian Dûf (û) admits a complete set of eigenvectors.
Since the opertor P(ρ) is linear, the homogeneous part of the non-conservative formulation can be rewritten as Here, we have used the equalityρ * (v * ∂ x ) = (ρ * v) * ∂ x , which is satisfied provided that the assumptions (A1) -(A3) hold, but not for general gPC bases, since the Galerkin product is typically not associative [50,7]. Likewise, the relaxation term in the conservative formulation is obtained by multiplying P(ρ), i.e. by applying the Galerkin product to the relaxation term of the non-conservative form. Therefore, the two formulations (C) and (N ) are equivalent. To state the eigenvalues, we rewrite the first equation of the non-conservative form as where we have used the symmetry of the Galerkin product. By subtracting equation (15) from the second equation in the non-conservative form and by using property (A2), i.e. an eigenvalue decomposition with constant, orthonormal eigenvectors V T = V −1 , we obtain for 0 ∈ R 2(K+1) and O ∈ R (K+1)×(K+1) . Due to the sparsity structure in the quasilinear form (16) a complete set of eigenvectors exists and eigenvaluesλ 1 ,λ 2 are obtained.

Stability analysis of the inhomogeneous ARZ model
The parameter τ > 0 determines the relaxation of the velocityv(ρ,ẑ), given by equation (14) as auxiliary variable, towards the gPC modes V eq (ρ) of the equilibrium velocity, which is a function of the density alone. We study in this section small, but positive values of the relaxation paramter τ > 0, when the ARZ model is close to the equilibrium model with Jacobian Dρ ρ * V eq (ρ) = P V eq (ρ) + P(ρ)Dρ V eq (ρ).
We observe from the Jacobian (18) that an eigenvalue decomposition of the equilibrium velocity of the form should be assumed such that all waves of the equilibrium model propagate at the characteristic speeds λ eq (ρ) := D V eq (ρ) + D(ρ)D V eq (ρ) not exceeding the equilibrium velocity. This is identified by the eigenvalues of the matrix P V eq (ρ) . Analogously to the analysis in [5,47,23,24], we use a Chapman-Enskog-type expansion that allows to study the behaviour of first-order perturbations of the equilibrium velocity. This yields a diffusion correction as stated in the following theorem.
Theorem 4.1. Let a gPC expansion with the properties (A1) -(A3), a stochastic Galerkin formulation of a hesitation functionĥ(ρ) and a Galerkin formulation of an equilibrium velocity V eq (ρ) be given. Assume further that the Jacobians can be written as with constant eigenvectors. The first-order correction to the local equilibrium approximation reads Furthermore, it is dissipative if and only if the sub-characteristic condition Proof. We apply a Chapman-Enskog expansion The linearity P(α + β) = P(α) Hence, in the non-conservative formulation we obtain The symmetry of the Galerkin product and the equilibrium model (17) yield which implies the claim Theorem 4.1 gives conditions to properly choose a hesitation function h(ρ) and an equilibrium velocity V eq (ρ). In the deterministic case, various choices have been investiated to model also phantom traffic jams and stop-and-go waves by introducing a negative diffusion coefficient [18,11,24]. Here, we investigate states close to the equilibrium and choose a hesitation function h(ρ) and an equilibrium velocity V eq (ρ) such that sub-characteristic condition is fulfilled. The following corollary extends a widely used class, which includes the Greenshields flux, see e.g. [20,19,48] for the deterministic case, to the derived stochastic Galerkin formulation.

Corollary 4.2. Let an equilibrium velocity and a hesitation function of the form
with strictly positive constants v max , ρ max , γ be given. Under the assumptions of Theorem 3.1 and Theorem 4.1 the sub-characteristic condition ( SC) is satisfied for the stochastic Galerkin formulations with unit vector e 1 = (1, 0, . . . , 0) T .

Numerical results
The introduction of the gPC modesv as auxiliary variable also allows for an efficient numerical evaluation of the flux function (13), the relaxation term (14) and the computation of eigenvalues by the numerically cheap and stable matrix vector multiplications Hence, the computational complexity grows like K 2 , which is relatively low compared to approaches with entropy and Roe variables [42,15]. The price is the restriction to gPC bases that satisfy the assumptions (A1) -(A3). Here, we use the Haar sequence [22,37,43] with level J ∈ N 0 that generates a gPC basis S K with K + 1 = 2 J+1 elements by Using a lexicographical order we identify the gPC basis φ 0 = 1, φ 1 = ψ, φ 2 = ψ 1,0 , φ 3 = ψ 1,1 , etc.
An equidistant space discretization ∆x > 0 is used to divide the space interval [0, x end ] into N cells such that ∆xN = x end with centers x j := j + 1 2 ∆x and edges x j− 1 /2 := j∆x. The discrete time steps are denoted by t k := k∆t for k ∈ N 0 . Due to the eigenvalue estimates λ 1 (ρ,ẑ) ≤ λ 2 (ρ,ẑ) = D v(ρ,ẑ) a local Lax-Friedrichs flux [35] is efficiently evaluated as For numerical purposes the relaxation term is expressed as Since the term Sẑ(û) depends also in the stochastic Galerkin formulation on the unknownẑ ∈ R K+1 in a linear way, a first-order IMEX scheme [40,23,44], which treats the advection part explicitly and the possibly stiff relaxation implicitly, can be employed: In the sequel, we consider a linear hesitation function and a relaxation towards the LWR model, i.e    ∂ t ρ + ∂ x (ρv) = 0, and normalized density in the equilibrium model. According to Corollary 4.2 the subcharacteristic condition is fulfilled and solutions to the ARZ model are expected to be close to the LWR model if the relaxation parameter τ > 0 is sufficiently small. Moreover reference solutions are provided, where a Monte-Carlo method is applied to the analytical solution with M = 10 6 uniformly distributed samples ρ (ξ) ∼ U for either of the following Riemann problems:

Homogeneous case
This section illustrates the hyperbolic character of the derived stochastic Galerkin formulation, in particular the statement of Theorem 3.1. Figure 1 and 2 illustrate the solution to the stochastic Galerkin formulation to the Haar basis with level J. The mean of the density is given by the modeρ 0 (t, x) and plotted as blue line. The confidence region to the truncated gPC expansion is black shaded. Furthermore, the Monte-Carlo confidence region is shown as black dotted line and the reference mean as green dashed line. We observe from Figure 1 for the rarefaction wave that the confidence region is already well captured for level J = 0 and the mean for J = 3. Likewise, Figure 2 shows the approximation for the shock case, when each realization admits a discontinuity. The mean, however, is smooth as an average of discontinuous functions. The stochastic Galerkin formulation approximates the mean as step functions (blue line). This behaviour is typical and has been observed also for continuous input distributions [8,41,15,14].

Inhomogeneous case
This section is devoted to the stability analysis in Section 4. We investigate the guaranteed dissipativity condition of Theorem 4.1 and Corollary 4.2, which presume a relaxation to a first-order model. Figure 3 and 4 show the behaviour of the inhomogeneous ARZ model for various relaxation parameter, including the limit τ = 0. The left panels show the results for the level J = 2 without relaxation and the exact confidence regions are plotted in the remaining panels for comparison. Indeed, we observe a convergence towards the LWR model according to Corollary 4.2.

Summary
A stochastic Galerkin formulation of the Aw-Rascle-Zhang (ARZ) model has been presented.
In particular, hyperbolicity has been shown for a special class of wavelet-based expansions. The analysis is based on a non-conservative formulation. This allows a stability analysis for the inhomogeneous ARZ with stiff relaxation, when solutions are expected to be close to an equilibrium velocity that satisfies a scalar conservation law. Due to the non-conservative formulation, the derived theoretical results hold only for smooth solutions. However, a relationship to a conservative form has been established. This allows for a numerical discretization with an IMEX scheme.