category: literaturenote citekey: brianinumericalcomputationgeneralizedwasserstein2024 title: Numerical computation of generalized Wasserstein distances with applications to traffic model analysis authors: "Briani, Maya; Cristiani, Emiliano; Franzina, Giovanni; Ignoto, Francesca L." year: 2024 date: 2024-00-00 2024 url: "https://arxiv.org/abs/2410.11441v2" zotero_key: VCF7ANB9 zotero_storage: W5ENIJ5G collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Briani, Maya"

status: converted

Numerical Computation of Generalized Wasserstein Distances with Applications to Traffic Model Analysis

Maya Briani^∗ , Emiliano Cristiani^∗,‡ , Giovanni Franzina^∗ , Francesca L. Ignoto^∗,†

March 18, 2025

Abstract

Generalized Wasserstein distances allow to quantitatively compare two continuous or atomic mass distributions with equal or different total mass. In this paper, we propose four numerical methods for the approximation of three different generalized Wasserstein distances introduced in the last years, giving some insights into their physical meaning. After that, we explore their usage in the context of the sensitivity analysis of differential models for traffic flow. The quantification of models sensitivity is obtained by computing the generalized Wasserstein distances between two (numerical) solutions corresponding to different inputs, including different boundary conditions.

MSC: 76A30; 90C31; 90C90; 35L50; 35L53.

| 1 | | Introduction | 2 | | | | | | |---|--------------------------------------------------------------|-----------------------|--------|--|--|--|--|--| | 2 | The Wasserstein distance in the context of vehicular traffic | | | | | | | | | | 2.1 | Motivations
| 5
5 | | | | | | | | 2.2 | Basic theory
| | | | | | | | | 2.3 | Numerical setting
| 6
7 | | | | | |

^∗ Istituto per le Applicazioni del Calcolo, Consiglio Nazionale delle Ricerche, Rome, Italy.

^†Dipartimento di Scienze di Base e Applicate per l'Ingegneria, Sapienza Universit`a di Roma, Rome, Italy.

^‡Corresponding author (e.cristiani@iac.cnr.it).

| 3 | Basic theory for Generalized Wasserstein distances | | | | | | | | |---|----------------------------------------------------|-----------------------------------------------------------------------------|----------|--|--|--|--|--| | | 3.1 | Figalli & Gigli approach [12]
| 10 | | | | | | | | 3.2 | Piccoli & Rossi approach [21]
| 11 | | | | | | | | 3.3 | Gaussian Hellinger–Kantorovich approach [16]
| 12 | | | | | | | | 3.4 | Savar´e & Sodini approach [25]
| 14 | | | | | | | 4 | | Discrete and computational approaches for Generalized Wasserstein distances | | | | | | | | | 4.1 | Figalli & Gigli approach
| 15 | | | | | | | | | 4.1.1
Discrete approach
| 15 | | | | | | | | | 4.1.2
Computational approach
| 15 | | | | | | | | | 4.1.3
Numerical tests
| 17 | | | | | | | | 4.2 | Piccoli & Rossi approach
| 19 | | | | | | | | | 4.2.1
Discrete approach
| 19 | | | | | | | | | 4.2.2
Computational approach
| 19 | | | | | | | | | 4.2.3
Numerical tests
| 20 | | | | | | | | 4.3 | Gaussian Hellinger–Kantorovich approach
| 22 | | | | | | | | | 4.3.1
Discrete approach
| 22 | | | | | | | | | 4.3.2
Computational approach
| 22 | | | | | | | | | 4.3.3
Numerical test
| 22 | | | | | | | | 4.4 | Savar´e & Sodini approach
| 23 | | | | | | | | | 4.4.1
Discrete approach
| 23 | | | | | | | | | 4.4.2
Computational approach
| 23 | | | | | | | | | 4.4.3
Numerical tests
| 26 | | | | | | | 5 | | Applications to traffic model analysis | 29 | | | | | | | | 5.1 | Mathematical modelling
| | | | | | | | | 5.2 | Numerical tests
| 29
31 | | | | | | | | | 5.2.1
Test 0: comparison of the three GWDs
| 31 | | | | | | | | | 5.2.2
Test 1: sensitivity to boundary conditions
| 32 | | | | | | | | | 5.2.3
Test 2: sensitivity to traffic light cycle
| 33 | | | | | | | | | 5.2.4
Test 3: sensitivity to model order
| 34 | | | | | | | | | | | | | | | | | 6 | | Discussion | 35 | | | | | |

1 Introduction

Context. In this paper, we are concerned with the numerical computation of Wasserstein distance (WD) and Generalized Wasserstein distance (GWD) in the context of traffic flow models. As for all dynamical systems, the study of the sensitivity of the models can be realized measuring the 'distance' between two solutions obtained with different inputs (like, e.g., any of the model parameters, initial conditions, boundary conditions, etc.). This allows one to understand their impact on the final solution, and ultimately quantify the degree of chaoticity of the system. The question arises which distance is more suitable to this kind of investigation. It is by now well understood that L ^p distances do not catch the natural concept of distance among traffic (vehicle) densities, see, e.g., the discussion in [9, Sec. 7.1] and Sec. 2 in this paper, while the WD appears more natural, at least in the context of traffic flow, see, e.g., [4, 9, 10]. The drawback is that WD is limited to balanced mass distribution (equal masses), while real traffic problems often need to consider scenarios with a different amount of vehicles, especially because of different inflow/outflow at boundaries. This suggests to move towards GWDs, which allows to deal with unbalanced mass distributions, but serious computational problems arise: all existing definitions are quite abstract and it is not immediate to derive suitable numerical approximations.

Relevant literature. A common feature of various optimal transport distances between finite positive Borel measures on a separable complete metric space (X, D) is that their primal definition requires the minimization of a transportation cost including an additive contribution of the form

$$\int_{X \times X} c(x, y) \, d\gamma \,, \tag{1.1}$$

where c : X × X → [0, +∞] is a fixed lower semi-continuous cost function.

In the balanced problem, that one considers to define the classical Wasserstein distance between two probability measures µ ^s and µ d , the minimization occurs upon transport plans γ, i.e., probability measures on X × X with marginals µ ^s and µ d . Namely, one considers the projection π s (resp., π d ) onto the first (resp., the second) factor of the Cartesian product X × X, and takes as admissible competitors the finite positive measures γ for which

$$\pi{\sharp}^{\mathrm{S}} \gamma = \mu^{\mathrm{S}}, \qquad \pi{\sharp}^{\mathrm{D}} \gamma = \mu^{\mathrm{D}},$$

$$(1.2)$$

where π s γ (resp., π d γ) denotes the push-forward measure of γ under π s (resp., π d ).

In order to define distances between finite positive measures on a bounded open set Ω in a finite dimensional Euclidean space, there exist various approaches to similar, but unbalanced, problems, that correspond each to the appropriately constrained minimization of a specific objective functional. In this paper we focus on four of such approaches.

1. The approach by Figalli & Gigli [12], which consists in minimizing the integral cost (1.1) among all finite positive measures γ on Ω × Ω whose marginals induce on Ω the prescribed finite positive Borel measures.
2. The approach by Piccoli & Rossi [21], which includes in the definition of the functional the additive contributions

$$|\pi{\sharp}^{s}\gamma - \mu^{s}|(\Omega) + |\pi{\sharp}^{D}\gamma - \mu^{D}|(\Omega), \qquad (1.3)$$

in addition to (1.1), with the minimization occurring upon finite positive Borel measures on Ω × Ω.

A third approach relies on the logarithmic entropy-transport distance that was independently introduced by Kondratyev, Monsaingeon & Vorotnikov in [15], by Chizat, Peyr´e, Schmitzer & Vialard in [7, 8], and by Liero, Mielke & Savar´e in [16], called Hellinger–Kantorovich or Wasserstein–Fisher–Rao distance. This approach also relies on the minimization, among finite positive Borel measures on Ω × Ω, of a perturbation of (1.1), except that the penalization (1.3) is replaced by

$$\int{\Omega} \left( \sigma^{\mathrm{S}} \log(\sigma^{\mathrm{S}}) - \sigma^{\mathrm{S}} + 1 \right) d\mu^{\mathrm{S}} + \int{\Omega} \left( \sigma^{\mathrm{D}} \log(\sigma^{\mathrm{D}}) - \sigma^{\mathrm{D}} + 1 \right) d\mu^{\mathrm{D}}, \tag{1.4}$$

if d(π s γ) = σ ^s dµ^s and d(π d γ) = σ ^d dµ^d for appropriate densities σ ^s and σ d , and +∞ otherwise. (Incidentally, a similar objective, in which an entropic regularization is added to (1.1) and (1.4), in considered in [19], where the appropriate dual formulation of the problem is addressed by Sinkhorn iterations.)

A fourth approach relates to the contribution by Savar´e & Sodini [25], in which the logarithmic entropy transport distance is characterized, in its various equivalent formulations, by convex relaxation. In particular, their approach can be considered the natural one to numerically approximate the (Gaussian) Hellinger–Kanrotovich distance in its homogeneous formulation.

We point out that the last two approaches only differ numerically, as in principle the corresponding unbalanced optimal transport formulations are analytically equivalent. Both extend WD to the cone of non-negative Borel measures of finite mass that are different from zero (see Rem. 5 below). On the contrary, the first two approaches do give an extension of WD, despite defining a distance on the cone (see Rems. 3 and 4 below). We recall that no extension of WD to the set of signed measures, possibly zero, can define a distance [18] (non-canonical optimal transport distance have however been successfully defined, see [1, 20]).

Main contribution. This paper has multiple goals: first of all, the four approaches are recalled and compared at the theoretical level, so to have a common and exhaustive reference for all recent developments about GWDs.

Second, four new ad hoc computational methods are devised in order to practically compute the four GWDs (actually they are three GWDs, since two of them coincide under suitable assumptions). The methods are individually studied and their main properties are discussed through some explanatory numerical tests.

Third, the four approaches are compared on the same academic test in order to understand which one is the most suitable in the context of traffic flow models.

Finally, the three most promising approaches are further compared on three realistic traffic problems.

Anticipating the conclusions, we can say that Figalli & Gigli approach appears to be preferable among all, mainly for two reasons: (i) it is the most manageable from the computational point of view (easy to code, fast computation), and (ii) it is the most informative on the modeling side.

Plan of the paper. The paper is organized as follows. In Sect. 2 we recall the main motivations of the paper and we set up notations for linear programming and discrete-to-discrete optimal mass transport. In Sect. 3 we recall the definition and the basic properties of some known unbalanced optimal transport-type distances between positive Borel measures of finite mass. In Sect. 4 we revisit the distances in the discrete-to-discrete case and, more important, we propose computational approaches for the approximation of the distances. We also study their interplay with the numerical boundary, and their capability to detect significant patterns, through numerical tests in which the exact value of the distances is available analytically. In Sect. 5 these conclusions are applied in the context of traffic analysis in order to achieve various goals: (i) to assess the impact of boundary conditions, (ii) to evaluate the sensitivity to traffic lights frequency variations, and (iii) to compare forecasts depending on whether or not the model takes inertial effects into account. The paper is concluded with a final discussion in Sect. 6.

2 The Wasserstein distance in the context of vehicular traffic

2.1 Motivations

Usually the quantification of the closeness of two time-dependent distributions is performed by means of the L 1 , L ² or L ^∞ distance (in space at final time or in both space and time). Although this can be satisfactory for nearly equal outputs or for convergence results, it appears inadequate for measuring the distance of largely different outputs. To see this, let us focus on the example depicted in Fig. 1, where three density functions ρ i , i = 1, 2, 3, corresponding to the same total

Figure 1: Three density functions with disjoint supports.

mass, say M, are plotted. It is plain that the L ¹ distances between ρ¹ and ρ² and between ρ¹ and ρ³ are both equal to 2M. Similarly, all L ^p distances are blind with respect to variation of the densities once the supports of them are disjoint. Our perception of distance suggests instead that ∥ρ³ − ρ1∥ > ∥ρ² − ρ1∥, and this is exactly what Wasserstein distance guarantees, as we will see in the next section.

2.2 Basic theory

Let us denote by (X, D) a locally compact complete and separable metric space with distance D, and by B(X) the Borel σ-algebra of (X, D). Let us also denote by M⁺(X) the set of non-negative finite measures on (X, B(X)). Let µ d (d standing for demand or destination) and µ s (s standing for supply or source) be two Radon measures in M⁺(X) such that µ d (X) = µ s (X) (same total mass).

Definition 1 (Wasserstein distance). For any p ∈ [1, +∞), the L p -Wasserstein distance between µ ^s and µ d is

$$Wp(\mu^{\mathrm{S}}, \mu^{\mathrm{D}}) := \left(\inf{\gamma \in \Gamma(\mu^{\mathrm{S}}, \mu^{\mathrm{D}})} \int_{X \times X} \mathfrak{D}(x, y)^p \ d\gamma(x, y)\right)^{1/p} \tag{2.1}$$

where Γ is the set of transport plans connecting µ s to µ d , i.e.

$$\Gamma(\mu^{\mathsf{s}}, \mu^{\mathsf{D}}) := \left{ \gamma \in \mathcal{Q}^+(X \times X) \text{ s.t. } \gamma(A \times X) = \mu^{\mathsf{s}}(A), \ \gamma(X \times B) = \mu^{\mathsf{D}}(B), \ \forall \ A, B \in \mathcal{B}(X) \right}.$$

Remark 1. The distance defined by (2.1) makes sense for any two given measures µ ^s and µ ^d on a finite dimensional normed space, whether they are absolutely continuous or not with respect to the Lebesgue measure.

It is well known that the notion of Wasserstein distance can be put in relation with the Monge– Kantorovich optimal mass transfer problem [28]: a pile of, say, soil, has to be moved to an excavation with same total volume. Moving a unit quantity of mass has a cost which equals the distance between the source and the destination point. We consequently are looking for a way to rearrange the first mass onto the second which requires minimum cost.

Remark 2. In our framework, the mass to be moved corresponds to that of vehicles. We therefore measure the distance between two traffic density distributions by computing the minimal cost to rearrange vehicles from the first configuration onto the second one. We also stress that, since we are considering macroscopic models, vehicles are indistinguishable so we are not able to pinpoint the exact same vehicle in the two scenarios, cf. [10].

In the particular case of 1D optimal transport various characterizations give alternative, more manageable, definitions of WD. For example, if (X, D) is R with the euclidean distance | · |, then, in the case of two Dirac delta functions δx, δ^y we simply have Wp(δx, δy) = D(x, y) := |y − x|. That coincides with the desired distance quantification in the case of only two vehicles on a road, located in x and y respectively. Instead, in the case of two absolutely continuous measures with dµ^s = ρ ^sdx and dµ^d = ρ ^ddx, by [27, Rem. 2.19] we have

$$W1(\mu^{\rm S}, \mu^{\rm D}) = \int{\mathbb{R}} |F^{\rm S}(x) - F^{\rm D}(x)| dx = \int_0^1 |(F^{\rm S})^{-1}(t) - (F^{\rm D})^{-1}(t)| dt, \tag{2.2}$$

with

$$F^{\mathrm{s}}(x) := \int{-\infty}^{x} \rho^{\mathrm{s}}(z)dz, \quad F^{\mathrm{d}}(x) := \int{-\infty}^{x} \rho^{\mathrm{d}}(z)dz \tag{2.3}$$

for p = 1, whereas for p ≥ 1,

$$Wp(\mu^{\rm S}, \mu^{\rm D}) = \left( \int{\mathbb{R}} |T^*(x) - x|^p \, \rho^{\rm S}(x) dx \right)^{1/p},$$ (2.4)

where T ∗ : R → R satisfies

$$\int{-\infty}^{x} \rho^{\mathbf{s}}(z)dz = \int{-\infty}^{T^{*}(x)} \rho^{\mathbf{D}}(z)dz \quad \forall x \in \mathbb{R}.$$ (2.5)

Eqs. (2.4)-(2.5) translate the fundamental property of monotone rearrangement of the mass [27, Rem. 2.19], i.e. the optimal strategy consists in transferring the mass starting from the left, see Fig. 2.

Figure 2: Property of monotone rearrangement.

2.3 Numerical setting

Let us introduce here the numerical setting that will accompany us throughout the paper, see Fig. 3. For simplicity, we will only consider the one-dimensional interval X = [xmin, xmax], divided in N equispaced cells C1, . . . , C^N (N is our first parameter of discretization, and we expect to vanish the numerical error for N → +∞). We also denote by xⁱ the barycenter of the cell Cⁱ , i = 1, . . . , N and by ∆x := xi+1 − xⁱ the step of the space discretization. To measure the distance between two cells Cⁱ and C^j , we will always use the Euclidean norm D(x^j , xk) := |x^k − x^j |.

Moving to the masses, for any i = 1, . . . , N we denote by

$$m_i^{\mathrm{S}} := \rho^{\mathrm{S}}(x_i)\Delta x, \qquad m_i^{\mathrm{D}} := \rho^{\mathrm{D}}(x_i)\Delta x,$$ (2.6)

the supply and the demand mass located in Cⁱ , respectively. Roughly speaking, this corresponds to concentrate the masses on specific points, thus getting purely atomic density distributions ˆρ s :=

Figure 3: Numerical setting for a generic one-dimensional problem. Two copies of the domain X with supply mass only (left) and demand mass only (right).

P ⁱ ^m^s i δxi and ˆρ d := P ⁱ ^m^d i δ^xⁱ where δ denotes the Dirac delta. Note that, in our context, the vectors {m^s ⁱ}ⁱ and {m^d ⁱ }ⁱ will be given. Moreover, to avoid any confusion between theoretical and numerical boundary conditions we will always assume that

$$m_1^{\rm S} = m_N^{\rm S} = m_1^{\rm D} = m_N^{\rm D} = 0,$$ (2.7)

in such a way that boundaries remain out of play.

For some reasons which will be clear later on, it is also useful defining the remaining supply mass at Cⁱ , denoted by ˚m^s i , and the remaining demand mass at Cⁱ , denoted by ˚m^d i , i = 1, . . . , N.

In order to get a numerical approximation of the WD in R 1 , it is convenient to resort to the definition (2.2). Indeed, any quadrature formula is suitable to get the desired value with arbitrary precision, provided ρ s , ρ^d are sufficiently smooth and N is sufficiently large. For example, the naive composite mid-point formula gives

$$W_1(\rho^{\rm S}, \rho^{\rm D}) \approx \tilde{W}1(m^{\rm S}, m^{\rm D}) := \Delta x \sum{j=1}^N \left| \sum_{i=1}^j mi^{\rm S} - \sum{i=1}^j m_i^{\rm D} \right|.$$ (2.8)

The case of more general spaces, including networks embedded in R , is more difficult. Typically the problem is solved by reformulating it in a fully discrete setting as a Linear Programming (LP) problem. The interested reader can find in the books by Santambrogio [24, Sect. 6.4.1] and Sinha [26, Chapt. 19] the complete procedure, recently used also in the context of traffic flow on networks [4]. Even if we will stick with the one-dimensional case, we recall now the crucial steps of this approach because it will be useful in the following.

Following the mass transport interpretation, let cjk be the cost of shipping a unit quantity of mass from the source C^j , j = 1, . . . , N, to the destination Ck, k = 1, . . . , N. In general, cjk is usually defined as the length of the shortest path joining x^j and xk. Let γjk be the (unknown) quantity shipped from the source j to the destination k. The problem is then formulated as

minimize $$\sum{j=1}^{N} \sum{k=1}^{N} c{jk} \gamma{jk}$$ subject to $$\sum{k} \gamma{jk} = m_{j}^{s}, \quad \forall j$$

$$\sum{j} \gamma{jk} = m_{k}^{p}, \quad \forall k$$

$$\gamma_{jk} \ge 0.$$ (2.9)

The solution gives an approximation W˜ ¹ of W1. Note that the solution satisfies γjk ≤ min{m^s j , m^d k } since one cannot move more than m^s j from any source j and it is useless to bring more than m^d k to any destination k. From (2.9), it is easy to recover a standard LP problem with equality constraints

minimize $$\mathbf{c}^{\mathsf{T}}\mathbf{x}$$

subject to $\mathbf{A}\mathbf{x} = \mathbf{b}$ (2.10)
$\mathbf{x} \ge 0$ ,

simply defining

$$\mathbf{x} := (\gamma{11}, \gamma{12}, \dots, \gamma{1N}, \gamma{21}, \gamma{22}, \dots, \gamma{2N}, \dots, \gamma{N1}, \dots, \gamma{NN})^{\mathsf{T}} \mathbf{c} := (c{11}, c{12}, \dots, c{1N}, c{21}, c{22}, \dots, c{2N}, \dots, c{N1}, \dots, c{NN})^{\mathsf{T}} \mathbf{b} := (m_1^{\mathsf{S}}, \dots, m_N^{\mathsf{S}}, m_1^{\mathsf{D}}, \dots, m_N^{\mathsf{D}})^{\mathsf{T}}$$ (2.11)

and A as the 2N × N² block matrix

$$\mathbf{A} := \begin{bmatrix} \mathbf{1}{N} & \mathbf{0}{N} & \mathbf{0}{N} & \cdots & \mathbf{0}{N} \ \mathbf{0}{N} & \mathbf{1}{N} & \mathbf{0}{N} & \cdots & \mathbf{0}{N} \ \mathbf{0}{N} & \mathbf{0}{N} & \mathbf{1}{N} & \cdots & \mathbf{0}{N} \ \vdots & \vdots & \vdots & \ddots & \vdots \ \mathbf{0}{N} & \mathbf{0}{N} & \mathbf{0}{N} & \cdots & \mathbf{1}{N} \ I{N} & I{N} & I{N} & I{N} & I_{N} \end{bmatrix},$$

$$(2.12)$$

where I^N is the N × N identity matrix, 1^N := (1 1 · · · 1 | {z } N times ), and 0^N := (0 0 · · · 0 | {z } N times ).

3 Basic theory for Generalized Wasserstein distances

In this section we recall some definitions of metrics on the space of positive finite Borel measures on the locally compact complete and separable metric space (X, D), inducing the weak-narrow convergence, all inspired to the distance defined by (2.1) between probability measures.

3.1 Figalli & Gigli approach [12]

In [12], an optimal transport type distance is introduced, with respect to which the heat equation under Dirichlet conditions at the boundary of an open bounded domain Ω can be described as the gradient flow of an entropy, similarly to the case of the heat equation in free space. In order to encode the correct boundary conditions in the minimizing movement scheme, the GWD between two finite non-negative (Borel) measures µ s , µ^d on Ω is defined by

$$Wp^{\text{FG}}(\mu^{\text{S}}, \mu^{\text{D}}) = \left{ \inf{\gamma} \int_{\overline{\Omega} \times \overline{\Omega}} |x - y|^p d\gamma(x, y) \right}^{\frac{1}{p}}, \tag{3.1}$$

where infimum is over

$$\left{ \gamma \in \mathcal{M}^{+}(\overline{\Omega} \times \overline{\Omega}) \mid \left( \pi{\sharp}^{s} \gamma \right) \sqcup \Omega = \mu^{s}, \ \left( \pi{\sharp}^{D} \gamma \right) \sqcup \Omega = \mu^{D} \right}. \tag{3.2}$$

In the latter, we used the notation (¯µ Ω)(A) = ¯µ(A ∩ Ω) for the restriction to Ω of a measure ¯µ on Ω, and for all given γ we denoted by π s γ (resp., π d γ) the push-forward of γ under the projection defined on Ω × Ω by π s (x, y) = x (resp., π d (x, y) = y). Equivalently, a non-negative finite Radon measure γ on Ω × Ω is an admissible competitor for the minimization in (3.1) if and only if

$$\gamma(A \times \overline{\Omega}) = \mu^{s}(A), \ \gamma(\overline{\Omega} \times B) = \mu^{d}(B), \quad \text{for all } A, B \in \mathcal{B}(\Omega).$$ (3.3)

Let us note that (3.2) coincides with the set

$$\bigcup{(\overline{\mu}^{S},\overline{\mu}^{D})} \Big{ \gamma \in \mathcal{M}^{+}(\overline{\Omega} \times \overline{\Omega}) \; \Big| \; \pi^{S}{\sharp} \gamma = \bar{\mu}^{S} \, , \; \pi^{D}_{\sharp} \gamma = \bar{\mu}^{D} \Big},$$

where the union is over all pairs (¯µ s , µ¯ d ) ∈ M⁺(Ω) × M⁺(Ω) such that

$$\bar{\mu}^{\mathrm{S}} \, \sqcup \, \Omega = \mu^{\mathrm{S}} \,, \ \bar{\mu}^{\mathrm{D}} \, \sqcup \, \Omega = \mu^{\mathrm{D}} \,, \ \bar{\mu}^{\mathrm{S}}(\overline{\Omega}) = \bar{\mu}^{\mathrm{D}}(\overline{\Omega}) \,.$$

Hence, an equivalent variational definition of the unbalanced optimal transport distance (3.1) is

$$\min \left{ W_p(\bar{\mu}^{\mathrm{S}}, \bar{\mu}^{\mathrm{D}}) \mid \bar{\mu}^{\mathrm{S}} \sqcup \Omega = \mu^{\mathrm{S}}, \ \bar{\mu}^{\mathrm{D}} \sqcup \Omega = \mu^{\mathrm{D}}, \ \bar{\mu}^{\mathrm{S}}(\partial \Omega) - \bar{\mu}^{\mathrm{D}}(\partial \Omega) = \mu^{\mathrm{S}}(\Omega) - \mu^{\mathrm{D}}(\Omega) \right}, \tag{3.4}$$

where the classical WD is defined as in (2.1) with D(x, y) = ∥x − y∥ for probability measures, and extended by homogeneity to non-negative Radon measures with equal masses. With this definition, if ¯µ s , µ¯ ^d ∈ M⁺(Ω) are fixed, then Wp(tµ¯ s , tµ¯ d ) is arbitrarily large if t > 0 is large enough. Hence the minimum in (3.4) is indeed achieved. Notice that mass unbalance µ s (Ω) − µ d (Ω) is fixed if so are the marginals µ ^s and µ d

The authors of [12] deal explicitly with the case of a quadratic unitary cost, i.e., c(x, y) = ∥x−y∥ 2 in (3.1), but their main results also hold for a more general cost, e.g., c(x, y) = ∥x−y∥ p , with p ≥ 1. As intermediate steps towards their main results, they prove lower semicontinuity with respect to the weak-∗ convergence on Ω and provide a gluing lemma implying the triangle inequality, so that (3.1) defines a distance on M⁺(Ω), and the infimum is achieved by an appropriate optimal plan γ.

The idea of minimizing (3.1) among all γ as in (3.2) is that ∂Ω models an infinite mass reservoir and the transportation policy described by convenient competitors γ might include the option to cut out part of the material from the actual transference plan, by transferring it to ∂Ω, as well as the possibility to integrate the transportation plan by transferring material from ∂Ω to the final destination, provided that the transportation cost to and from ∂Ω is accounted for in the overall integral cost together with that of moving the mass that is authentically transported. A dual formulation of the problem leads to a generalization of Kantorovich potentials and Kantorovich norms, which had been already investigated especially in the case of a linear cost (p = 1), see [13, 14].

Remark 3. It is clear from (3.4) that the distance between positive finite Radon measures introduced in [12] is not an extension of the classical Wasserstein distance: if the two prescribed marginals are concentrated close to the boundary with supports far apart from each other, then while solving (3.4) there will be a convenience for ¯µ ^s and ¯µ d in charging appropriately two corresponding disjoint boundary regions, after which the Wasserstein to be minimized will be achieved by transport plans that only move mass from or to the boundary. In fact, the goal of [12] was to generalize Otto's calculus to PDEs under Dirichlet conditions, rather than extending the optimal transport metrics to finite positive Radon measures.

3.2 Piccoli & Rossi approach [21]

In [21], the set of positive finite Borel measures on Ω is endowed with a metrics, defined by

$$Wp^{\text{PR}}(\mu^{\text{S}}, \mu^{\text{D}}; a, b) := \inf{\gamma} \left{ a \left| \mu^{\text{S}}(\Omega) - \pi{\sharp}^{\text{S}} \gamma(\Omega) \right| + a \left| \mu^{\text{D}}(\Omega) - \pi{\sharp}^{\text{D}} \gamma(\Omega) \right| + b \int_{\Omega \times \Omega} |x - y|^p \, d\gamma \right}^{\frac{1}{p}}, \quad (3.5)$$

where a, b > 0 are fixed parameters and the infimum is among all finite positive Borel measures γ on Ω × Ω.

Note that (3.5) can be equivalently written in the form

$$W_p^{\mathrm{PR}}(\mu^{\mathrm{S}}, \mu^{\mathrm{D}}; a, b) = \inf \left{ a \left( \left| \mu^{\mathrm{S}}(\Omega) - \mathring{\mu}^{\mathrm{S}}(\Omega) \right| + \left| \mu^{\mathrm{D}}(\Omega) - \mathring{\mu}^{\mathrm{D}}(\Omega) \right| \right) + b W_p \left( \mathring{\mu}^{\mathrm{S}}, \mathring{\mu}^{\mathrm{D}} \right) \, \middle| \, \mathring{\mu}^{\mathrm{S}}(\Omega) = \mathring{\mu}^{\mathrm{D}}(\Omega) \right}, \tag{3.6}$$

where µ ^s and µ d form an arbitrary pair of positive finite Borel measure on Ω with equal mass. In the right hand-side of (3.6), the term Wp(˚µ s ,˚µ d ) is understood as the value taken by the extension by homogeneity to pairs of finite positive measures of equal mass of the Wasserstein distance between probability measures.1

¹Namely, we are setting Wp(˚µ ,˚µ ) = mWp(˚µ ^s/m,˚µ ^d/m), where m = ˚µ d (Ω) = ˚µ d (Ω) and ˚µ s ,˚µ ^d are thought of as measures on Ω that vanish on all measurable subsets of the boundary ∂Ω = Ω \ Ω.

The infimum in (3.6), and hence in (3.5) too, is indeed achieved and it defines a distance [21, Prop. 1].

The mass transfer may include either physical transportation or creation and destruction. Adding mass to and removing it from µ s , as well as from µ d , are both allowed actions, with the implied costs being weighted by the parameter a. After these actions have been taken, the supplies eligible for transport, as well as the actual demands, are established and the corresponding optimal transfer policy contributes to the overall cost, via a classical WD weighted by b.

The constraints ˚µ s (Ω) ≤ µ s (Ω) and ˚µ d (Ω) ≤ µ d (Ω) may be added in the minimization without changing the minimum value, see [22, Prop. 2].

Remark 4. Heuristically, all combinations of the above actions are admissible for the minimization: an arbitrary fraction of mass may be removed from the marginals, along with transporting the remaining mass. Accordingly, (3.5) does not extend the Wasserstein distance because in the case of two marginals with equal masses the allowed actions include the possibility of removing mass from both. On the other hand, in the case p = 1 with a = b = 1 it still offers a natural generalization of the 1-WD, being equivalent to the flat metric.

The Piccoli & Rossi distance (3.5) differs from the Figalli & Gigli one (3.1). Recalling (3.4), it is clear that, in the latter, the mass added to the marginals is located at the boundary and that no other pay-off is required except for the transportation cost between all locations, including those at the boundary. By contrast, the boundary geometry plays no role in (3.5).

On the other hand, at least in the case of a quadratic cost (p = 2), it is known (see [7, Prop. 1.1]) that (3.5) is equivalent to the partial optimal transport problem [6, 11], i.e.,

$$\inf{\gamma} \int{\Omega \times \Omega} \left( |x - y|^2 - \lambda \right) d\gamma \,,$$

where the infimum is among all non-negative Radon measures on Ω × Ω, whose left and right marginals are dominated by µ ^s and µ d , that are subject to have a fixed overall mass m, and the Lagrange multiplier for this constraint satisfies λ = 2a ²/b² .

Furthermore, the distance (3.5) is inspired by the dynamical formulation of the Monge–Kantorovich problem: by generalizing the Benamou and Brenier's approach [3], the supply and the demand distributions are thought of as the initial and the final prescribed configuration for a fluid motion and the minimal value of (3.5) is obtained by minimizing the total kinetic energy necessary to connect the two configurations so as to balance mass conservation with creation and destruction of mass, thanks to the presence of an additional source term in a differential constraint for the macroscopic fluid variables.

3.3 Gaussian Hellinger–Kantorovich approach [16]

Another distance between positive finite Borel measures on Ω is defined by

$$W2^{\text{GHK}}(\mu^{\text{S}}, \mu^{\text{D}}; a, b) := \inf \left{ a \mathcal{U}\left(\pi{\sharp}^{\text{S}} \gamma \mid \mu^{\text{S}}\right) + a \mathcal{U}\left(\pi{\sharp}^{\text{D}} \gamma \mid \mu^{\text{D}}\right) + b \int{\Omega \times \Omega} |x - y|^2 \, d\gamma \right}^{\frac{1}{2}}, \tag{3.7a}$$

with a, b > 0 being fixed parameters, and the infimum is among all positive finite measures γ on Ω × Ω. Here,

$$\mathcal{U}(\tilde{\mu} \mid \mu) = \begin{cases} \int U(\sigma) d\mu, & \text{if } d\tilde{\mu} = \sigma d\mu, \ +\infty, & \text{otherwise,} \end{cases}, \quad \text{where } U(\sigma) = \sigma \log(\sigma) - \sigma + 1. \tag{3.7b}$$

Such a distance is called in [16] the Gaussian Hellinger–Kantorovich distance in view of the occurrence of a Gaussian function in its homogeneous formulation, see Sect. 3.4 below. We note that (3.7) can be equivalently written in the form

$$W_2^{\text{GHK}}(\mu^{\text{S}}, \mu^{\text{D}}; a, b) = \inf \left{ a \mathcal{U}(\hat{\mu}^{\text{S}} \mid \mu^{\text{S}}) + a \mathcal{U}(\hat{\mu}^{\text{D}} \mid \mu^{\text{D}}) + b W_2^2(\hat{\mu}^{\text{S}}, \hat{\mu}^{\text{D}}) \right}^{\frac{1}{2}}, \tag{3.8}$$

where the infimum is among all positive finite measures ˆµ ^s and ˆµ ^d of equal mass on Ω. In order to define the last term in the right hand-side of (3.8), we have implicitly understood both ˆµ ^s and ˆµ ^d as multiples of probability measures on Ω vanishing on ∂Ω and we have defined their WD accordingly, with the obvious extension by homogeneity to pairs of finite positive measures of equal mass on Ω of the WD between probability measures.

When the given measures µ ^s and µ d themselves happen to have equal mass, they form an admissible pair for the minimization: then, the right hand side of (3.8) is the classical (quadratic) WD, whose the value is a computed by solving a balanced optimal transport problem, because in view of (3.7b) we have U(1) = 0. In the unbalanced case, the two additional non-negative contributions in the objective penalize the deviation of the competitors from the prescribed measures.

The variational principle (3.7) belongs to the wide class of optimal entropy-transport problems considered in [16], in which the penalization is proportional to the Csisz´ar–Morimoto divergence, defined in terms of proper convex functions U : [0, +∞) → [0, +∞], called the entropy functions, possibly different from (3.7b). That large class covers the Piccoli & Rossi distance (3.6) for p = 1, with entropy function U(σ) = |σ − 1|, and is comprised of several models, alternative to [21] but equally inspired by dynamical formulations with source terms in the differential constraints, which have appeared in the literature since.

In particular, the model case of the logarithmic optimal entropy-transport problem (3.7) has been investigated under various perspectives: its value coincides with that of the metrics studied in [7, 15] (also called Wasserstein–Fisher–Rao distance, as it interpolates between the quadratic Wasserstein and the Fisher–Rao distance), for which a static definition had been provided in [8].

Remark 5. It is known, e.g. by [16, Part II], that the logarithmic entropy-transport minimization (3.7) defines a distance that metrizes the weak/narrow convergence of non-negative finite measure, extending the quadratic WD (otherwise defined only for probability measures).

3.4 Savar´e & Sodini approach [25]

In [16, Part II] the logarithmic entropy-transport minimization (3.7) is proved to be equivalent to an optimal transport type problem in the metric cone C[X] over X. More precisely, one defines C[X] as the quotient of X × [0, +∞) obtained by collapsing the whole fiber X × {0} to a single element, endowed with a distance that is defined, once an appropriate function g is chosen, by setting

$$d_{\mathfrak{C},g}([x^{\scriptscriptstyle \mathrm{S}},r^{\scriptscriptstyle \mathrm{S}}],[x^{\scriptscriptstyle \mathrm{D}},r^{\scriptscriptstyle \mathrm{D}}]) = \sqrt{H_g(x^{\scriptscriptstyle \mathrm{S}},r^{\scriptscriptstyle \mathrm{S}};x^{\scriptscriptstyle \mathrm{D}},r^{\scriptscriptstyle \mathrm{D}})}\,,$$

where

$$H_g(x, r; y, s) = r + s - 2\sqrt{rs} \ g(|x - y|) \text{ and } g(t) = e^{-\frac{t^2}{2}}.$$ (3.9)

Then, one considers the minimization of

$$\int{\mathfrak{C}[X]\times\mathfrak{C}[X]} d{\mathfrak{C},g}^2([x^{\scriptscriptstyle \mathrm{S}},r^{\scriptscriptstyle \mathrm{S}}],[x^{\scriptscriptstyle \mathrm{D}},r^{\scriptscriptstyle \mathrm{D}}])\,d\alpha$$

among positive measures α on C[X] × C[X] with homogeneous marginals µ ^s and µ d , i.e., such that

$$\mu^{\mathrm{S}}(A) = \int{\mathfrak{C}[A] \times \mathfrak{C}[X]} r^{\mathrm{S}} \, d\alpha([x^{\mathrm{S}}, r^{\mathrm{S}}], [x^{\mathrm{D}}, r^{\mathrm{D}}]) \,, \quad \mu^{\mathrm{D}}(B) = \int{\mathfrak{C}[X] \times \mathfrak{C}[B]} r^{\mathrm{D}} \, d\alpha([x^{\mathrm{S}}, r^{\mathrm{S}}], [x^{\mathrm{D}}, r^{\mathrm{D}}]) \,.$$

for all Borel sets A, B ⊂ X.

The constrained minimum is achieved, and it coincides with the Gaussian Hellinger–Kantorovich distance defined by (3.7) in the case a = b = 1.2

4 Discrete and computational approaches for Generalized Wasserstein distances

In this section, we first revise the GWDs previously introduced in the special case of atomic measures (sum of Dirac delta's). This allows us to naturally introduce a discretization in the space domain, which is the first step towards a numerical computation. Secondly, we propose original computational approaches to get algorithms ready to be implemented on a computer. Finally, we test the algorithms on some academic examples which are devised to highlight the specific features of each GWD.

²A match is possible also for different choices of parameters a, b, provided that the definition of H^g changes accordingly in (3.9). Moreover, the choice g(t) = cos t ² ∧ 2 is also possible in (3.9), leading to the Hellinger– Kantorovich distance, which equals (3.7a) if (3.7b) is replaced by the entropy function U(σ) = σ − 1 − log σ.

4.1 Figalli & Gigli approach

4.1.1 Discrete approach

In the special case of two purely atomic marginals with finite supports in a one-dimensional domain, i.e.,

$$\mu^{\mathrm{S}} = \sum_{i=2}^{N-1} mi^{\mathrm{S}} \delta{xi} \,, \qquad \mu^{\mathrm{D}} = \sum{i=2}^{N-1} mi^{\mathrm{D}} \delta{x_i} \label{eq:multiple}$$

we get an interpretation of (3.4) as a linear programming problem. Namely, we consider the problem

$$\min_{(m_1^{\mathrm{S}}, m_N^{\mathrm{D}}, m_N^{\mathrm{D}}, m_N^{\mathrm{D}}) \in \mathbb{R}^4} \left{ Wp^p \Big( \sum{i=1}^N mi^{\mathrm{S}} \delta{xi} \; , \sum{i=1}^N mi^{\mathrm{D}} \delta{x_i} \Big) \; \Big| \; m_1^{\mathrm{S}} + m_N^{\mathrm{S}} - (m_1^{\mathrm{D}} + m_N^{\mathrm{D}}) = \Delta M \right} \; ,$$

where

$$\Delta M := \sum_{i=2}^{N-1} mi^{\scriptscriptstyle{\mathrm{S}}} - \sum{i=2}^{N-1} m_i^{\scriptscriptstyle{\mathrm{D}}},$$

and, after setting cjk = ∥x^j − xk∥ p , for all j, k ∈ {1, . . . , N}, we write it in the equivalent form

$$\tilde{W}{1}^{\text{FG}}(m^{\text{S}}, m^{\text{D}}) := \min \left{ \sum{j,k=1}^{N} c{jk} \gamma{jk} \mid \gamma \in [0, +\infty)^{N \times N}, \sum{k=1}^{N} \gamma{jk} = m{j}^{\text{S}}, \text{ for } j \in {2, \dots, N-1}, \ \sum{j=1}^{N} \gamma{jk} = m{k}^{\text{D}}, \text{ for } k \in {2, \dots, N-1}, \ \sum{k=1}^{N} \left( \gamma{1k} + \gamma{Nk} \right) - \sum{j=1}^{N} \left( \gamma{j1} + \gamma{jN} \right) = \Delta M \right}^{\frac{1}{p}}$$

That amounts to the minimization of a linear objective in a convex polytope obtained by intersecting 2N − 3 affine half-spaces in R ^N×^N with [0, +∞) ^N×^N . Alternatively, the last constraint may be converted into an additional contribution to the total cost involving suitable Lagrange multipliers.

4.1.2 Computational approach

The minimization problem introduced in the previous section can be formulated as a LP problem in the standard form (2.10) defining

$$\begin{aligned} \mathbf{x} &:= (\gamma{11}, \gamma{12}, \dots, \gamma{1N}, \gamma{21}, \gamma{22}, \dots, \gamma{2N}, \dots, \gamma{N1}, \dots, \gamma{NN})^{\mathsf{T}} \ \mathbf{c} &:= (c{11}, c{12}, \dots, c{1N}, c{21}, c{22}, \dots, c{2N}, \dots, c{N1}, \dots, c{NN})^{\mathsf{T}} \ \mathbf{b} &:= (m{2}^{\mathsf{S}}, \dots, m{N-1}^{\mathsf{S}}, m{2}^{\mathsf{D}}, \dots, m{N-1}^{\mathsf{D}}, \Delta M)^{\mathsf{T}} \end{aligned}$$

and A as the (2N − 3) × N² block matrix

$$\mathbf{A} := \begin{bmatrix} 0{N} & \mathbb{1}{N} & 0{N} & \cdots & 0{N} & 0{N} \ 0{N} & 0{N} & \mathbb{1}{N} & \cdots & 0{N} & 0{N} \ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \ 0{N} & 0{N} & 0{N} & \cdots & \mathbb{1}{N} & 0{N} \ J{N} & J{N} & J{N} & \cdots & J{N} & J{N} \ U{N} & Z{N} & Z{N} & \cdots & Z{N} & U_{N} \end{bmatrix},$$

$$(4.1)$$

where I^N is the N×N identity matrix, 1^N := (1 1 · · · 1 | {z } N times ), 0^N := (0 0 · · · 0 | {z } N times ), Z^N := (−1 0 0 · · · 0 | {z } N−2 times −1), U^N := (0 1 1 · · · 1 | {z } N−2 times 0), and J^N is the (N − 2) × N block matrix [0 ⊺ N−2 IN−² 0 ⊺ N−2 ].

The LP formulation described above can be reformulated in an almost equivalent form, which is actually more similar to the problem (2.9) and also allows for a nice, further, generalization. Still keeping the simplifying assumption (2.7), the vector x of the unknowns is modified adding four new unknowns M^s 1 , M^s ^N , M^d 1 , M^d ^N , corresponding to the masses taken from the left/right boundaries, and the masses brought to the left/right boundaries, respectively. Note that these masses are potentially unlimited. Moreover, four input parameters c s 1 , c s ^N , c d 1 , c d ^N related to the four new unknowns are added. For example, the quantity c s 1M^s ¹ will give the cost to create the mass M^s 1 at the left boundary (N.B., just creating it, not transporting it inside the domain), while the quantity c d NM^d ^N will give the cost to destroy the mass M^d ^N at the right boundary.

The PL formulation becomes:

$$\tilde{W}{1}^{\text{FC}'}(m^{\text{S}}, m^{\text{D}}) := \min \left{ \sum{j,k=1}^{N} c{jk} \gamma{jk} + c{1}^{\text{S}} M{1}^{\text{S}} + c{N}^{\text{S}} M{N}^{\text{S}} + c{1}^{\text{D}} M{1}^{\text{D}} + c{N}^{\text{D}} M{N}^{\text{D}} \, \Big| \, \gamma \in [0, +\infty)^{N \times N}, \ \sum{k=1}^{N} \gamma{jk} = m{j}^{\text{S}}, \text{ for } j \in {2, \dots, N-1}, \quad \sum{j=1}^{N} \gamma{jk} = m{k}^{\text{D}}, \text{ for } k \in {2, \dots, N-1}, \ \sum{k=1}^{N} x{jk} = M{j}^{\text{S}}, \text{ for } j \in {1, N}, \quad \sum{j=1}^{N} x{jk} = M{k}^{\text{D}}, \text{ for } k \in {1, N}, \ M{1}^{\text{S}}, M{N}^{\text{S}}, M{1}^{\text{D}}, M{N}^{\text{D}} \ge 0, \quad c{11} = c{1N} = c{N1} = c{NN} = +\infty \right}^{\frac{1}{p}}, \tag{4.2}$$

which corresponds to (2.10) with

$$\begin{aligned} \mathbf{x} &:= (x{11}, x{12}, \dots, x{1N}, x{21}, x{22}, \dots, x{2N}, \dots, x{N1}, \dots, x{NN}, M_1^{\mathrm{S}}, M_N^{\mathrm{S}}, M_1^{\mathrm{D}}, MN^{\mathrm{D}})^{\mathsf{T}} \ \mathbf{c} &:= (c{11}, c{12}, \dots, c{1N}, c{21}, c{22}, \dots, c{2N}, \dots, c{N1}, \dots, c_{NN}, c_1^{\mathrm{S}}, c_N^{\mathrm{S}}, c_N^{\mathrm{D}}, c_N^{\mathrm{D}})^{\mathsf{T}} \ \mathbf{b} &:= (0, m2^{\mathrm{S}}, \dots, m{N-1}^{\mathrm{S}}, 0, 0, m2^{\mathrm{D}}, \dots, m{N-1}^{\mathrm{D}}, 0)^{\mathsf{T}} \end{aligned}$$

and A is a 2N × (N² + 4) block matrix, obtained from matrix (2.12) by extending it with four columns. The extra columns have a single non-zero element each, namely

$$\mathbf{A}{1,N^2+1} = -1, \quad \mathbf{A}{N,N^2+2} = -1, \quad \mathbf{A}{N+1,N^2+3} = -1, \quad \mathbf{A}{2N,N^2+4} = -1,$$ (4.3)

due to the four extra unknowns and constraints of the problem (4.2). The last constraints of (4.2) are needed to assure that it never convenient creating mass at boundaries just to keep it there or to transport it on the other boundary.

4.1.3 Numerical tests

Test 1

In this test we consider the case of atomic balanced masses in the domain X = [0, 5]. We set N = 50 so that ∆x = 0.1. Total supply mass is 0.2, divided in two equal parts and concentrated in two points. Total demand mass is also 0.2, concentrated in one point only.

1.1. Here supply mass is located at x = 1 and x = 2.5, while demand mass is located at x = 3, see Fig. 4(left). Optimal mass transport map prescribes to move supply mass onto demand mass, without resorting to mass at boundaries. We have W˜ fg = 0.25.
1.2. Here supply mass is located at x = 0.5 and x = 2.5, while demand mass is located at x = 4, see Fig. 4(right). Optimal mass transport map prescribes to move leftmost supply mass to the left boundary, rightmost supply mass onto the demand mass, and taking mass from right boundaries to bring them onto the demand mass. We have W˜ fg = 0.3.

Test 2

In this test we consider the case of continuous unbalanced distributions in the domain X = [−5, 5],

$$\rho^{\text{S}}(x) = e^{-x^2}, \qquad \rho^{\text{D}}(x) = \frac{e^{-(x-3)^2}}{4},$$

see Fig. 5. We set N = 50 so that ∆x = 0.2.

In this case the optimal mass transport map prescribes to move a part of the supply mass to the let boundary and a part of the demand mass to the right boundary, then transport the remaining mass. We have W˜ fg = 6.84.

Figure 4: Test 1. Case 1.1 where boundaries are not used (left), and Case 1.2 where mass is both taken from and brought to boundaries (right).

Figure 5: Test 2. Some supply (resp., demand) mass is brought to left (resp., right) boundary, the rest is transported.

4.2 Piccoli & Rossi approach

4.2.1 Discrete approach

In the case of two purely atomic marginals with finite supports

$$\mu^{\mathrm{S}} = \sum_{i=1}^N mi^{\mathrm{S}} \delta{xi} \,, \qquad \mu^{\mathrm{D}} = \sum{i=1}^N mi^{\mathrm{D}} \delta{x_i} \,, \label{eq:mu_spectrum}$$

the distance (3.5) is computed by minimizing

$$a \sum{i=1}^{N} \left( m{i}^{\mathrm{S}} - \mathring{m}{i}^{\mathrm{S}} \right) + a \sum{i=1}^{N} \left( m{i}^{\mathrm{D}} - \mathring{m}{i}^{\mathrm{D}} \right) + b \sum{j,k=1}^{N} c{jk} \gamma_{jk} ,$$

where cjk = ∥x^j − xk∥ p , among all γ ∈ [0, +∞) ^N×^N such that

$$\sum{k=1}^{N} \gamma{jk} = \mathring{m}_{j}^{s}, \text{ for } j \in {1, \dots, N},$$

$$\sum{j=1}^{N} \gamma{jk} = \mathring{m}_{k}^{D}, \text{ for } k \in {1, \dots, N},$$

and among all ˚m^s 1 , . . . , ˚m^s ^N , ˚m^d 1 , . . . , ˚m^d ^N with ˚m^s ⁱ ≤ m^s i and ˚m^d ⁱ ≤ m^d i , i ∈ {1, . . . , N}. By taking into account such conditions, the objective is to minimize

$$a\sum_{i=1}^{N} \left(m_i^{\rm S} + mi^{\rm D}\right) + \sum{j,k=1}^{N} \left(b c{jk} - 2a\right) \gamma{jk}$$

among all γ ∈ [0, +∞) ^N×^N such that

$$\sum{k=1}^{N} \gamma{jk} \le m_j^{\mathrm{s}}, \quad \text{for all } j = 1, \dots, N,$$

$$\sum{j=1}^{N} \gamma{jk} \le m_k^{\mathrm{d}}, \quad \text{for all } k = 1, \dots, N.$$

4.2.2 Computational approach

This time, the problem can be solved through a linear programming problem with inequality constraints. First, one solves the problem

minimize $$\hat{\mathbf{c}}^{\dagger}\mathbf{x}$$

subject to $\mathbf{A}\mathbf{x} \leq \mathbf{b}$
$\mathbf{x} \geq 0$ , (4.4)

where ˆc := bc − 2a1N² , x, c, b are defined as in (2.11), and A is defined as in (2.12). Then, once the solution x ^∗ = x ∗ (a, b) of (4.4) is found, the final GWD is given by

$$\tilde{W}1^{\text{\tiny PR}}(m^{\text{\tiny S}},m^{\text{\tiny D}};a,b) := \mathbf{x}^*(a,b) + a\sum{i=1}^N \left(m_i^{\text{\tiny S}} + m_i^{\text{\tiny D}}\right).$$

4.2.3 Numerical tests

Test 1

Consider the simple example proposed in [21, Sect. 2.2], where authors consider the two constant distributions

$$\rho^{\text{S}}(x) = \chi{[-1,0]}(x), \qquad \rho^{\text{D}}(x) = \chi{[\xi,1+\xi]}(x), \qquad x \in \mathbb{R}$$

with parameter ξ ≥ 0. The exact value of the GWD (with a = b = p = 1) is easily computed and gives

$$W_1^{\text{PR}}(\rho^{\text{S}}, \rho^{\text{D}}; 1, 1) = \begin{cases} 1 + \xi - \xi^2 / 4, & \text{if } \xi \in [0, 2], \ 2, & \text{if } \xi \ge 2. \end{cases}$$

Fig. 6 shows the result of the numerical computation in the domain X = [−2, 5] divided in N = 200 cells, for ξ ∈ [0, 3].

Figure 6: Test 1. Exact solution compared with numerical solution (N = 200).

Test 2

In this test we run the algorithm in a nontrivial case where exact solution is not known. More in detail, we consider X = [−4, 4], N = 100, a = b = p = 1, and

$$\rho^{\mathsf{S}}(x) = \frac{e^{-x+1}}{5} \chi_{[-2,0]}(x), \qquad \rho^{\mathsf{D}}(x) = e^{-(x-1)^2}.$$

Fig. 7 shows the result in terms of remaining mass after minimization. One can see that here the

Figure 7: Test 2. Initial distributions m^s and m^d (left) and remaining masses ˚m^s and ˚m^d corresponding to the optimal solution (right).

optimal solution (W˜ pr ¹ = 4.36) requires to destroy only a part of the mass.

4.3 Gaussian Hellinger–Kantorovich approach

4.3.1 Discrete approach

When computing the distance (3.7) of two purely atomic measures with finite supports

$$\mu^{\scriptscriptstyle \mathrm{S}} = \sum_{i=1}^N mi^{\scriptscriptstyle \mathrm{S}} \delta{xi} \,, \qquad \mu^{\scriptscriptstyle \mathrm{D}} = \sum{i=1}^N mi^{\scriptscriptstyle \mathrm{D}} \delta{x_i}$$

we may resort to a dual formulation. Namely, by adapting the proof of [16, Theorem 6.3] to the case in which a, b are not necessarily both equal to 1, one sees that the distance equals the maximum of

$$a \sum{i=1}^{N} m{i}^{\mathrm{S}} \delta{x{i}} \left( 1 - e^{-b\varphi{i}/a} \right) + a \sum{i=1}^{N} m{i}^{\mathrm{D}} \delta{x{i}} \left( 1 - e^{-b\psi{i}/a} \right)$$

among all (φ1, . . . , φ^N , ψ1, . . . , ψ^N ) ∈ R ^N × R ^N such that

$$\varphi_j + \psik \le c{jk} := |x_j - x_k|^2, \quad \text{for all } j, k \in {1, \dots, N}.$$ (4.5)

4.3.2 Computational approach

In this case the optimization problem can be written as a nonlinear maximization problem with linear constraints. To this end, we define the function f : R ²^N → R,

$$f(\mathbf{z}) := a \sum_{i=1}^{N} m_i^{\mathrm{S}} \left( 1 - e^{-bzi/a} \right) + a \sum{i=N+1}^{2N} m_{i-N}^{\mathrm{D}} \left( 1 - e^{-bz_i/a} \right), \qquad \mathbf{z} \in \mathbb{R}^{2N}.$$

Then, we aim at finding

$$\tilde{W}2^{\mathrm{GHK}}(m^{\mathrm{S}},m^{\mathrm{D}}) := \max{\mathbf{z}} f(\mathbf{z})$$

under constraints

$$\mathbf{A}^{\mathsf{T}}\mathbf{z} \leq \mathbf{c},$$

where c is defined as in (2.11) (with cjk as in (4.5)) and A is the N² × 2N matrix defined as

$$\mathbf{A} := \left[ \begin{array}{ccccc} \mathbf{e}_1 \ \mathbf{e}_1 \ \mathbf{e}_1 \ \mathbf{e}_2 \cdots \mathbf{e}_N & \mathbf{e}_2 \cdots \mathbf{e}_N & \mathbf{e}_N \cdots \mathbf{e}_N \ \mathbf{e}_1 \ \mathbf{e}_2 \cdots \mathbf{e}_N & \mathbf{e}_1 \ \mathbf{e}_2 \cdots \mathbf{e}_N & \cdots & \mathbf{e}_1 \ \mathbf{e}_2 & \cdots \mathbf{e}_N \end{array} \right],$$

with eⁱ = (0 0 · · · 1 |{z} i th · · · 0 0)^⊺ ∈ R ^N being, as usual, the ith element of the canonical base of R N .

4.3.3 Numerical test

Numerical tests for the GWD discussed here are in common with those of the following GWD, see Sect. 4.4.3.

4.4 Savar´e & Sodini approach

4.4.1 Discrete approach

For the purpose of empiric approximation, it may be convenient to resort to the convex relaxation approach introduced in [25]. In particular, in view of [25, Th. 3.5], in the case of two atomic homogeneous marginals

$$\mu^{\mathrm{S}} = \sum_{i=1}^N mi^{\mathrm{S}} \delta{xi} \,, \qquad \mu^{\mathrm{D}} = \sum{i=1}^N mi^{\mathrm{D}} \delta{x_i} \,$$

the Gaussian Hellinger–Kantorovich distance can be obtained by minimizing

$$S := \sum{j,h,k,\ell} \gamma{(jh,k\ell)} H_g(xj, r{jh}^{\scriptscriptstyle S}; xk, r{k\ell}^{\scriptscriptstyle D})$$

$$\tag{4.6}$$

under the constraints γ(jh,kℓ) ≥ 0, P jhkℓ γ(jh,kℓ) = 1, and for the homogeneous marginals that

$$m_j^{\scriptscriptstyle \mathrm{S}} = \sumh r{jh}^{\scriptscriptstyle \mathrm{S}} \sum{k\ell} \gamma{(jh,k\ell)} \,, \qquad mk^{\scriptscriptstyle \mathrm{D}} = \sum\ell r{k\ell}^{\scriptscriptstyle \mathrm{D}} \sum{jh} \gamma_{(jh,k\ell)} \,,$$

and

$$\sumh r{jh}^{\scriptscriptstyle \mathrm{S}} = mj^{\scriptscriptstyle \mathrm{S}} \,, \qquad \sum\ell r_{k\ell}^{\scriptscriptstyle \mathrm{D}} = m_k^{\scriptscriptstyle \mathrm{D}} \,,$$

for all j, k ∈ {1, . . . , N} (see also Fig. 8). We point out that this latter approach differs from that of the previous one only numerically, because the distance defined in this way in principle coincides with the one discussed in Sect. 4.3.1.

4.4.2 Computational approach

Starting again from the basic numerical setting introduced in Sect. 2.3, we introduce the additional "vertical" discretization at any fixed grid point xⁱ and for each mass distribution, see Fig. 8.

Figure 8: Additional discretization of the concentrated masses.

At each grid point xⁱ , the mass m^s i is divided in n s ⁱ ≥ 0 intervals (reasonably, n s ⁱ = 0 iff m^s ⁱ = 0). Analogously, the mass m^d i is divided in n d ⁱ ≥ 0 intervals (again, m^d ⁱ = 0 iff m^d ⁱ = 0). Note that n {s,d} i acts as additional discretization parameter, and we expect the numerical error to vanish only for both n {s,d} ⁱ → +∞ and N → +∞.

Now, at each subdivision of the supply mass is associated the additional variable r s jh, j = 1, . . . , N, h = 1, . . . , n^s j , and at each subdivision of the demand mass is associated the additional variable r d kℓ, ^k = 1, . . . , N, ^ℓ = 1, . . . , n^d k .

Define also

$$\Gamma \in \mathbb{R}^d$$ , with $d = \sum_{i=1}^N ni^{\scriptscriptstyle{\mathrm{S}}} \cdot \sum{i=1}^N n_i^{\scriptscriptstyle{\mathrm{D}}}$

as the vector of all the weights γ(jh,kℓ) 's associated to the pairs (r s jh, r^d kℓ), ^j = 1, . . . , N, ^h = 1, . . . , n^s j , k = 1, . . . , N, ℓ = 1, . . . , n^d k .

The optimization problem to be solved reads as

minimize $$S(\Gamma, r^{s}, r^{D})$$

subject to $\gamma{(jh,k\ell)} \geq 0$ , $\forall j, h, k, \ell$ , $\sum{jhk\ell} \gamma{(jh,k\ell)} = 1$ ,
$m{j}^{s} = \sum{h=1}^{n{j}^{s}} \left(r{jh}^{s} \sum{k\ell} \gamma{(jh,k\ell)}\right)$ , $\forall j$ ,
$m{k}^{D} = \sum{\ell=1}^{n{D}^{D}} \left(r{k\ell}^{D} \sum{jh} \gamma_{(jh,k\ell)}\right)$ , $\forall k$ ,

where S is defined in (4.6), and the sums P jhkℓ, P kℓ, P jh run on all over the index combinations. It is important to note that the linear systems for r s , r ^d appearing in (4.7) can be actually underdetermined for a given Γ. Analogously to the previous cases, we will denote the solution of the optimization problem by W˜ ss 2 .

In order to solve the optimization problem (4.7), we have devised two strategies: the first one is based on a slow but exhaustive search, which guarantees to find the global minimum. The second one is based on a nondifferential descend method and allows for faster computations, although a priori limited to the search of a local minimum.

Algorithm SS-A: Exhaustive search

1. Introduce a grid with Q nodes in the unit interval [0, 1], namely 0 = q¹ < q² < . . . , q^Q = 1. Note that Q is yet another approximation parameter (the associated error is expected to vanish for Q → +∞).
2. List all the possible vectors {Γ˜ ∈ R ^d} such that each element belongs to the set {q1, . . . , qQ}, and the sum of the elements equals 1. A good algorithm to do this is the following (with d = 5),

although some generated Γ's lead to linear systems with no solution (we have discarded them). ˜

1: for i1 = 0, . . . , Q do
2: for i2 = 0, . . . , Q − i1 do
3: for i3 = 0, . . . , Q − i1 − i2 do
4: for i4 = 0, . . . , Q − i1 − i2 − i3 do
5: i5 = Q − i1 − i2 − i3 − i4
6: Γ = ( ˜ i1, i2, i3, i4, i5)/Q
7: end for
8: end for
9: end for
10: end for

3. For each Γ, find ˜ {r s jh}jh and {r d kℓ}kℓ (depending on that Γ) solving the underdetermined linear ˜ systems, selecting the solution with minimal norm.
4. For each Γ (and corresponding ˜ {r s jh}jh, {r d kℓ}kℓ), evaluate ^S(Γ). ˜
5. Find minΓ˜ ^S(Γ) by exhaustive comparison. ˜

Algorithm SS-B: Random guess and random descend

List a largest-as-possible number of vectors Γ˜ ∈ R d such that each element belongs to the interval [0, 1] and the sum of the elements equals 1. A simple algorithm to do this is the following.
```
1: for z = 1, . . . , d do
2: Extract a uniform random variable ωz ∈ {1, . . . , Q}
3: end for
4: Γ = ( ˜ ω1, . . . , ωd)/
                  P
                     z ωz
```

2. For each Γ, find ˜ {r s jh}jh and {r d kℓ}kℓ (depending on that Γ) solving the underdetermined linear ˜ systems, selecting the solution with minimal norm.
3. For each Γ (and corresponding ˜ {r s jh}jh, {r d kℓ}kℓ), evaluate ^S(Γ). ˜
4. Find Γ˜^∗ := arg minΓ˜ ^S(Γ) by exhaustive comparison. ˜
5. Starting from Γ˜^∗ as initial guess, explore the control space perturbing the best Γ obtained so far with small perturbation ±ε on random pairs of its elements (thus keeping the sum equal to 1). If the perturbation leads to a lower value, the best Γ^∗ is updated with the perturbed one, otherwise the attempt is discarded. As usual, duly diminishing ε while the algorithm runs can improve the accuracy (simulated annealing).

4.4.3 Numerical tests

We present here numerical tests for both the approaches described in Sects. 4.3.2 and 4.4.2. We recall that the two approaches, although quite different from the computational point of view, are theoretically equivalent (with a = b = 1, which will be always the case from now on). This is also confirmed by our numerical experiments.

Test 1

This test is mainly conceived to check the correctness of the code on simple examples with known exact solution. Let us consider the simplest possible scenario, with only one concentrated supply mass δ ^s and only one concentrated demand mass δ d . To simplify the notations, we denote by m^s , x ^s and by m^d , x d the mass and the position of the supply mass and the demand mass, respectively. In this case the exact solution is W˜ ghk ² = W˜ ss ² = Hg(x s , m^s ; x d , m^d ), where H^g is defined in (3.9).

Fig. 9 shows the results of the algorithm for W˜ ghk as well as both algorithms for W˜ ss (SS-A and SS-B), alongside the exact solution. More in detail, Fig. 9(left) shows the GWD in the case the two distributions are balanced (m^s = m^d = ¯m) and their distance progressively increases from 0 to 6. One can see that the GWD rapidly saturates at 2 ¯m.

Fig. 9(right) shows instead the GWD in the case the two distributions are located at the same point and the difference of their masses progressively increases from 0 to 7.

Figure 9: Test 1. GWD between two equal atomic mass distributions with mass ¯m as a function of their distance ∆ξ (left), and GWD between two overlapping atomic mass distributions as a function of the difference of mass (right). Dotted black line is the exact solution.

Test 2

The goal of this test is mainly comparing the algorithms SS-A and SS-B. We consider another simple scenario with only two concentrated masses for supply, and only one concentrated mass for demand, see Fig. 10. We assume that supply masses are located at x = 0 and x = 1, have masses

Figure 10: Test 2: two concentrated masses for supply mass and one for demand mass.

m = 2 and m = 1, and have n = 2, n = 1 additional divisions, respectively. Demand mass is instead located at x = 3, has mass m = 4, and has n = 4 additional divisions. In this scenario it results d = 12.

Fig. 11(top) shows the history of all the evaluations of W˜ ss ² done by the exhaustive algorithm SS-A, in the same order as they are computed, with Q = 6. One can observe a recursive pattern due to the way Γ's are generated (Step 2). Minimum is already found among the first evaluations and it is W˜ ss ² = 6.45.

Fig. 11(bottom) shows the tentative values of W˜ ss 2 found in the first exploratory phase of the algorithm SS-B (Steps 1-3), followed by the minimum values found during the random descend, with Q = 1,000. The algorithm returns again W˜ ss ² = 6.45. Convergence seems to be well guaranteed and the two algorithms give the same value for W˜ ss 2 (W˜ ghk 2 is also the same).

For these reasons, in the following we will use algorithm SS-B only, which is obviously faster.

Figure 11: Test 2. Algorithm SS-A, complete view of all the evaluations of W˜ ss 2 (top), and Algorithm SS-B, 1,000 random searches followed by about 1,000 descend iterations (bottom). Dotted red line indicates the minimum value found in the first phase of SS-B.

5 Applications to traffic model analysis

In this section we will try to understand to what extent the GWDs previously studied are suitable in the context of traffic flow models. It is useful to recall here that in [4] authors have already extensively studied the sensitivity of the LWR model [17, 23] using the WD. More precisely, many simulations were performed varying some key parameters like the initial vehicle density, the road capacity (maximal flux), the critical density (vehicle density associated to the maximal flux), the size of the road network, and the distribution coefficients at junctions (how traffic distributes across a junction with more than one outgoing road). Interestingly, authors found that the parameter the model is most sensitive to is the last one, namely the distribution coefficients at junctions. That said, the usage of WD have limited the authors to the comparison of balanced vehicle distributions only. This means, in practice, that the number of vehicles should be the same in the two considered scenarios.

In the following we try to complete the sensitivity analysis including the cases in which the mass distributions are unbalanced.

5.1 Mathematical modelling

In a macroscopic setting, traffic flow is described by means of average quantities only, like density ρ (vehicles per unit of space, at a given time), flux f (vehicles per unit of time crossing a given point, at a given time), and velocity v. Models can be first-order (i.e., velocity based), if velocity is given as a function of the density and acceleration is neglected, or second-order (i.e., acceleration based) if, instead, velocity is an independent variable and acceleration is given as a function of both density and velocity.

In the following we will assume the road is a finite segment X = [xmin, xmax] with a single lane and a single class of vehicles.

LWR model

The LWR model, introduced independently by Lighthill & Whitham [17] and by Richards [23], describes the evolution of the vehicle density ρ(x, t) by means of the following IBVP involving an hyberbolic PDE

$$\begin{cases} \partial{t}\rho + \partial{x}f(\rho) = 0, & t \in [t{0}, T], \quad x \in [x{\text{MIN}}, x{\text{MAX}}], \ \rho(x, 0) = \rho{0}(x), & x \in [x{\text{MIN}}, x{\text{MAX}}], \ \rho(x{\text{MIN}}, t) = \rho{\text{IN}}(t), & t \in [t{0}, T], \ \rho(x{\text{MAX}}, t) = \rho{\text{OUT}}(t), & t \in [t{0}, T], \end{cases}$$ (5.1)

for some initial time t0, final time T > t0, and concave flux function f : [0, ρmax] → R ⁺ (called fundamental diagram in the context of traffic flow) such that f(0) = f(ρmax) = 0, for some maximal density ρmax. Note that we have naturally guaranteed the constraint 0 ≤ ρ ≤ ρmax at any time, provided 0 ≤ ρ0, ρin, ρout ≤ ρmax.

Recalling the well known physical law f = ρv, which is always valid for any x and t, we get that the LWR model implicitly assumes that the vehicle velocity v only depends on the density ρ (since f = f(ρ), therefore v = v(ρ)). From the modeling point of view, this means that the velocity adapts instantaneously to the changes of the density, i.e. acceleration is infinite. For this reason the LWR model is said to be first-order.

As usual in the mathematical literature, in the following we will define

$$f(\rho) = 1 - \rho,$$

implicitly assuming that the maximal density ρmax and the maximal velocity of vehicles are normalized and equal to 1. This is done for simplicity, since a qualitative analysis is sufficient to highlight the differences among the GWDs.

Eq. (5.1) can be numerically approximated by means of, e.g., the Godunov scheme. Using again the numerical setting introduced in Sect. 2.3, this scheme reads as

$$\rho_i^{n+1} = \rho_i^n - \frac{\Delta x}{\Delta t} \Big( \mathcal{G}(\rhoi^n, \rho{i+1}^n) - \mathcal{G}(\rho_{i-1}^n, \rho_i^n) \Big), \qquad i = 2, \dots, N-1, \quad n = 0, \dots, \left\lceil \frac{T - t_0}{\Delta t} \right\rceil, \quad (5.2)$$

where ρ n i represents the approximated value of ρ at the center of the cell Cⁱ at time t⁰ + n∆t, and G is the Godunov's numerical flux defined by

$$\mathcal{G}(\rho{-}, \rho{+}) := \begin{cases} \min{z \in [\rho{-}, \rho{+}]} f(z), & \text{if } \rho{-} \leq \rho{+}, \ \max{z \in [\rho{+}, \rho{-}]} f(z), & \text{if } \rho{-} \geq \rho{+}. \end{cases}$$ (5.3)

Remark 6. Boundary conditions can be also given in form of incoming and outgoing flux, rather than as left and right densities ρin, ρout. This happens, e.g., if boundary data are given by fixed sensors counting vehicles passing through x = xmin and x = xmax in a time interval ∆t. If this is the case, values ρ n ¹ = ρ n in and ρ n ^N = ρ n out are dropped, and incoming/outgoing fluxes F n in, F n out are directly injected in the scheme (5.2) taking the place of G(ρ n i−1 , ρⁿ i ) if i = 1, and of G(ρ n i , ρⁿ ⁱ+1) if i = N, respectively.

ARZ model

The ARZ model, introduced independently by Aw & Rascle [2] and by Zhang [29], describes the evolution of the vehicle density ρ(x, t) and the velocity v(x, t) by means of the following IBVP involving two hyberbolic PDEs

$$\begin{cases} \partial{t}\rho + \partial{x}(\rho v) = 0, & t \in [t{0}, T], \quad x \in [x{\text{MIN}}, x{\text{MAX}}], \ \partial{t}(\rho w) + \partial{x}(v\rho w) = \rho \frac{v{eq}(\rho) - v}{\tau}, & t \in [t{0}, T], \quad x \in [x{\text{MIN}}, x{\text{MAX}}], \ \rho(x, 0) = \rho{0}(x), & x \in [x{\text{MIN}}, x{\text{MAX}}], \ v(x, 0) = v{0}(x), & x \in [x{\text{MIN}}, x{\text{MAX}}], \ \rho(x{\text{MIN}}, t) = \rho{\text{IN}}(t), & t \in [t{0}, T], \ \rho(x{\text{MAX}}, t) = \rho{\text{OUT}}(t), & t \in [t{0}, T], \ v(x{\text{MIN}}, t) = v{\text{IN}}(t), & t \in [t{0}, T], \ v(x{\text{MAX}}, t) = v{\text{OUT}}(t), & t \in [t_{0}, T], \end{cases}$$

$$(5.4)$$

where w(x, t) := v(x, t) + P(ρ(x, t)), P(ρ) := ^vref γ ρ ^ρmax ^γ , and vref , γ, an τ are strictly positive parameters. This time, the velocity v is an independent unknown and the acceleration ^veq(ρ)−^v τ appears explicitly in the model, depending on ρ and v. For this reason, the ARZ model is said to be second-order.

In order to fairly compare the two models we define veq(ρ) := ^f(ρ) ρ , so that the desired equilibrium velocity of the ARZ model coincides with the instantaneous velocity of the LWR model. Also in the boundary conditions, we set vin = veq(ρin) and vout = veq(ρout).

Eq. (5.4) is numerically approximated by means of the Lax–Friedrichs scheme.

5.2 Numerical tests

5.2.1 Test 0: comparison of the three GWDs

In this preliminary test, we compare the four computational approaches for the three GWDs individually studied above, in order to understand which one is the most suitable in the context of traffic flow modeling. This is done by devising a numerical test which can be run by all algorithms. From now on, the four algorithms will be referred to by FG, PR, SS (variant B only), and GHK, with obvious meaning.

Denoting by χ the indicator function, we consider the case of two time-dependent distributions with parameters α, β, η > 0,

$$\rho^{S}(x,t) = \alpha \chi{[-t-\eta, -t+\eta]}(x), \qquad \rho^{D}(x,t) = \beta \chi{[t-\eta, t+\eta]}(x), \qquad 0 \le t \le T$$ (5.5)

in the domain X = [−T −η, T +η]. The two step functions start with perfectly overlapping support, then supply mass moves leftward while demand mass moves rightward, until they both reach the boundary of the domain, see Fig. 12.

Figure 12: Moving distributions ρ ^s and ρ ^d as defined in (5.5).

In particular, for this comparison test we set α ∈ {1, 2}, β = 2, η = 1, T = 4. Moreover, we have the following algorithm-specific parameters:

• PR: a = 2.5, b = 1.
GHK: a = b = 1.
SS: a = b = 1, 4 atomic masses per distribution, 4 additional divisions per mass, Q = 1,000, ϵ = 0.001, iterations 5,000 (random searches) + 40,000 (random descend), average on 10 runs.

Fig. 13 shows the behavior of the three GWDs between ρ ^s and ρ ^d as defined in (5.5), as a function of time. WD is also shown in the balanced case α = β.

Figure 13: Comparison of the three GWDs between ρ ^s and ρ ^d as defined in (5.5), as a function of time t (β = 2, η = 1, T = 4). α = 2 (balanced), including also classical WD (left), and α = 1 (unbalanced, right).

This comparison test allows us to sketch some preliminary conclusions. First of all, SS seems to be computationally unfeasible. It is impossible to deal with a reasonable number of concentrated masses and additional divisions, therefore it will be discarded in the following analysis.

Among the others, PR and GHK share with L ^p an important drawback: the distances saturate when distributions do not overlap and they are far enough from each other. This does not seem to be a nice feature in the context of traffic flow models.

Lastly, FG has an important feature which deserves a discussion: the distance decreases whenever distributions are close to the boundaries, even if they are very far from each other. On first glance, this seems to be an issue, but it can be also be seen as a good point: basically, it implies that what happens near the boundaries is less important than what happens on the rest of the road. Considering the fact that boundary conditions are typically unknown in traffic flow modeling, it could make sense to give priority to that part of the road which is less or not at all influenced by boundary conditions.

5.2.2 Test 1: sensitivity to boundary conditions

In this test we investigate the sensitivity to boundary conditions. We consider the case described in Rem. 6, which is a scenario often met in real life [5]. In the same spirit of [5], we assume that the time interval is divided in two subintervals: [t0, 0] and (0, T], where t = 0 corresponds to the current time (the 'now'), i.e. the time the simulation is actually performed. Until time 0, we can rely on the sensors' measurements Fin, Fout, and we use it as boundary conditions (see Rem. 6). After that time, we have no such a data because they are not yet available, but we assume to be able to (exactly) forecast the average value of them,

$$\hat{F}_{\text{in}} = \frac{1}{T} \int0^T F{\text{in}}(t) dt, \qquad \hat{F}_{\text{out}} = \frac{1}{T} \int0^T F{\text{out}}(t) dt.$$

In short: we denote by ρ s the "exact" solution obtained by using boundary fluxes Fin(t), Fout(t) in the whole interval [t0, T], and by ρ d the "predicted" solution obtained by using flux boundary data Fin(t), Fout(t) until time t = 0 and then by using the predicted average values Fˆ in, Fˆ out until time T. At any time t, the vehicle distributions ρ s (·, t) and ρ d (·, t) are compared computing the GWDs.

Initial conditions are ρ s (x, t0) = ρ d (x, t0) = 0.2 for any x. Boundary fluxes are randomly generated: Fin has uniform distribution in [0, 4 ], Fout has uniform distribution in [0, G(0.2, 0.2)]. Other parameters are: xmin = 0, xmax = 40, t⁰ = −10, T = 30, ∆x = 0.4 (corresponding to N = 100), and ∆t = 0.25.

Fig. 14 shows two screenshots of the simulations at times t = 0 and t = 12.

Figure 14: Test 1 on sensitivity to boundary conditions. Two screenshots at t = 0 ('now'), i.e. when solutions cease to be equal (left), and t = 12 (right).

Fig. 15 shows the total mass on the road and the GWDs between ρ s (·, t) and ρ d (·, t) at any time.

5.2.3 Test 2: sensitivity to traffic light cycle

In this test we consider a shorter road with respect to the previous test, and we assume there is a traffic light in the middle of the road, at the interface between two numerical cells. The duration

Figure 15: Test 1 on sensitivity to boundary conditions. Difference of mass between the two simulations (left) and comparison of GWDs (right).

of the red phase equals that of the green phase. Technically, the red traffic light is modeled by imposing null flux at the right (resp., left) boundary of the cell before (resp., after) the traffic light. The aim is to investigate the sensitivity to the duration of the traffic light cycle.

Initial conditions are ρ s (x, t0) = ρ d (x, t0) = 0 for any x. Left Dirichlet boundary conditions are ρ s in = ρ d in = 0.4. Right Dirichlet boundary conditions are ρ s out = ρ d out = 0. Other parameters are: xmin = 0, xmax = 4, t⁰ = 0, T = 20, ∆x = 0.04 (corresponding to N = 100), and ∆t = 0.025. Both red and green phases last 50 time steps for supply distribution and 40 time steps for demand distribution.

Fig. 16 shows two screenshots of the simulations at times t = 3.7 and t = 7.4.

Fig. 17 shows the total mass on the road in the two scenarios and the GWDs between ρ s (·, t) and ρ d (·, t) at any time.

5.2.4 Test 3: sensitivity to model order

In this test we quantify the difference between two simulations obtained by means of the LWR model (supply, first-order) and the ARZ model (demand, second-order).

Initial conditions are ρ s (x, t0) = ρ d (x, t0) = 0.7 for x ∈ [1.6, 2.4] and 0.1 elsewhere. Left and right Dirichlet boundary conditions are null. Other parameters are: τ = 0.05, γ = 2, vref = 1, xmin = 0, xmax = 4, t⁰ = 0, T = 3, ∆x = 0.02 (corresponding to N = 200), and ∆t = 0.01.

Fig. 18 shows two screenshots of the simulations at times t = 0.3 and t = 1.

Fig. 19 shows the total mass on the road in the two scenarios and the GWDs between ρ s (·, t) and ρ d (·, t) at any time. Here we have chosen a = 0.5 and b = 1 as parameters for computing W˜ pr 1

Figure 16: Test 2 on sensitivity to traffic light cycle. Two screenshots at t = 3.7 with red traffic light for demand and green traffic light for supply (left), and t = 7.4 with red traffic light for both distributions (right).

Figure 17: Test 2 on sensitivity to traffic light cycle. Difference of mass between the two simulations (correctly masses differ only in a certain time interval, left), and comparison of GWDs (right).

6 Discussion

Observing the results of the last numerical tests, we tend to confirm what we have anticipated in Sect. 5.2.1, namely that the Figalli & Gigli approach seems to be the most suitable GWD to be

Figure 18: Test 3 on sensitivity to model order. Two screenshots at t = 0.3 (left), and t = 1 (right).

Figure 19: Test 3 on sensitivity to model order. Difference of mass between the two simulations (left) and comparison of GWDs (right).

used in the context of traffic flow modeling, mainly because of how boundaries are regarded.

In this respect, the most significant test is probably the first one, in Sect. 5.2.2: here the discrepancy between the exact and predicted solution slowly propagates from the boundaries to the interior of the road, impacting more and more on the reliability of the traffic prediction. This behavior of the error is clearly better caught by FG (see Fig. 17-right), which is the last distance that saturates.

Although less evident, the same conclusions apply for the third test, in Sect. 5.2.4: here the two solutions start differing immediately, both at the boundary and inside the road. Moreover, the difference increases as the time goes by. Among the three GWDs, FG is the one which increases more rapidly at the initial time, and then keeps increasing all the time.

Conversely, the second test does not show much difference among the three GWDs, having them basically the same cycling behavior. That behavior is correct since the two solutions cyclically move away from each other and then come closer again.

Overall, we have proposed four algorithms for computing, in a discrete context, four GWDs (two of them theoretically equivalent to each other). One of them (SS) results to be numerically unfeasible, while the other three (namely FG, PR, GHK) are all relative simple to implement once the computational procedure is arranged.

In a future work, we will aim at exploring the possibility to compute the considered GWDs in more complicated domain like R 2 , R 3 , and (road) networks.

Funding. M.B. and E.C. would like to thank the Italian Ministry of University and Research (MUR) to support this research with funds coming from PRIN Project 2022 PNRR entitled "Heterogeneity on the road - Modeling, analysis, control", No. 2022XJ9SX.

E.C. would like to thank the Italian Ministry of University and Research (MUR) to support this research with funds coming from PRIN Project 2022 entitled "Optimal control problems: analysis, approximation and applications", No. 2022238YY5.

G.F. would like to thank the Italian Ministry of University and Research (MUR) to support this research with funds coming from PRIN Project 2022 entitled "Variational analysis of complex systems in materials science, physics and biology", No. 2022HKBF5C, CUP B53D23009290006.

This study was carried out within the Spoke 7 of the MOST – Sustainable Mobility National Research Center and received funding from the European Union Next-Generation EU (PIANO NAZIONALE DI RIPRESA E RESILIENZA (PNRR) – MISSIONE 4 COMPONENTE 2, INVES-TIMENTO 1.4 – D.D. 1033 17/06/2022, CN00000023).

This manuscript reflects only the authors' views and opinions. Neither the European Union nor the European Commission can be considered responsible for them.

E.C. and F.L.I. are funded by INdAM–GNCS Project, CUP E53C23001670001, entitled "Numerical modeling and high-performance computing approaches for multiscale models of complex systems".

M.B., E.C., and F.L.I. are members of the e INdAM research group GNCS, while G.F. is member of the e INdAM research group GNAMPA.

Acknowledgments. Authors want to thank Benedetto Piccoli, Francesco Rossi, Giuseppe Savar´e, and Giacomo Enrico Sodini for the useful discussions, and Anna Th¨unen for the contribution given to Sect. 4.1.2.

References

[1] L. Ambrosio, Mainini E., and Serfaty S. Gradient flow of the chapman-rubinstein-schatzman model for signed vortices. Ann. Inst. H. Poincar´e C Anal. Non Lin´eaire, 28(2), 2011.
[2] A. Aw and M. Rascle. Resurrection of "second order" models of traffic flow? SIAM J. Appl. Math., 60:916–938, 2000.
[3] J.-D. Benamou and Y. Brenier. A computational fluid mechanics solutions to the Monge– Kantoriovich mass transfer problem. Numer. Math., 84(3):375–393, 2000.
[4] M. Briani, E. Cristiani, and E. Iacomini. Sensitivity analysis of the LWR model for traffic forecast on large networks using Wasserstein distance. Commun. Math. Sci., 16(1):123–144, 2018.
[5] M. Briani, E. Cristiani, and E. Onofri. Inverting the fundamental diagram and forecasting boundary conditions: How Machine Learning can improve macroscopic models for traffic flow. To appear on Advances in Computational Mathematics. Preprint arXiv:2303.12740.
[6] L. Caffarelli and R. J. McCann. Free boundaries in optimal transport and Monge-Amp`ere obstacle problems. Annals of Mathematics, 171(2):673–730, 2010.
[7] L. Chizat, G. Peyr´e, B. Schmitzer, and F.-X. Vialard. An interpolating distance between optimal transport and Fisher-Rao metrics. Foundation of Computational Mathematics, 18:1– 44, 2018.
[8] L. Chizat, G. Peyr´e, B. Schmitzer, and F.-X. Vialard. Unbalanced optimal transport: Dynamic and Kantorovich formulations. Journal of Functional Analysis, 274:3090–3123, 2018.
[9] E. Cristiani, B. Piccoli, and A. Tosin. Multiscale modeling of pedestrian dynamics, volume 12 of Series Modeling, Simulation & Applications. Springer, 2014.
[10] E. Cristiani and M. C. Saladino. Comparing comparisons between vehicular traffic states in microscopic and macroscopic first-order models. Mathematical Methods in the Applied Sciences, 42(3):918–934, 2019.
[11] A. Figalli. The optimal partial transport problem. Archive for Rational Mechanics and Analysis, 195(2):533–560, 2010.
[12] A. Figalli and N. Gigli. A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions. Journal de Math´ematiques Pures et Appliqu´ees, 94(2):107–130, 2010.
[13] L. G. Hanin. Kantorovich-Rubinstein norm and its application in the theory of Lipschitz spaces. Proceedings of the American Mathematical Society, 115(2):345–352, 1992.
[14] L. Kantorovich and G. S. Rubinstein. On a space of totally additive functions. Vestnik Leningradskogo Universiteta, 13:52–59, 1958.
[15] S. Kondratyev, L. Monsaingeon, and D. Vortnikov. A new optimal transport distance on the space of finite Radon measures. Adv. Differential Equations, 21(11/12):1117–1164, 2016.
[16] M. Liero, A. Mielke, and G. Savar´e. Optimal entropy-transport problems and a new Hellinger-Kantorovich distance between positive measures. Invent. Math., 211(3):969–1117, 2018.
[17] M. J. Lighthill and G. B. Whitham. On kinematic waves II. A theory of traffic flow on long crowded roads. Proc. R. Soc. Lond. Ser. A, 229:317–345, 1955.
[18] L. Lombardini and F. Rossi. Obstructions to extension of Wasserstein distances for variable masses. Proceedings of the American Mathematical Society, 150(11):4879–4890, 2022.
[19] Z. Ma, X. Wei, X. Hong, H. Lin, Y. Qiu, and Y. Gong. Learning to count via unbalanced optimal transport. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 2319–2327, 2021.
[20] E. Mainini. A description of transport cost for signed measures. J. Math. Sci. (N.Y.), 181(6), 2012.
[21] B. Piccoli and F. Rossi. Generalized Wasserstein distance and its application to transport equations with source. Archive for Rational Mechanics and Analysis, 211:335–358, 2014.
[22] B. Piccoli and F. Rossi. On properties of the generalized Wasserstein distance. Arch. Rational Mech. Anal., 222:1339–1365, 2016.
[23] P. I. Richards. Shock waves on the highway. Oper. Res., 4:42–51, 1956.
[24] F. Santambrogio. Optimal Transport for Applied Mathematicians. Calculus of Variations, PDEs, and Modeling, volume 87 of Series Progress in Nonlinear Differential Equations and Their Applications. Birkh¨auser, 2015.
[25] G. Savar´e and G. E. Sodini. A relaxation viewpoint to Unbalanced Optimal Transport: Duality, optimality and Monge formulation. Journal de Math´ematiques Pures et Appliqu´ees, 188:114– 178, 2024.
[26] S. M. Sinha. Mathematical Programming. Theory and Methods. Elsevier, 2006.
[27] C. Villani. Topics in Optimal Transportation, volume 58 of Series Graduate Studies in Mathematics. American Mathematical Society, 2003.
[28] C. Villani. Optimal Transport. Old and New, volume 338 of Series Grundlehren der mathematischen Wissenschaften. Springer, 2009.

| [29] | | H. M. Zhang. A non-equilibrium traffic model devoid of gas-like behavior. | | Transportation Res. | | |------|---------------------------|---------------------------------------------------------------------------|--|---------------------|--| | | Part B, 36:275–290, 2002. | | | | |

brianinumericalcomputationgeneralizedwasserstein2024.md 94 KB Историја Датотека

status: converted

Numerical Computation of Generalized Wasserstein Distances with Applications to Traffic Model Analysis

Abstract

Contents

1 Introduction

2 The Wasserstein distance in the context of vehicular traffic

2.1 Motivations

2.2 Basic theory

2.3 Numerical setting

3 Basic theory for Generalized Wasserstein distances

3.1 Figalli & Gigli approach [12]

3.2 Piccoli & Rossi approach [21]

3.3 Gaussian Hellinger–Kantorovich approach [16]

3.4 Savar´e & Sodini approach [25]

4 Discrete and computational approaches for Generalized Wasserstein distances

4.1 Figalli & Gigli approach

4.1.1 Discrete approach

4.1.2 Computational approach

4.1.3 Numerical tests

Test 1

Test 2

4.2 Piccoli & Rossi approach

4.2.1 Discrete approach

4.2.2 Computational approach

4.2.3 Numerical tests

Test 1

Test 2

4.3 Gaussian Hellinger–Kantorovich approach

4.3.1 Discrete approach

4.3.2 Computational approach

4.3.3 Numerical test

4.4 Savar´e & Sodini approach

4.4.1 Discrete approach

4.4.2 Computational approach

Algorithm SS-A: Exhaustive search

Algorithm SS-B: Random guess and random descend

4.4.3 Numerical tests

Test 1

Test 2

5 Applications to traffic model analysis

5.1 Mathematical modelling

LWR model

ARZ model

5.2 Numerical tests

5.2.1 Test 0: comparison of the three GWDs

5.2.2 Test 1: sensitivity to boundary conditions

5.2.3 Test 2: sensitivity to traffic light cycle

5.2.4 Test 3: sensitivity to model order

6 Discussion

References

brianinumericalcomputationgeneralizedwasserstein2024.md 94 KB

Историја Датотека