2.1. Gas flow in single pipes

In a given time interval $\left [0, T\right ]$ with $T \gt 0$ , on each edge $e \in \mathcal {E}$ , the physical properties of gas flow can be modelled by the isothermal Euler equations for compressible and inviscid fluids. In order to simplify the Euler equation system, we assume the pipe walls and the gas to have the same temperature $\mathcal {T}$ , which allows us to omit the energy conservation equation (see [Reference Domschke, Hiller, Lang, Mehrmann, Morandin and Tischendorf11]). For a mass density $\rho _e: \left [0, L_e\right ] \times \left [0, T\right ] \longrightarrow \mathbb {R}_{\geq 0}$ and a velocity $v_e: \left [0, L_e\right ] \times \left [0, T\right ] \longrightarrow \mathbb {R}$ , the system of equations reads as

(ISO1) \begin{align} \begin{split} \frac {\partial \rho _e}{\partial t} + \frac {\partial }{\partial x}\left (\rho _e v_e\right ) & = 0,\\ \frac {\partial (\rho _e v_e)}{\partial t} + \frac {\partial }{\partial x} \left (p_e \left (\rho _e\right )+ \rho _e v_e^2 \right ) & = -\frac {\lambda _e}{2\mathcal {D}_e} \rho _e v_e \left \lvert v_e\right \rvert - g \rho _e \sin {\left (\omega _e\right )},\\ \text {for} \, \left (x, t\right ) & \in \left (0, L_e\right ) \times \left (0, T\right ). \end{split} \end{align}

The first equation (continuity equation) encodes conservation of mass and the second equation ensures conservation of momentum. Here, $\lambda _e \geq 0$ denotes the pipe friction coefficient, $\mathcal {D}_e \gt 0$ the pipe diameter, $g \approx 6.7 \cdot 10^{-11} \, N m^2 / kg^2$ the gravitational constant and $\omega _e \in \left [0, 2\pi \right ]$ the inclination level of the pipe $e$ . The pressure $p_e: \left [0, L_e\right ] \times \left [0, T\right ] \longrightarrow \mathbb {R}_{\geq 0}$ is implicitly defined by the state of real gases equation

(RGE) \begin{equation} p_e\left (\rho _e\right ) = R \rho _e \mathcal {T} z, \end{equation}

with gas constant $R \approx 8.3 \, J / \left (K \cdot mol\right )$ , temperature $\mathcal {T} \gt 0$ of the gas and pipe walls and compressibility factor $z \geq 0$ , which we assume to be constant, but generally depends on $p_e$ and in the non-isothermal case also on $\mathcal {T}$ . The state of real gases equation constitutes the pressure law for ideal gases in the case of $z = 1$ . For each edge $e \in \mathcal {E}$ , we furthermore define the mass flux

\begin{equation*} j_e: \left [0, L_e\right ] \times \left [0, T\right ] \longrightarrow \mathbb {R} \quad \left (x, t\right ) \mapsto j_e\left (x, t\right ) :\!= \rho _e\left (x, t\right ) v_e\left (x, t\right ). \end{equation*}

In order to further simplify the (ISO1) model, we assume small flow rates $|v_e|$ , which enables us to eliminate the non-linearity on the left side of the momentum equation (leading to the (ISO2) model in [Reference Domschke, Hiller, Lang, Mehrmann, Morandin and Tischendorf11]), as shown in [Reference Osiadacz26]. Furthermore, we assume a friction-dominated flow, which corresponds to the friction on the pipe walls dominating the heat conduction effects. If we also neglect the gravitational force in the asymptotic considerations, then we obtain, by using the scaling approaches of [Reference Brouwer, Gasser and Herty7], an even further simplified left side of the momentum equation and overall the model then corresponds to the (ISO3) model in [Reference Domschke, Hiller, Lang, Mehrmann, Morandin and Tischendorf11], which is given by

(ISO3) \begin{align} \begin{split} \frac {\partial \rho _e}{\partial t} + \frac {\partial }{\partial x}\left (\rho _e v_e\right ) & = 0,\\ \frac {\partial }{\partial x} \left (p_e \left (\rho _e\right )\right ) & = -\frac {\lambda _e}{2\mathcal {D}_e} \rho _e v_e \left \lvert v_e\right \rvert - g \rho _e \sin {\left (\omega _e\right )}, \end{split} \end{align}

with the linear constitutive relation of $p_e(\rho _e)$ being based on (RGE).

2.2. Gas flow at interior vertices

For connecting the gas flow of the individual pipes to a feasible network, we need coupling conditions for the flux at vertices, which model pipe intersections. This applies to both interior vertices $\stackrel{\small\circ}{\mathcal {V}}$ and boundary vertices $\partial \mathcal {V}$ , where for interior vertices $\nu \in \stackrel{\small\circ}{\mathcal {V}}$ we additionally assume that it is possible to store gas at the individual vertices.

In order to formulate different coupling and boundary conditions in a unified manner, we introduce the vertex excess flux, based on a generalized version of Kirchhoff’s law, as

(VF KL) \begin{equation} f_{\nu }: \left [0, T\right ] \longrightarrow \mathbb {R} \quad t \mapsto f_{\nu }\left (t\right ) :\!= \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} j_e\left (L_e, t\right ) - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} j_e\left (0, t\right ). \end{equation}

For each vertex $\nu \in \mathcal {V}$ , it catches the difference between the total influx from all ingoing edges and the total outflux from all outgoing edges.

To model the possibility to store gas at interior vertices $\nu \in \stackrel{\small\circ}{\mathcal {V}}$ , we introduce vertex mass densities $\gamma _{\nu }: \left [0, T\right ] \longrightarrow \mathbb {R}_{\geq 0}$ for all interior vertices $\nu \in \stackrel{\small\circ}{\mathcal {V}}$ . At each interior vertex, we have the initial storage volume $\gamma _{\nu }\left (0\right ) \geq 0$ and for all $t \in \left (0, T\right )$ the vertex mass density has to satisfy

(1) \begin{equation} \frac {\partial \gamma _{\nu }}{\partial t} = f_{\nu }. \end{equation}

The combination of the vertex excess flux $f_{\nu }$ and the vertex mass density $\gamma _{\nu }$ ensures mass conservation at each interior vertex at all times, since gas can either flow through the vertex completely (indicating $f_{\nu } = 0$ , as accumulated inflow equals accumulated outflow) or gas is stored at the given vertex for $f_{\nu } \gt 0$ or gas is taken out of storage at the given vertex for $f_{\nu } \lt 0$ .

Not storing gas in the interior vertices can be achieved by setting $f_{\nu } \equiv 0$ . Hence, for all $t \in \left [0, T\right ]$ we obtain the classical version of Kirchhoff’s law

(C KL) \begin{equation} 0 = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} j_e\left (L_e, t\right ) - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} j_e\left (0, t\right ). \end{equation}

Observation 2.4 (No storage at interior vertices). After setting $f_{\nu } \equiv 0$ , the vertex mass density $\gamma _{\nu }$ can be omitted from further calculations, since $\frac {\partial \gamma _{\nu }}{\partial t} \equiv 0$ directly implies a constant amount of stored gas in each vertex $\nu \in \stackrel{\small\circ}{\mathcal {V}}$

\begin{equation*} \gamma _{\nu }\left (0\right ) = \gamma _{\nu }\left (T\right ) \equiv \gamma _{\nu }\left (t\right ), \qquad \forall t \in \left (0, T\right ). \end{equation*}

Note that, continuity of the velocity function $v$ or the mass density $\rho$ in the interior vertices is not covered by these coupling conditions, however it is also not expected in general.

2.3. Gas flow at boundary vertices

Boundary vertices $\nu \in \partial \mathcal {V}$ always either constitute source vertices or sink vertices of the network, which are responsible for the supply (inflow) and demand (outflow) of gas in the network. In contrast to interior vertices, we assume that storage of gas at boundary vertices is not possible. Therefore, we need a modified version of Kirchhoff’s law, which also includes the inflow and outflow of gas to and from outside of the network, to ensure conservation of mass at boundary vertices. Hence, we introduce time-dependent and time-independent boundary conditions for boundary vertices $\nu \in \partial \mathcal {V}$ .

Note that, in PDE literature this type of boundary conditions is typically called “inhomogeneous” and inspired by the coupling conditions at interior vertices, we can formulate these conditions as a system of equations.

Observation 2.6 (Time-dependent boundary conditions). With the generalized Kirchhoff’s law (VF KL), these boundary conditions can be written as

\begin{align*} \qquad s_{\nu }^G \left (t\right ) & = f_{\nu } \left (t\right ) \qquad & \forall \nu \in \partial ^+ \mathcal {V}, \, t \in \left [0, T\right ], \qquad \\ \qquad d_{\nu }^G \left (t\right ) & = f_{\nu } \left (t\right ) \qquad & \forall \nu \in \partial ^- \mathcal {V}, \, t \in \left [0, T\right ], \qquad \end{align*}

because the source vertex flux $\left \lvert s_{\nu }^G \right \rvert$ acts as a flow into vertex $\nu \in \partial ^+ \mathcal {V}$ (thus an outflow $\left .j_e\right |_{x = L_e}$ of an imaginary edge $e$ ) and the sink vertex flux $d_{\nu }^G$ can be seen as a flow out of vertex $\nu \in \partial ^- \mathcal {V}$ (thus an inflow $\left .j_e\right |_{x = 0}$ of an imaginary edge $e$ ). With this interpretation, we obtain

\begin{align*} \qquad 0 & = \Bigg ( \left \lvert s_{\nu }^G \left (t\right ) \right \rvert + \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} j_e\left (L_e, t\right ) \Bigg ) - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} j_e\left (0, t\right ) \qquad & \forall \nu \in \partial ^+ \mathcal {V}, \, t \in \left [0, T\right ], \qquad & \\ \qquad 0 & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} j_e\left (L_e, t\right ) - \Bigg ( d_{\nu }^G \left (t\right ) + \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} j_e\left (0, t\right ) \Bigg ) \qquad & \forall \nu \in \partial ^- \mathcal {V}, \, t \in \left [0, T\right ]. \qquad & \\ \end{align*}

Another option for conditions at the boundary vertices $\partial \mathcal {V}$ is only considering the (mandatory) total gas supply $\left \lvert S_{\nu }^G \right \rvert$ for each source vertex $\nu \in \partial ^+ \mathcal {V}$ for the time period $\left [0, T\right ]$ , which is available starting at time $t = 0$ , as well as the (mandatory) total gas demand $D_{\nu }^G$ for each sink vertex $\nu \in \partial ^- \mathcal {V}$ for the time period $\left [0, T\right ]$ , which has to be met at the latest at time $t = T$ . Since these conditions only have to be fulfilled at the end of the time period, we call them time-independent boundary conditions.

Definition 2.7 (Time-independent boundary conditions). Suppose for the time period $\left [0, T\right ]$ , the accumulated supplies $\left \lvert S_{\nu }^G \right \rvert$ (with $S_{\nu }^G \leq 0$ ) at source vertices $\nu \in \partial ^+ \mathcal {V}$ and accumulated demands $D_{\nu }^G \geq 0$ at sink vertices $\nu \in \partial ^- \mathcal {V}$ are given. Then our aim is to find source vertex flux functions

\begin{equation*} s_{\nu }: \left [0, T\right ] \longrightarrow \mathbb {R}_{\leq 0} \quad t \mapsto s_{\nu }\left (t\right ) :\!= f_{\nu } \left (t\right ) \end{equation*}

for each source vertex $\nu \in \partial ^+ \mathcal {V}$ , such that

(2) \begin{equation} S_{\nu }^G = \int _0^T s_{\nu }\left (t\right ) \, \mathrm d t, \end{equation}

and to find sink vertex flux functions

\begin{equation*} d_{\nu }: \left [0, T\right ] \longrightarrow \mathbb {R}_{\geq 0} \quad t \mapsto d_{\nu }\left (t\right ) :\!= f_{\nu } \left (t\right ) \end{equation*}

for each sink vertex $\nu \in \partial ^- \mathcal {V}$ , such that

(3) \begin{equation} D_{\nu }^G = \int _0^T d_{\nu }\left (t\right ) \, \mathrm d t. \end{equation}

Inspired by the coupling conditions at interior vertices, we can once again reformulate these conditions.

Observation 2.8 (Time-independent boundary conditions). If we define a source vertex mass density $S_{\nu }: \left [0, T\right ] \longrightarrow \mathbb {R}_{\geq 0}$ for all source vertices $\nu \in \partial ^+ \mathcal {V}$ , then constraint ( 2 ), can be rewritten as

\begin{align*} \qquad \frac {\partial S_{\nu } \left (t\right )}{\partial t} & = s_{\nu } \left (t\right ) \qquad & \forall \nu \in \partial ^+ \mathcal {V}, \, t \in \left (0, T\right ), \qquad & \end{align*}

with initial and final condition $S_{\nu }\left (0\right ) = \left \lvert S_{\nu }^G \right \rvert$ and $S_{\nu }\left (T\right ) = 0$ .

In a similar manner, by introducing a sink vertex mass density $D_{\nu }: \left [0, T\right ] \longrightarrow \mathbb {R}_{\geq 0}$ for each sink vertex $\nu \in \partial ^- \mathcal {V}$ , constraint ( 3 ) can be rewritten as

\begin{align*} \qquad \frac {\partial D_{\nu } \left (t\right )}{\partial t} & = d_{\nu } \left (t\right ) \qquad & \forall \nu \in \partial ^- \mathcal {V}, \, t \in \left (0, T\right ), \qquad & \end{align*}

with initial and final condition $D_{\nu }\left (0\right ) = 0$ and $D_{\nu }\left (T\right ) = D_{\nu }^G$ .

Note that in the case of time-independent boundary conditions the non-positivity property of the source vertex flux and the non-negativity of the sink vertex flux are not automatically fulfilled.

The source vertex mass can also be defined in the case of time-dependent boundary conditions as

\begin{equation*} S_{\nu } \left (t\right ) :\!= \big \lvert {\underbrace {S_{\nu }^G}_{\leq \, 0}} \vphantom {S_{\nu }^G} \big \rvert - \bigg \lvert \int _0^t {\underbrace {s_{\nu }^G \left (\tilde {t}\right )}_{\leq \, 0}}\vphantom {s_{\nu }^G \left (\tilde {t}\right )} \, \mathrm d \tilde {t} \bigg \rvert \quad \text {with} \quad S_{\nu }^G :\!= \int _0^T s_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t} \end{equation*}

for source vertices $\nu \in \partial ^+ \mathcal {V}$ , giving the remaining supply at time $t$ , which has not entered the network yet. Similarly, the sink vertex mass can be defined as

\begin{equation*} D_{\nu } \left (t\right ) :\!= \int _0^t d_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t} \quad \text {with} \quad D_{\nu }^G :\!= \int _0^T d_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t} \end{equation*}

for sink vertices $\nu \in \partial ^- \mathcal {V}$ , encoding the amount of gas which has already exited the network until time $t$ .

Proposition 2.10 (Source vertex mass density). The source vertex mass density $S_{\nu } \left (t\right )$ of a source vertex $\nu \in \partial ^+ \mathcal {V}$ at any time $t \in \left [0, T\right ]$ can be calculated by

\begin{equation*} S_{\nu }\left (t\right ) = \left \lvert \int _t^T s_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t} \right \rvert \end{equation*}

for time-dependent boundary conditions, or with $s_{\nu }$ instead of $s_{\nu }^G$ for time-independent boundary conditions.

With the vertex flux functions and the vertex mass densities at boundary vertices, we can also calculate the total gas supply $\left \lvert S \right \rvert$ (with $S \leq 0$ ) entering the network and the total gas demand $D$ exiting the network in the time period $\left [0, T\right ]$ .

Observation 2.11 (Total supply and total demand). In the case of time-dependent boundary conditions, the total supply and total demand of the network for the time period $\left [0, T\right ]$ are given by

\begin{align*} S & = \sum _{\nu \in \partial ^+ \mathcal {V}} - S_{\nu }\left (0\right ) = \sum _{\nu \in \partial ^+ \mathcal {V}} S_{\nu }^G = \sum _{\nu \in \partial ^+ \mathcal {V}} \int _0^T s_{\nu }^G \left (t\right ) \, \mathrm d t, \\ D & = \sum _{\nu \in \partial ^- \mathcal {V}} D_{\nu } \left (T\right ) = \sum _{\nu \in \partial ^- \mathcal {V}} D_{\nu }^G = \sum _{\nu \in \partial ^- \mathcal {V}} \int _0^T d_{\nu }^G \left (t\right ) \, \mathrm d t. \end{align*}

Alternatively, exactly the same equations hold true with $s_{\nu }$ instead of $s_{\nu }^G$ and $d_{\nu }$ instead of $d_{\nu }^G$ for time-independent boundary conditions.

Note that, dead-end vertices can in general not be considered a special case of sink or source vertices, as they are interior vertices and thus provide a storage opportunity, which is not given for boundary vertices. However, the two definitions coincide in the case of a no-flux boundary condition, so $f_{\nu } = s_{\nu }^G = d_{\nu }^G \equiv 0$ (or alternatively for $f_{\nu } \equiv 0$ and $S_{\nu }^G = D_{\nu }^G = 0$ ).

3.1. Measure spaces on graphs

By $\mathcal {M}_+\left (\Omega \right )$ we denote the set of Borel measures on a metric space $\left (\Omega, d\right )$ with $\Omega \subseteq \mathbb {R}^k$ and $k \in \mathbb {N}$ (and $\Sigma$ being the corresponding $\sigma$ -algebra of Borel sets of $\Omega$ ), which are non-negative and bounded, so $\forall X \subseteq \Omega$ we require $0 \leq \mu \left (X\right )$ , and also $\mu \left (\Omega \right ) \lt \infty$ . Based on this notation, we choose for each edge $e \in \mathcal {E}$ the mass density $\rho _e$ from the set of non-negative bounded measures on edge $e$ for the time period $\left [0, T\right ]$ with

\begin{equation*} \mathcal {M}_+(e) :\!= \mathcal {M}_+\left (\left [0, L_e\right ] \times \left [0, T\right ]\right ). \end{equation*}

Furthermore, we choose the vertex mass density $\gamma _{\nu }$ for each interior vertex $\nu \in \stackrel{\small\circ}{\mathcal {V}}$ , the source vertex mass $S_{\nu }$ for each $\nu \in \partial ^+ \mathcal {V}$ and the sink vertex mass $D_{\nu }$ for each $\nu \in \partial ^- \mathcal {V}$ , from the set of non-negative bounded measures on vertex $\nu$ for the time period $\left [0, T\right ]$ , so we define

\begin{equation*} \mathcal {M}_+\left (\nu \right ) :\!= \mathcal {M}_+\left (\left \{\nu \right \} \times \left [0, T\right ]\right ). \end{equation*}

For a fixed time $t \in \left [0, T\right ]$ , we define in a similar manner $\mathcal {M}_+^t(e) :\!= \mathcal {M}_+\left (\left [0, L_e\right ]\right )$ and $\mathcal {M}_+^t\left (\nu \right ) :\!= \mathcal {M}_+\left (\left \{\nu \right \}\right ) = a \delta _{\nu }$ , for an edge $e \in \mathcal {E}$ and a vertex $\nu \in \mathcal {V}$ , where $\delta _{\nu }$ denotes the corresponding Dirac measure with $a \in \mathbb {R}_{\geq 0}$ .

For the velocity functions $v_e$ , the mass flux functions $j_e$ of edges $e \in \mathcal {E}$ , and for the vertex excess flux function $f_{\nu }$ of interior vertices $\nu \in \stackrel{\small\circ}{\mathcal {V}}$ , we do not need the non-negativity of the defined measure spaces. Hence, by $\mathcal {M}\left (\Omega \right )$ we denote the set of bounded Borel measures $\mu$ on a metric space $\left (\Omega, d\right )$ , and thus we define

\begin{equation*} \mathcal {M}(e) :\!= \mathcal {M}\left (\left [0, L_e\right ] \times \left [0, T\right ]\right ) \quad \text {and} \quad \mathcal {M}\left (\nu \right ) :\!= \mathcal {M}\left (\left \{\nu \right \} \times \left [0, T\right ]\right ) \end{equation*}

for any edge $e \in \mathcal {E}$ and any vertex $\nu \in \mathcal {V}$ .

With the set of measures being defined on a single edge or a single vertex of the graph, we can now also introduce coupled measures on a subset of edges or a subset of vertices. These coupled measure are defined on the domain of a cartesian product of $\mathcal {M}_+(e)$ , $\mathcal {M}_+^t(e)$ , $\mathcal {M} (e)$ , $\mathcal {M}_+\left (\nu \right )$ , $\mathcal {M}_+^t\left (\nu \right )$ or $\mathcal {M} \left (\nu \right )$ for edges $e \in \mathcal {E}$ and vertices $\nu \in \mathcal {V}$ . A formal definition can be found in the appendix A.

3.2. Formulation of the optimal transport problems

Depending on the the existence of boundary vertices and the type of boundary conditions, as well as the coupling conditions at interior vertices, we obtain different formulations of the optimal transport problem on the graph $\mathcal {G} = \left (\mathcal {V}, \mathcal {E}\right )$ .

The aim of this subsection is to introduce a variation of dynamical optimal transport problems in the sense of Benamou and Brenier, [Reference Benamou and Brenier5], i.e. we consider optimisation problems on the graph over curves satisfying the continuity equation together with a version of Kirchhoff’s law as a coupling condition at the vertices, which ensures global mass conservation of our transport. These conditions will serve as side constraints to an optimization problem with a (for now) general cost functional $c$ , which can be thought of as dynamical costs that are chosen depending on the application. We do not yet specify the exact transport costs, as we first focus on the modelling of mass transfer between supply, demand, storage and the network. As we will see in section 5, building a minimizing movement scheme using these optimal transport problems as the distance, we will indeed recover the (ISO3) model.

Assuming the case with no boundary vertices and classical Kirchhoff’s law as coupling condition (as in [Reference Erbar, Forkert, Maas and Mugnolo14]), the corresponding optimal transport problem can be written as:

\begin{alignat*} {3} && &\inf _{\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, v \in \mathcal {M}\left (\mathcal {E}\right )} \quad c\left (\rho, v\right ) \qquad \text {subject to}\\&& 0 & = \frac {\partial \rho _e}{\partial t} + \frac {\partial }{\partial x}\left (\rho _e v_e\right ) & \forall e \in \mathcal {E}, \, x \in \left [0, L_e\right ], \, t \in \left [0, T\right ] \\ && 0 & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . \left (\rho _e v_e\right ) \right |_{x = L_e} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . \left (\rho _e v_e\right ) \right |_{x = 0} \qquad \qquad \quad & \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}, \, t \in \left [0, T\right ] \\ && \left .\rho _e\right |_{t = 0} & = \left (\rho _0\right )_e, \, \left .\rho _e\right |_{t = T} = \left (\rho _T\right )_e & \forall e \in \mathcal {E} \end{alignat*}

If we rewrite the presented optimal transport problem with the mass flux $j_e$ instead of the product $\rho _e v_e$ for each edge $e \in \mathcal {E}$ , we obtain a convex problem with linear side constraints. Based on this form, we obtain the following extensions for different combinations of coupling- and boundary conditions.

Definition 3.1 (Optimal transport problem). If there are no boundary vertices ( $\left \lvert \partial \mathcal {V}\right \rvert = 0$ so $\mathcal {V} = \stackrel{\small\circ}{\mathcal {V}}$ ) and thus no boundary conditions, we can formulate the optimal transport problem with coupling conditions based on the generalized Kirchhoff’s law (VF KL) as:

(OT VF KL) \begin{align} &\inf _{\substack {\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, j \in \mathcal {M}\left (\mathcal {E}\right ), \\ \gamma \in \mathcal {M}_+ \left (\stackrel{\small\circ}{\mathcal {V}}\right ), \, f \in \mathcal {M}\left (\stackrel{\small\circ}{\mathcal {V}}\right )}} \quad c\left (\rho, j, \gamma, f\right ) \qquad \text {subject to}\end{align}

(CE 1) \begin{align} 0 &= \frac {\partial \rho _e}{\partial t} + \frac {\partial j_e}{\partial x} & \forall e \in \mathcal {E}, \, x \in \left [0, L_e\right ], \, t \in \left [0, T\right ]\end{align}

(VF KL 1) \begin{align} f_{\nu }\left (t\right ) & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . j_e \right |_{x = L_e} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . j_e \right |_{x = 0} & \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}, \, t \in \left [0, T\right ]\end{align}

(VF KL 2) \begin{align}\frac {\partial \gamma _{\nu }}{\partial t} & =\ f_{\nu } & \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}, \, t \in \left [0, T\right ]\end{align}

(CE 2) \begin{align}\left .\rho _e\right |_{t = 0} & = \left (\rho _0\right )_e, \, \left .\rho _e\right |_{t = T} = \left (\rho _T\right )_e & \forall e \in \mathcal {E}\end{align}

(VF KL 3) \begin{align} \left .\gamma _{\nu }\right |_{t = 0} & = \left (\gamma _0\right )_{\nu }, \, \left .\gamma _{\nu }\right |_{t = T} = \left (\gamma _T\right )_{\nu } & \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}\end{align}

When using the classical version of Kirchhoff’s law (C KL) at the interior vertices (corresponding to no storage of gas at interior vertices), then this simplifies to:

(OT C KL) \begin{align} &\inf _{\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, j \in \mathcal {M}\left (\mathcal {E}\right )} \quad c\left (\rho, j\right ) \qquad \text {subject to} \end{align}

(CE 1) \begin{align} 0 & = \frac {\partial \rho _e}{\partial t} + \frac {\partial j_e}{\partial x} & \forall e \in \mathcal {E}, \, x \in \left [0, L_e\right ], \, t \in \left [0, T\right ] \end{align}

(C KL 1) \begin{align} 0 & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . j_e \right |_{x = L_e} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . j_e \right |_{x = 0} & \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}, \, t \in \left [0, T\right ] \end{align}

(CE 2) \begin{align} \left .\rho _e\right |_{t = 0} & = \left (\rho _0\right )_e, \, \left .\rho _e\right |_{t = T} = \left (\rho _T\right )_e & \forall e \in \mathcal {E} \end{align}

The optimal transport problem (OT VF KL) is exactly the same as the model studied in [Reference Burger, Humpert and Pietschmann8], while (OT C KL) was analysed in [Reference Erbar, Forkert, Maas and Mugnolo14]. In both cases, no boundary vertices are present and thus mass transport into or out of the network is not possible. In view of our application, we hence extend these problems by the following boundary conditions.

Remark 3.2 (Extension with boundary conditions). If boundary vertices are present ( $\left \lvert \partial \mathcal {V}\right \rvert \gt 0$ so $\mathcal {V} = \partial \mathcal {V} \cup \stackrel{\small\circ}{\mathcal {V}} \supset \stackrel{\small\circ}{\mathcal {V}}$ ), then depending on the type of boundary conditions, different constraints need to be added to the optimal transport problem. Assuming time-dependent boundary conditions, then the following constraints need to be included:

(TD BC 1) \begin{align} 0 & = \Bigg ( \left \lvert s_{\nu }^G \right \rvert + \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . j_e \right |_{x = L_e} \Bigg ) - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . j_e \right |_{x = 0} & \forall \nu \in \partial ^+ \mathcal {V}, \, t \in \left [0, T\right ] \end{align}

(TD BC 2) \begin{align} 0 & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . j_e \right |_{x = L_e} - \Bigg ( d_{\nu }^G + \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . j_e \right |_{x = 0} \Bigg ) & \forall \nu \in \partial ^- \mathcal {V}, \, t \in \left [0, T\right ]\end{align}

Similarly, for time-independent boundary conditions, the following constraints, as well as initial and final conditions need to be included in the optimal transport problem instead:

(TI BC 1) \begin{align} s_{\nu } & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . j_e \right |_{x = 0} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . j_e \right |_{x = L_e} & \forall \nu \in \partial ^+ \mathcal {V}, \, t \in \left [0, T\right ] \end{align}

(TI BC 2) \begin{align} d_{\nu } & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . j_e \right |_{x = L_e} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . j_e \right |_{x = 0} & \forall \nu \in \partial ^- \mathcal {V}, \, t \in \left [0, T\right ] \end{align}

(TI BC 3) \begin{align} \frac {\partial S_{\nu }}{\partial t} & = s_{\nu } & \forall \nu \in \partial ^+ \mathcal {V}, \, t \in \left [0, T\right ] \end{align}

(TI BC 4) \begin{align} \frac {\partial D_{\nu }}{\partial t} & = d_{\nu } & \forall \nu \in \partial ^- \mathcal {V}, \, t \in \left [0, T\right ] \end{align}

(TI BC 5) \begin{align} \left .S_{\nu }\right |_{t = 0} & = \left \lvert S_{\nu }^G \right \rvert, \, \left .S_{\nu }\right |_{t = T} = 0 & \forall \nu \in \partial ^+ \mathcal {V} \end{align}

(TI BC 6) \begin{align} \left .D_{\nu }\right |_{t = 0} & = 0, \, \left .D_{\nu }\right |_{t = T} = D_{\nu }^G & \forall \nu \in \partial ^- \mathcal {V} \end{align}

In order to avoid $s_{\nu }\left (t\right ) \gt 0$ or $d_{\nu }\left (t\right ) \lt 0$ for all $t \in \left [0, T\right ]$ , we require $-s \in \mathcal {M}_+ \left (\partial ^+ \mathcal {V}\right )$ and $d \in \mathcal {M}_+ \left (\partial ^- \mathcal {V}\right )$ . In addition, for time-independent boundary conditions, the cost functional changes to

\begin{equation*} \inf _{\substack {\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, j \in \mathcal {M}\left (\mathcal {E}\right ), \\ \gamma \in \mathcal {M}_+ \left (\stackrel{\small\circ}{\mathcal {V}}\right ), \, f \in \mathcal {M}\left (\stackrel{\small\circ}{\mathcal {V}}\right ), \\ S, \, -s \in \mathcal {M}_+ \left (\partial ^+ \mathcal {V}\right ), \\ D, \, d \in \mathcal {M}_+ \left (\partial ^- \mathcal {V}\right )}} c\left (\rho, j, \gamma, f, S, s, D, d\right ), \qquad \text {or} \qquad \inf _{\substack {\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, j \in \mathcal {M}\left (\mathcal {E}\right ), \\ S, \, -s \in \mathcal {M}_+ \left (\partial ^+ \mathcal {V}\right ), \\ D, \, d \in \mathcal {M}_+ \left (\partial ^- \mathcal {V}\right )}} c\left (\rho, j, S, s, D, d\right ), \end{equation*}

depending on the type of coupling conditions at the interior vertices.

Examples for possible (action or) cost functionals will be given in subsection 4.2, aiming at calculating the necessary kinetic energy for the transport, which is sensible for gas networks as well as being a typical type of action functional used for dynamic optimal transport.

3.3. Feasibility of the optimal transport problem

Depending on whether we allow storage of gas at interior vertices $\stackrel{\small\circ}{\mathcal {V}}$ or not, and depending on the existence of boundary vertices and the type of boundary conditions, necessary conditions for the feasibility of the optimal transport problem differ slightly.

So far, the continuity equation and the Kirchhoff’s law respectively only ensure local mass conservation; however, we are naturally interested in analysing global mass conservation for the entire network.

Observation 3.4 (Total mass). The contextual interpretation of the different terms in the total mass $\Sigma \left (t\right )$ is given in table 1 . Moreover, note that for time-dependent boundary conditions, the source vertex mass density $S_{\nu }$ and sink vertex mass density $D_{\nu }$ at time $t \in \left [0, T\right ]$ can easily be calculated by utilizing the given source vertex flux $s_{\nu }^G$ and sink vertex flux $d_{\nu }^G$ with

Table 1. Contributions to the total mass

\begin{equation*} S_{\nu } \left (t\right ) = \left \lvert \int _t^T s_{\nu }^G (\tilde {t}\,) \, \mathrm d \tilde {t} \right \rvert, \qquad \text {and} \qquad D_{\nu } (t) = \int _0^t d_{\nu }^G \left (\tilde {t\,}\right ) \, \mathrm d \tilde {t}. \end{equation*}

The condition of global mass conservation can thus be written as the following.

Definition 3.5 (Global mass conservation). The global mass conservation condition in a graph $\mathcal {G} = \left (\mathcal {V}, \mathcal {E}\right )$ is given by

\begin{equation*} \frac {\partial \Sigma \left (t\right )}{\partial t} = 0 \qquad \forall t \in \left (0, T\right ). \end{equation*}

Remark 3.6 (Global mass conservation). For the measures $\rho _e$ , $\gamma _{\nu }$ , $S_{\nu }$ and $D_{\nu }$ , we will always assume that they are chosen in such a way that the total mass in the network is constant for all $t \in \left [0, T\right ]$ , so

\begin{equation*} \Sigma \left (t\right ) \equiv C \in \mathbb {R}_{\geq 0}. \end{equation*}

Hence, we especially demand that the initial and final conditions fulfill

\begin{equation*} \Sigma \left (0\right ) = \Sigma \left (T\right ) = C, \end{equation*}

in order to obtain a feasible optimal transport problem.

The global mass conservation conditions can be used to define feasible upper bounds for the total demand of all sink vertices for the time period $\left [0, T\right ]$ .

Proposition 3.7 (Upper bound for the demand with time-independent boundary conditions). In the case of generalized Kirchhoff’s law at interior vertices and the existence of boundary vertices with time-independent boundary conditions, the accumulated demand $D = \sum _{\nu \in \partial ^- \mathcal {V}} D_{\nu } \left (T\right )$ of the time period $\left [0, T\right ]$ is bounded from above through the inequality

\begin{equation*} \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \left \lvert S \right \rvert + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x \geq D, \end{equation*}

where $\left \lvert S \right \rvert = \sum _{\nu \in \partial ^- \mathcal {V}} S_{\nu } \left (0\right )$ corresponds to the accumulated supply for the time period $\left [0, T\right ]$ .

Proof. If we compare the total mass $\Sigma \left (t\right )$ for any $t \in \left [0, T\right ]$ to the initial total mass $\Sigma \left (0\right )$ , then we obtain the equation

\begin{align*} & \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} \underbrace {S_{\nu } \left (0\right )}_{= \, \left \lvert S_{\nu }^G \right \rvert } + \sum _{\nu \in \partial ^- \mathcal {V}} \underbrace {D_{\nu } \left (0\right )}_{= \, 0} + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x = \\ & \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (t\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} S_{\nu } \left (t\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} D_{\nu } \left (t\right ) + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, t\right ) \, \mathrm d x, \end{align*}

which can be rewritten to

(GMC) \begin{align} \begin{split} & \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} \underbrace {\left (\left \lvert S_{\nu }^G \right \rvert - S_{\nu } \left (t\right )\right )}_{= \, \left \lvert \int _0^t s_{\nu } \left (\tilde {t}\right ) \, \mathrm d \tilde {t} \right \rvert } + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x = \\ & \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (t\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} \underbrace {D_{\nu } \left (t\right )}_{= \, \int _0^t d_{\nu } \left (\tilde {t}\right ) \, \mathrm d \tilde {t}} + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, t\right ) \, \mathrm d x. \end{split} \end{align}

The left side of the global mass conservation equation (GMC) calculates the amount of gas stored in the network at time $t = 0$ (in interior vertices or in the pipes) plus the gas, which has entered the network through source vertices in the time period $\left [0, t\right ]$ . The right side encodes the amount of gas in the interior vertices or the pipes at time $t$ together with the gas, which has exited the network through sink vertices in the time period $\left [0, t\right ]$ .

For $t = T$ we obtain from the GMC

\begin{align*} & \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} \left \lvert \int _0^T s_{\nu } \left (t\right ) \, \mathrm d t \right \rvert + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x = \\ & \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (T\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} \int _0^T d_{\nu } \left (t\right ) \, \mathrm d t + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, T\right ) \, \mathrm d x, \end{align*}

which is equivalent to

\begin{equation*} \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \left \lvert S \right \rvert + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x = \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (T\right ) + D + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, T\right ) \, \mathrm d x. \end{equation*}

By utilizing the non-negativity of the vertex mass density $\gamma$ and of the mass density $\rho$

\begin{equation*} \gamma _{\nu } \in \mathcal {M}_+ \left (\stackrel{\small\circ}{\mathcal {V}}\right ), \, \rho _e \in \mathcal {M}_+ \left (\mathcal {E}\right ) \quad \Longrightarrow \quad \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (T\right ), \, \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, T\right ) \, \mathrm d x \geq 0, \end{equation*}

we obtain the upper bound for the accumulated demand of all sink vertices $\nu \in \partial ^- \mathcal {V}$

\begin{equation*} \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \left \lvert S \right \rvert + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x \geq D. \end{equation*}

Corollary 3.8 (Upper bound for the demand with time-dependent boundary conditions). If the graph has boundary vertices with time-dependent boundary conditions and uses generalized Kirchhoff’s law as coupling condition at interior vertices, the accumulated demand at all times $t \in \left [0, T\right ]$ is bounded above by

\begin{equation*} \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} \left \lvert \int _0^t s_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t}\right \rvert + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x \geq \sum _{\nu \partial ^- \mathcal {V}} \int _0^t d_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t}, \end{equation*}

which can be used to test the given source vertex flux $s_{\nu }^G$ and sink vertex flux $d_{\nu }^G$ for plausibility.

Proof. The non-negativity of the vertex mass density $\gamma$ and the mass density $\rho$ on the right side of the GMC

\begin{equation*} \gamma _{\nu } \in \mathcal {M}_+ \left (\stackrel{\small\circ}{\mathcal {V}}\right ), \, \rho _e \in \mathcal {M}_+ \left (\mathcal {E}\right ) \quad \Longrightarrow \quad \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (\tilde {t}\right ), \, \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, \tilde {t}\right ) \, \mathrm d x \geq 0, \end{equation*}

holds true for all $\tilde {t} \in \left [0, T\right ]$ . Therefore, we obtain the inequality

\begin{equation*} \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} \left (\left \lvert S_{\nu }^G \right \rvert - S_{\nu } \left (t\right )\right ) + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x \geq \sum _{\nu \partial ^- \mathcal {V}} D_{\nu } \left (t\right ). \end{equation*}

Plugging in the definitions of the source vertex mass densities $S_{\nu }^G$ and $S_{\nu }$ , as well as the sink vertex mass densities $D_{\nu }$ , yields

\begin{equation*} \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \gamma _{\nu } \left (0\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} \left \lvert \int _0^t s_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t}\right \rvert + \sum _{e \in \mathcal {E}} \int _0^{L_e} \rho _e \left (x, 0\right ) \, \mathrm d x \geq \sum _{\nu \partial ^- \mathcal {V}} \int _0^t d_{\nu }^G \left (\tilde {t}\right ) \, \mathrm d \tilde {t} \end{equation*}

for all $t \in \left [0, T\right ]$

The global mass conservation equation (GMC) together with the upper bound for the accumulated demand can be used to test given initial, final and boundary conditions for possibly leading to an infeasible optimal transport problem. However, note that these tests do not detect instances which for example require unfeasibly fast velocities $v$ , violating the necessary upper bound for the flow rates, when going from (ISO1) to (ISO3).

3.4. Strong and weak solutions of the continuity equation constraints

Assuming feasibility of the introduced optimal transport problem, we now want to define strong and weak solutions of the continuity equation constraints (CE 1), (CE 2) combined with the constraints from the coupling conditions and possibly the boundary conditions.

For the introduction of strong solutions, we need a variety of conditions:

1. (S1) $\forall e \in \mathcal {E}: \left (x, t\right ) \mapsto \rho _e \left (x, t\right )$ continuous on $\left [0, L_e\right ] \times \left [0, T\right ]$ $\forall e \in \mathcal {E}, \, \forall x \in \left [0, L_e\right ]: t \mapsto \rho _e \left (x, t\right )$ continuously differentiable on $\left (0, T\right )$

2. (S2) $\forall e \in \mathcal {E}: \left (x, t\right ) \mapsto j_e \left (x, t\right )$ continuous on $\left [0, L_e\right ] \times \left [0, T\right ]$ $\forall e \in \mathcal {E}, \, \forall t \in \left [0, T\right ]: x \mapsto j_e \left (x, t\right )$ continuously differentiable on $\left (0, L_e\right )$

3. 1. (S3) i) $\rho \in \mathcal {M}_+ \left (\mathcal {E}\right )$ , $j \in \mathcal {M} \left (\mathcal {E}\right )$ fulfill the constraints (CE 1), (CE 2) and (C KL 1)

   2. ii) $\rho \in \mathcal {M}_+ \left (\mathcal {E}\right )$ , $j \in \mathcal {M} \left (\mathcal {E}\right )$ , $\gamma \in \mathcal {M}_+ \left (\stackrel{\small\circ}{\mathcal {V}}\right )$ and $f \in \mathcal {M} \left (\stackrel{\small\circ}{\mathcal {V}}\right )$ fulfill the side constraints (CE 1), (CE 2) and (VF KL 1) - (VF KL 3)

4. (S4) $\forall \nu \in \stackrel{\small\circ}{\mathcal {V}}: t \mapsto \gamma _{\nu } \left (t\right )$ continuous on $\left [0, T\right ]$ , and continuously differentiable on $\left (0, T\right )$ (which directly implies $\forall \nu \in \stackrel{\small\circ}{\mathcal {V}}: t \mapsto f_{\nu } \left (t\right )$ continuous on $\left (0, T\right )$ )

5. (S5) $\forall \nu \in \partial ^+ \mathcal {V}: t \mapsto S_{\nu } \left (t\right )$ continuous on $\left [0, T\right ]$ , and continuously differentiable on $\left (0, T\right )$ (which directly implies $\forall \nu \in \partial ^+ \mathcal {V}: t \mapsto s_{\nu } \left (t\right )$ continuous on $\left (0, T\right )$ )

6. (S6) $\forall \nu \in \partial ^- \mathcal {V}: t \mapsto D_{\nu } \left (t\right )$ continuous on $\left [0, T\right ]$ , and continuously differentiable on $\left (0, T\right )$ (which directly implies $\forall \nu \in \partial ^- \mathcal {V}: t \mapsto d_{\nu } \left (t\right )$ continuous on $\left (0, T\right )$ )

7. (S7)

   1. i) $S, -s \in \mathcal {M}_+\left (\partial ^+ \mathcal {V}\right )$ and $D, d \in \mathcal {M}_+\left (\partial ^- \mathcal {V}\right )$ fulfill the constraints (TD BC 1), (TD BC 2)

   2. ii) $S, -s \in \mathcal {M}_+\left (\partial ^+ \mathcal {V}\right )$ and $D, d \in \mathcal {M}_+\left (\partial ^- \mathcal {V}\right )$ fulfill the constraints (TI BC 1) - (TI BC 6)

For the definition of weak solutions, we need suitable test functions defined on the domain of the graph

\begin{equation*} \Omega _{\mathcal {G}} :\!= \left (\times _{\nu \in \mathcal {V}} \left \{\nu \right \}\right ) \times \left (\times _{e \in \mathcal {E}} \left [0, L_e\right ]\right ). \end{equation*}

The following conditions for the definitions of weak solutions can be derived from the presented strong solutions. Here, $\mathcal {L}$ denotes the standard $1$ -dimensional Lebesgue measure on $\left [0, T\right ]$ .

1. (W1) $\rho \in \mathcal {M}_+ \left (\mathcal {E}\right )$ and $j \in \mathcal {M} \left (\mathcal {E}\right )$

2. (W2) $\rho$ fulfills the constraints (CE 2)

3. (W3)

   1. i) $j$ fulfills the constraint

      \begin{equation*} \int _0^T \sum _{e \in \mathcal {E}} \int _0^{L_e} \left \lvert\,\, j_e \left (x, t\right ) \right \rvert \, \mathrm {d} x \, \mathrm d t \lt \infty \end{equation*}

   2. ii) $j$ , $s$ and $d$ fulfill the constraint

      \begin{equation*} \int _0^T \Bigg (\sum _{e \in \mathcal {E}} \int _0^{L_e} \left \lvert\,\, j_e \left (x, t\right ) \right \rvert \, \mathrm d x - \sum _{\nu \in \partial ^+ \mathcal {V}} s_{\nu } \left (t\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} d_{\nu } \left (t\right )\Bigg ) \, \mathrm d t \lt \infty \end{equation*}

   3. iii) $j$ and $f$ fulfill the constraint

      \begin{equation*} \int _0^T \Bigg ( \sum _{e \in \mathcal {E}} \int _0^{L_e} \left \lvert\,\, j_e \left (x, t\right ) \right \rvert \, \mathrm d x + \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \left \lvert f_{\nu } \left (t\right ) \right \rvert \Bigg )\, \mathrm d t \lt \infty \end{equation*}

   4. iv) $j$ , $f$ , $s$ and $d$ fulfill the constraint

      \begin{equation*} \int _0^T \Bigg (\sum _{e \in \mathcal {E}} \int _0^{L_e} \left \lvert\, j_e \left (x, t\right ) \right \rvert \, \mathrm d x + \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} \left \lvert f_{\nu } \left (t\right ) \right \rvert - \sum _{\nu \in \partial ^+ \mathcal {V}} s_{\nu } \left (t\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} d_{\nu } \left (t\right )\Bigg ) \, \mathrm d t \lt \infty \end{equation*}

4. (W4) $\gamma \in \mathcal {M}_+ \left (\stackrel{\small\circ}{\mathcal {V}}\right )$ and $f \in \mathcal {M} \left (\stackrel{\small\circ}{\mathcal {V}}\right )$

5. (W5) $\gamma$ fulfills the constraint (VF KL 3), $\forall \nu \in \stackrel{\small\circ}{\mathcal {V}}$ and for $\mathcal {L}$ -a.e. $t \in \left (0, T\right ): \, \frac {\partial \gamma _{\nu } \left (t\right )}{\partial t} = f_{\nu } \left (t\right )$

6. (W6) $S, -s \in \mathcal {M}_+ \left (\partial ^+ \mathcal {V}\right )$

7. (W7) $D, d \in \mathcal {M}_+ \left (\partial ^- \mathcal {V}\right )$

8. (W8) $S$ fulfills the constraint (TI BC 5), $\forall \nu \in \partial ^+ \mathcal {V}$ and for $\mathcal {L}$ -a.e. $t \in \left (0, T\right ): \, \frac {\partial S_{\nu } \left (t\right )}{\partial t} = s_{\nu } \left (t\right )$ , $D$ fulfills the constraint (TI BC 6), $\forall \nu \in \partial ^- \mathcal {V}$ and for $\mathcal {L}$ -a.e. $t \in \left (0, T\right ): \, \frac {\partial D_{\nu } \left (t\right )}{\partial t} = d_{\nu } \left (t\right )$

9. (W9)

   1. i) $\forall$ test functions $\varphi$ and $\mathcal {L}$ -a.e. $t \in \left (0, T\right )$ :

      \begin{align*} & \sum _{e \in \mathcal {E}} \int _0^{L_e} \frac {\partial \rho _e \left (x, t\right )}{\partial t} \varphi _e \left (x, t\right ) \, \mathrm d x = \sum _{e \in \mathcal {E}} \int _0^{Le} j_e \left (x, t\right ) \frac {\partial \varphi _e \left (x, t\right )}{\partial x} \, \mathrm d x \end{align*}

   2. ii) $\forall$ test functions $\varphi$ and $\mathcal {L}$ -a.e. $t \in \left (0, T\right )$ :

      \begin{align*} & \sum _{e \in \mathcal {E}} \int _0^{L_e} \frac {\partial \rho _e \left (x, t\right )}{\partial t} \varphi _e \left (x, t\right ) \, \mathrm d x + \sum _{\nu \in \partial ^+ \mathcal {V}} s_{\nu } \left (t\right ) \varphi _{\nu } \left (t\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} d_{\nu } \left (t\right ) \varphi _{\nu } \left (t\right ) = \\ & \sum _{e \in \mathcal {E}} \int _0^{Le} j_e \left (x, t\right ) \frac {\partial \varphi _e \left (x, t\right )}{\partial x} \, \mathrm d x \end{align*}

   3. iii) $\forall$ test functions $\varphi$ and $\mathcal {L}$ -a.e. $t \in \left (0, T\right )$ :

      \begin{equation*} \sum _{e \in \mathcal {E}} \int _0^{L_e} \frac {\partial \rho _e \left (x, t\right )}{\partial t} \varphi _e \left (x, t\right ) \, \mathrm d x + \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} f_{\nu } \left (t\right ) \varphi _{\nu } \left (t\right ) = \sum _{e \in \mathcal {E}} \int _0^{Le} j_e \left (x, t\right ) \frac {\partial \varphi _e \left (x, t\right )}{\partial x} \, \mathrm d x \end{equation*}

   4. iv) $\forall$ test functions $\varphi$ and $\mathcal {L}$ -a.e. $t \in \left (0, T\right )$ :

      \begin{align*} & \sum _{e \in \mathcal {E}} \int _0^{L_e} \frac {\partial \rho _e \left (x, t\right )}{\partial t} \varphi _e \left (x, t\right ) \, \mathrm d x + \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} f_{\nu } \left (t\right ) \varphi _{\nu } \left (t\right ) + \sum _{\nu \in \partial ^+ \mathcal {V}} s_{\nu } \left (t\right ) \varphi _{\nu } \left (t\right ) \\ &+ \sum _{\nu \in \partial ^- \mathcal {V}} d_{\nu } \left (t\right ) \varphi _{\nu } \left (t\right ) = \sum _{e \in \mathcal {E}} \int _0^{Le} j_e \left (x, t\right ) \frac {\partial \varphi _e \left (x, t\right )}{\partial x} \, \mathrm d x \end{align*}

4.1. General $p$ -Wasserstein distance

Given a metric space $\left (\Omega, d\right )$ , where we assume the domain $\Omega \subseteq \mathbb {R}^k$ with $k \in \mathbb {N}$ to be bounded, the $p$ -Wasserstein distance can be defined in the form of a static formulation based on the distance function $d$ , or it can also be given in a dynamic formulation using absolutely continuous curves $\mu : \left [0, T\right ] \longrightarrow \mathcal {M}_+ \left (\Omega \right )$ .

In literature, the study of the $p$ -Wasserstein distance mostly utilizes the distance function $d \left (x, y\right ) :\!= \left \lvert x - y\right \rvert$ , where $\left \lvert \, \cdot \, \right \rvert$ denotes the vector norm in $\mathbb {R}^k$ . The following static formulation of the $p$ -Wasserstein distance could analogously be formulated for a general distance function $d$ , however for the corresponding dynamic formulation of the $p$ -Wasserstein distance, the cost needs to be associated with an action, see for example [Reference Cédric29, Chapter 7] for details.

Definition 4.1 ( $p$ -Wasserstein distance [Reference Ambrosio, Gigli and Savaré2]). For two probability measures $\mu _0, \mu _T \in \mathbb {P}\left (\Omega \right )$ , the static formulation of the $p$ -Wasserstein distance for $p \in \left [1, \infty \right )$ is defined as

\begin{equation*} \widetilde {W_p} \left (\mu _0, \mu _T\right ) = \min _{\pi \in \Pi \left (\mu _0, \mu _T\right )} \left \{\int _{\Omega \times \Omega } \left \lvert x - y\right \rvert ^p \, \mathrm d \pi \left (x, y\right )\right \}^{\frac {1}{p}}, \end{equation*}

with $\Pi \left (\mu _0, \mu _T\right )$ being the set of all joint probability distributions $\pi$ on $\Omega \times \Omega$ , with the respective marginals $\mu _0$ and $\mu _T$ , i.e. $\mu _0\left (x\right ) = \int _{\Omega } \pi \left (x, y\right ) \, \mathrm d y$ , and $\mu _T\left (y\right ) = \int _{\Omega } \pi \left (x, y\right ) \, \mathrm d x$ .

Remark 4.2 ( $p$ -Wasserstein distance). As we assume $\Omega$ to be bounded, the finiteness of the $p$ -th moments of $\mu _0$ and $\mu _T$ are not required. However, for non-negative bounded measures $\mu _0, \mu _t \in \mathcal {M}_+\left (\Omega \right )$ instead of probability measures, the $p$ -Wasserstein distance is only able to assign costs to balanced data, i.e. $\mu _0$ and $\mu _T$ need to have the same total mass

\begin{equation*} \int _{\Omega } \mu _0\left (z\right ) \, \mathrm d z = \int _{\Omega } \mu _T\left (z\right ) \, \mathrm d z. \end{equation*}

In their seminal work [Reference Benamou and Brenier5], Benamou and Brenier introduced a dynamical formulation of the $2$ -Wasserstein distance that was also extended to the case $p\neq 2$ . In this formulation the task is to minimise an action functional, which corresponds to the kinetic energy of curves connecting the initial measure $\mu _0$ to the final measure $\mu _T$ . This minimization is subject to the initial and final conditions, $\left .\mu \right |_{t = 0} = \mu _0$ and $\left .\mu \right |_{t = T} = \mu _T$ , as well as the continuity equation $\frac {\partial \mu }{\partial t} + \frac {\partial j}{\partial x} = 0$ and therefore closely resembles the models introduced in Section 2.

Proposition 4.3 ( $p$ -Wasserstein distance). For $\Omega = \mathbb {R}^k$ and two probability measures $\mu_{0}, \mu_{T} \in \mathbb {P}\left (\Omega \right )$ , the dynamic formulation of the $p$ -Wasserstein distance for $p \in \left [1, \infty \right )$ is given by

\begin{equation*} \widetilde {W_p} \left (\mu_{0}, \mu_{T}\right ) = \inf _{\mu _t, v_t} \left \{\int _0^T \int _{\Omega } \left \lvert v_t\left (x\right )\right \rvert ^p \, \mathrm d \mu _t\left (x\right ) \, \mathrm d t\right \}^{\frac {1}{p}}, \end{equation*}

where the infimum runs over curves of measures $\mu _t \in \mathbb {P}\left (\Omega \right )$ and $v_t \in \mathcal {M}(\Omega ^k)$ for $t \in [0,T]$ , which satisfy the continuity equation in the sense of distributions and are such that $\left. \mu \right|_{t = 0} = \mu_0, \left. \mu \right|_{t = T} = \mu_T$ .

Note that, to increase readability, we will always denote the densities of measures with respect to a given reference measure such as the Lebesgue measure by the same symbol as the original measures themselves.

The dynamic formulation of the $p$ -Wasserstein distance can also be written in terms of mass density and mass flux (instead of mass density and velocity), resulting in a convex minimization problem with linear side constraints. For this formulation we will write $\mu$ instead of $\rho$ , since later on, the measure $\mu$ will not only consist of the mass density $\rho$ , but also the vertex mass $\gamma$ , and optionally the source vertex mass $S$ and the sink vertex mass $D$ .

Furthermore, we will restrict the set of feasible mass densities to $\mu \in \mathcal {M}_+ \left (\Omega \times \left [0, T\right ]\right )$ and in order to avoid a positivity constraint for $\mu$ , when reformulating $j :\!= v \mu$ to $v = \frac {j}{\mu }$ , we define

(4) \begin{equation} h: \mathbb {R} \times \mathbb {R}_+ \longrightarrow \mathbb {R}_{\geq 0} \quad \left (a, b\right ) \mapsto h\left (a, b\right ) = \left \{ \begin{array}{lr} \frac {\left \lvert a \right \rvert ^p}{b^{p - 1}} \quad & b \gt 0, \\ 0 \quad & b = a = 0,\\ +\infty \quad & b = 0, \, a \neq 0. \end{array} \right . \end{equation}

Applying all these reformulations results in the following dynamic formulation of the $p$ -Wasserstein distance.

Observation 4.5 (Dynamic formulation of the $p$ -Wasserstein distance). Using a reparametrization in time argument (see, e.g. [Reference Ambrosio, Gigli and Savaré2 , Lemma 1.1.4]), the $p$ -Wasserstein distance for two probability measures $\mu _0, \mu _T \in \mathbb {P} \left (\Omega \right )$ on a convex and compact domain $\Omega$ is given as

(DYN-T) \begin{equation} \widetilde {W}_p^p \left (\mu _0, \mu _T\right ) = \end{equation}

\begin{equation*} \min _{\substack {\mu \in \mathcal {M}_+ \left (\Omega \times \left [0, T\right ]\right ), \\ j \in \mathcal {M} \left (\Omega \times \left [0, T\right ]\right )}} \left \{T^{p - 1} \int _0^T \int _{\Omega } h\left (j\left (x, t\right ), \mu \left (x, t\right )\right ) \, \mathrm d x \, \mathrm d t \, \middle | \, \begin{array}{l} \frac {\partial \mu }{\partial t} + \frac {\partial j}{\partial x} = 0, \\ \left . \mu \right |_{t = 0} = \mu _0, \, \left . \mu \right |_{t = T} = \mu _T \end{array} \right \}, \end{equation*}

Note that, this formulation of the $p$ -Wasserstein distance is restricted to physically feasible solutions, as $\widetilde {W}_p^p \left (\mu _0, \mu _T\right ) \lt \infty$ is only possible, if for $\mathcal {L}$ -a.e. $x \in \Omega$ and $t \in \left [0, T\right ]$ either

\begin{equation*} \mu \left (x, t\right ) \gt 0 \quad \text {or} \quad \mu \left (x, t\right ) = j \left (x, t\right ) = 0, \end{equation*}

so a non-zero flux is only possible if the mass density is positive. (Here, $\mathcal {L}$ denotes the $k + 1$ -dimensional Lebesgue measure on $\Omega \times \left [0, T\right ] \subseteq \mathbb {R}^{k + 1}$ .)

4.2. $p$ -Wasserstein distance on networks

In order to utilize the $p$ -Wasserstein distance as the cost functional for the transport along an edge $e \in \mathcal {E}$ , we can set $\Omega :\!= \left [0, L_e\right ]$ as well as $k :\!= 1$ , and with the notation of the previously introduced optimal transport problems we can thus write for a given initial mass distribution of gas $\left (\rho _0\right )_e \in \mathbb {P} \left (\left [0, L_e\right ]\right )$ , a given final distribution $\left (\rho _T\right )_e \in \mathbb {P} \left (\left [0, L_e\right ]\right )$ and the mass flux $j_e$ along the edge $e$

\begin{equation*} W_p^p \left (\left (\rho _0\right )_e, \left (\rho _T\right )_e\right ) = \min _{\substack {\rho _e \in \mathcal {M}_+ (e), \\ j_e \in \mathcal {M} (e)}} \left \{T^{p - 1} \int _0^T \int _0^{L_e} h\left (j_e\left (x, t\right ), \rho _e\left (x, t\right )\right ) \, \mathrm d x \, \mathrm d t \, \middle | \, \begin{array}{l} \frac {\partial \rho _e}{\partial t} + \frac {\partial j_e}{\partial x} = 0, \\ \left . \rho _e \right |_{t = 0} = \left (\rho _0\right )_e, \, \left . \rho _e \right |_{t = T} = \left (\rho _T\right )_e \end{array} \right \}. \\ \end{equation*}

For a single interior vertex $\nu \in \stackrel{\small\circ}{\mathcal {V}}$ , in the case of generalized Kirchhoff’s law as coupling condition, we could also write the cost of transport from an initial vertex mass density $\left (\gamma _0\right )_{\nu } \in \mathbb {P} \left (\left \{\nu \right \}\right )$ to a final vertex mass density $\left (\gamma _T\right )_{\nu } \in \mathbb {P} \left (\left \{\nu \right \}\right )$ with vertex excess flux $f_{\nu }$ as

\begin{equation*} W_p^p \left (\left (\gamma _0\right )_{\nu }, \left (\gamma _T\right )_{\nu }\right ) =\min _{\substack {\gamma _{\nu } \in \mathcal {M}_+ \left (\nu \right ), \\ f_{\nu } \in \mathcal {M} \left (\nu \right )}} \left \{T^{p - 1} \int _0^T h\left (f_{\nu }\left (t\right ), \gamma _{\nu }\left (t\right )\right ) \, \mathrm d t \, \middle | \, \begin{array}{l} \frac {\partial \gamma _{\nu }}{\partial t} = f_{\nu }, \\ \left . \gamma _{\nu } \right |_{t = 0} = \left (\gamma _0\right )_{\nu }, \, \left . \gamma _{\nu } \right |_{t = T} = \left (\gamma _T\right )_{\nu } \end{array} \right \}, \end{equation*}

where instead of the standard continuity equation, we have $\frac {\partial \gamma _{\nu }}{\partial t} = f_{\nu }$ to ensure mass conservation in the vertex.

Observation 4.6 ( $p$ -Wasserstein distance at interior vertex). The application of the dynamic formulation of the $p$ -Wasserstein distance to an interior vertex $\nu \in \stackrel{\small\circ}{\mathcal {V}}$ is not particularly interesting as $\left (\gamma _0\right )_{\nu }, \left (\gamma _T\right )_{\nu } \in \mathbb {P} \left (\left \{\nu \right \}\right )$ directly implies

\begin{equation*} \left (\gamma _0\right )_{\nu } = \left (\gamma _T\right )_{\nu } = \delta _{\nu } \end{equation*}

and as $h\left (f_{\nu }\left (t\right ), \gamma _{\nu }\left (t\right )\right ) \geq 0$ for all $\gamma _{\nu } \in \mathcal {M}_+ \left (\left \{\nu \right \} \times \left [0, T\right ]\right )$ and $f_{\nu } \in \mathcal {M} \left (\left \{\nu \right \} \times \left [0, T\right ]\right )$ , the unique strong solution of the dynamic formulation of the $p$ -Wasserstein distance is given by

\begin{equation*} \left (\gamma _0\right )_{\nu } = \left (\gamma _T\right )_{\nu } = \gamma _{\nu } \left (t\right ) \equiv \delta _{\nu } \qquad \text {and} \qquad f_{\nu } \left (t\right ) \equiv 0, \end{equation*}

because $f_{\nu } \left (t\right ) \neq 0$ for some $t \in \left [0, T\right ]$ would lead to

\begin{equation*} h\left (f_{\nu }\left (t\right ), \gamma _{\nu }\left (t\right )\right ) \in \mathbb {R}_{\gt 0} \, \text {for} \, \gamma _{\nu }\left (t\right ) \gt 0 \qquad \text {and} \qquad h\left (f_{\nu }\left (t\right ), \gamma _{\nu }\left (t\right )\right ) = \infty \, \text {for} \, \gamma _{\nu }\left (t\right ) = 0. \end{equation*}

The continuity of $t \mapsto \gamma _{\nu } \left (t\right )$ and $t \mapsto f_{\nu }$ thus would lead to $T^{p - 1} \int _0^T h\left (f_{\nu }\left (t\right ), \gamma _{\nu }\left (t\right )\right ) \, \mathrm d t \gt 0$ if there exists any $t \in \left [0, T\right ]$ with $f_{\nu } \left (t\right ) \neq 0$ .

By extending the previous remarks about the dynamic formulation of the $p$ -Wasserstein distance on a single edge or an interior vertex, we can define the $p$ -Wasserstein metric on general oriented graphs modelling network structures.

Definition 4.7 ( $p$ -Wasserstein distance on networks without boundary vertices). For an oriented graph modelling a network, the dynamic formulation of the $p$ -Wasserstein distance, depending on the coupling condition at interior vertices, is given by

• • Classical Kirchhoff’s law and no boundary vertices, with $\mu _0 :\!= \rho _0$ and $\mu _T :\!= \rho _T$ :

  \begin{align*} & W_p^p \left (\mu _0, \mu _T\right ) = \min _{\substack {\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \\ j \in \mathcal {M} \left (\mathcal {E}\right )}} \left \{T^{p - 1} \int _0^T \sum _{e \in \mathcal {E}} \left (\int _0^{L_e} h\left (j_e\left (x, t\right ), \rho _e\left (x, t\right )\right ) \, \mathrm d x\right ) \, \mathrm d t \, \left .\vphantom {\begin{array}{l} \forall e \in \mathcal {E} : \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \mathcal {V} : \, 3 \, \left .ii\right ), \, \left .4\right ), \, \left .5\right ) \end{array}}\right | {\begin{array}{l} \forall e \in \mathcal {E}: \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \mathcal {V}: \, \left .3\right ) \, \left .i\right ) \end{array}} \right \} \end{align*}
  In [Reference Erbar, Forkert, Maas and Mugnolo14], this case for $p = 2$ is considered, including a proof that the $2$ -Wasserstein distance with classical Kirchhoff’s law and no boundary vertices does indeed define a metric on a given metric graph.

• • Generalised Kirchhoff’s law and no boundary vertices, with $\mu _0 :\!= \left (\rho _0, \gamma _0\right )$ and $\mu _T :\!= \left (\rho _T, \gamma _T\right )$ :

  \begin{align*} W_p^p \left (\mu _0, \mu _T\right ) = \min _{\substack {\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, j \in \mathcal {M} \left (\mathcal {E}\right ), \\ \gamma \in \mathcal {M}_+ \left (\mathcal {V}\right ), \, f \in \mathcal {M} \left (\mathcal {V}\right )}} &\left \{ T^{p - 1} \int _0^T \sum _{e \in \mathcal {E}} \left (\int _0^{L_e} h\left (j_e\left (x, t\right ), \rho _e\left (x, t\right )\right ) \, \mathrm d x\right ) \right .\\ & \left .+ \sum _{\nu \in \mathcal {V}} h\left (f_{\nu }\left (t\right ), \gamma _{\nu }\left (t\right )\right ) \, \mathrm d t \left .\vphantom {\begin{array}{l} \forall e \in \mathcal {E} : \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \mathcal {V} : \, 3 \, \left .ii\right ), \, \left .4\right ), \, \left .5\right ) \end{array}}\right | \, {\begin{array}{l} \forall e \in \mathcal {E} : \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \mathcal {V} : \, 3 \, \left .ii\right ), \, \left .4\right ), \, \left .5\right ) \end{array}} \right \} \end{align*}
  This case with $p = 2$ is analysed in [Reference Burger, Humpert and Pietschmann8], including a proof that the $2$ -Wasserstein distance with generalized Kirchhoff’s law and no boundary vertices does define a metric on a given metric graph.

The respective side constraints are given by

\begin{align*} \begin{array}{ll} \left .1\right )\,\, \frac {\partial \rho _e}{\partial t} + \frac {\partial j_e}{\partial x} = 0 &\qquad \qquad \left .4\right ) \frac {\partial \gamma _{\nu }}{\partial t} = f_{\nu }\\ \left .2\right )\,\, \left . \rho _e \right |_{t = 0} = \left (\rho _0\right )_e, \, \left . \rho _e \right |_{t = T} = \left (\rho _T\right )_e &\qquad \qquad \left .5\right )\,\, \left . \gamma _{\nu } \right |_{t = 0} = \left (\gamma _0\right )_{\nu }, \, \left . \gamma _{\nu } \right |_{t = T} = \left (\gamma _T\right )_{\nu }\\ 3 \left .i\right )\,\, 0 = \sum _{e \in \mathcal {E}: \, \delta ^E(e) = \nu } {\left . j_e \right |}_{x = L_e} - \sum _{e \in \mathcal {E}: \, \delta ^S(e) = \nu } {\left . j_e \right |}_{x = 0}&\\ 3 \left .ii\right )\,\, f_{\nu } = \sum _{e \in \mathcal {E}: \, \delta ^E(e) = \nu } {\left . j_e \right |}_{x = L_e} - \sum _{e \in \mathcal {E}: \, \delta ^S(e) = \nu } {\left . j_e \right |}_{x = 0}&\\ \end{array} \end{align*}

This definition of the $p$ -Wasserstein distance allows for a specific type of unbalanced optimal transport, as mass can be moved from edges (the network) to vertices (storage opportunity), and vice versa, for additional costs, i.e. kinetic energy. In literature, this is known as the Wasserstein–Fisher–Rao metric or the Hellinger-Kantorovich distance (cf. [Reference Monsaingeon24]). The concept of unbalanced optimal transport can be further generalized to also include boundary vertices and therefore the possibility to move gas from interior vertices and edges to boundary vertices (supply and demand).

Remark 4.8 ( $p$ -Wasserstein distance on networks with boundary vertices). For an oriented graph with boundary vertices, the dynamic formulation of the $p$ -Wasserstein distance, depending on the coupling condition at interior vertices, is given by

• • Classical Kirchhoff’s law and time-dependent or time-independent boundary conditions, with $\mu _0 :\!= \left (\rho _0, S_0. D_0\right )$ and $\mu _T :\!= \left (\rho _T, S_T, D_T\right )$ :

  \begin{align*} & W_p^p \left (\mu _0, \mu _T\right ) = \min _{\substack {\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, j \in \mathcal {M} \left (\mathcal {E}\right ), \\ S, -s \in \mathcal {M}_+ \left (\partial ^+ \mathcal {V}\right ), \\ D, d \in \mathcal {M}_+ \left (\partial ^- \mathcal {V}\right )}} \left \{T^{p - 1} \int _0^T \sum _{e \in \mathcal {E}} \left (\int _0^{L_e} h\left (j_e\left (x, t\right ), \rho _e\left (x, t\right )\right ) \, \mathrm d x\right ) \right .\\ &\left . + \sum _{\nu \in \partial ^+ \mathcal {V}} h_s\left (s_{\nu }\left (t\right ), S_{\nu }\left (t\right )\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} h_d\left (d_{\nu }\left (t\right ), D_{\nu }\left (t\right )\right ) \, \mathrm d t \, \left .\vphantom {\begin{array}{l} \forall e \in \mathcal {E} : \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}: \, 3 \, \left .i\right ) \\ \forall \nu \in \partial ^+ \mathcal {V}: \, \left .6\right ), \, \left .7\right ), \, \left .8\right ) \\ \forall \nu \in \partial ^- \mathcal {V} : \, \left .9\right ), \, \left .10\right ), \, \left .11\right ) \end{array}}\right | \, \begin{array}{l} \forall e \in \mathcal {E} : \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}: \, 3 \, \left .i\right ) \\ \forall \nu \in \partial ^+ \mathcal {V}: \, \left .6\right ), \, \left .7\right ), \, \left .8\right ) \\ \forall \nu \in \partial ^- \mathcal {V} : \, \left .9\right ), \, \left .10\right ), \, \left .11\right ) \end{array} \right \} . \end{align*}

• • Generalized Kirchhoff’s law and time-dependent or time-independent boundary conditions, with $\mu _0 :\!= \left (\rho _0, \gamma _0, S_0, D_0\right )$ and $\mu _T :\!= \left (\rho _T, \gamma _T, S_T, D_T\right )$ :

  \begin{align*} & W_p^p \left (\mu _0, \mu _T\right ) = \\ & \min _{\substack {\rho \in \mathcal {M}_+ \left (\mathcal {E}\right ), \, j \in \mathcal {M} \left (\mathcal {E}\right ), \\ \gamma \in \mathcal {M}_+ \left (\stackrel{\small\circ}{\mathcal {V}}\right ), \, f \in \mathcal {M} \left (\stackrel{\small\circ}{\mathcal {V}}\right ), \\ S, -s \in \mathcal {M}_+ \left (\partial ^+ \mathcal {V}\right ), \\ D, d \in \mathcal {M}_+ \left (\partial ^- \mathcal {V}\right )}} \left \{T^{p - 1} \int _0^T \sum _{e \in \mathcal {E}} \left (\int _0^{L_e} h\left (j_e\left (x, t\right ), \rho _e\left (x, t\right )\right ) \, \mathrm d x\right ) + \sum _{\nu \in \stackrel{\small\circ}{\mathcal {V}}} h\left (f_{\nu }\left (t\right ), \gamma _{\nu }\left (t\right )\right ) \right .\\ &\left . + \sum _{\nu \in \partial ^+ \mathcal {V}} h_s\left (s_{\nu }\left (t\right ), S_{\nu }\left (t\right ),t\right ) + \sum _{\nu \in \partial ^- \mathcal {V}} h_d\left (d_{\nu }\left (t\right ), D_{\nu }\left (t\right ),t\right ) \, \mathrm d t \, \left .\vphantom {\begin{array}{l} \forall e \in \mathcal {E} : \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}: \, 3 \, \left .i\right ) \\ \forall \nu \in \partial ^+ \mathcal {V}: \, \left .6\right ), \, \left .7\right ), \, \left .8\right ) \\ \forall \nu \in \partial ^- \mathcal {V} : \, \left .9\right ), \, \left .10\right ), \, \left .11\right ) \end{array}}\right | \, \begin{array}{l} \forall e \in \mathcal {E} : \, \left .1\right ), \, \left .2\right ) \\ \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}: \, 3 \, \left .ii\right ), \, \left .4\right ), \, \left .5\right ) \\ \forall \nu \in \partial ^+ \mathcal {V}: \, \left .6\right ), \, \left .7\right ), \, \left .8\right ) \\ \forall \nu \in \partial ^- \mathcal {V} : \, \left .9\right ), \, \left .10\right ), \, \left .11\right ) \end{array} \right \}, \end{align*}

The function $h$ is defined in ( 4 ). As for the supply and demand contributions, we allow for different costs, $h_s$ and $h_d$ , that may also depend on time and have to be chosen depending on the application. For example, they could include effects of changing gas prices. Moreover, the additional constraints are given by

\begin{align*} \begin{array}{ll} \left .6\right )\,\, s_{\nu } = \sum _{e \in \mathcal {E}: \, \delta ^E(e) = \nu } \left . j_e \right |_{x = L_e} - \sum _{e \in \mathcal {E}: \, \delta ^S(e) = \nu } \left . j_e \right |_{x = 0} & \qquad \qquad\left .9\right )\,\, d_{\nu } = \sum _{e \in \mathcal {E}: \, \delta ^E(e) = \nu } \left . j_e \right |_{x = L_e} - \sum _{e \in \mathcal {E}: \, \delta ^S(e) = \nu } \left . j_e \right |_{x = 0}\\ \left .7\right )\,\, \frac {\partial S_{\nu }}{\partial t} = s_{\nu }&\qquad \qquad \left .10\right )\,\, \frac {\partial D_{\nu }}{\partial t} = d_{\nu }\\ \left .8\right )\,\, \left . S_{\nu } \right |_{t = 0} = \left \lvert S_{\nu }^G \right \rvert, \, \left . S_{\nu } \right |_{t = T} = 0 &\qquad \qquad \left .11\right )\,\, \left . D_{\nu } \right |_{t = 0} = 0, \, \left . D_{\nu } \right |_{t = T} = D_{\nu }^G\\ \end{array} \end{align*}

5.1. $p$ -Wasserstein gradient flow on networks

In the following, let us derive a more concrete structure for Wasserstein gradient flows on networks (for the case of classical Kirchhoff’s law at interior vertices and no boundary vertices). For this sake we consider energies of the natural form

\begin{equation*} E(\mu ) = \sum _{e\in \mathcal {E}} E_e \left (\rho _e\right ), \end{equation*}

i.e. an additive composition of functionals on each edge $e \in \mathcal {E}$ .

With this choice, our next goal is to derive optimality conditions for (5). A formal derivation can be obtained from assuming the existence of a saddle-point of the associated Lagrange functional, whose first variations can be set to zero. We first note that

\begin{equation*} \frac {1}{p \tau ^{p-1}}W_p^p \left (\left (\rho _t\right )_e, \left (\rho _{t+\tau }\right )_e\right ) = \min _{\rho _e, j_e } \max _{\phi _e} \sum _{e \in \mathcal {E}} \int _t^{t+\tau } \int _0^{L_e} \left ( \frac {\left \lvert j_e \right \rvert ^p}{p \rho _e^{p-1}} + \frac {\partial \rho _e}{\partial t} \phi _e + \frac {\partial j_e}{\partial x} \phi _e \right ) \, \mathrm d x \, \mathrm d t \end{equation*}

for $\tau \gt 0$ , $t \in \left [0, T - \tau \right ]$ and test functions $\phi _e$ , which, in case of the test function being a maximizing argument, correspond to the Lagrange multipliers. Here, the minimization is carried out among those $\rho _e, j_e$ , which fulfill the given initial and terminal values and satisfy (9). Thus, the Lagrange functional for (5) is given as

\begin{align*} L(\rho _e, j_e, \phi _e) = \sum _{e \in \mathcal {E}} \int _t^{t+\tau } \int _0^{L_e} \left ( \frac {\left \lvert j_e \right \rvert ^p}{p \rho _e^{p-1}} + \frac {\partial \rho _e}{\partial t} \phi _e + \frac {\partial j_e}{\partial x} \phi _e \right ) \, \mathrm d x \, \mathrm d t + \sum _{e\in \mathcal {E}} E_e \left (\rho _e\right ). \end{align*}

Then, the variation with respect to $j_e$ in direction $k_e$ yields

\begin{equation*} \sum _{e \in \mathcal {E}} \int _t^{t+\tau } \int _0^{L_e} \left ( \frac {j_e \left \lvert j_e\right \rvert ^{p - 2}}{\rho _e^{p-1}} k_e + \frac {\partial k_e}{\partial x} \phi _e \right ) \, \mathrm d x \, \mathrm d t = 0. \end{equation*}

Now, let us integrate by parts and use that each feasible variation $k_e$ still needs to satisfy Kirchhoff’s law, which yields

\begin{align*} 0 =& \sum _{\nu \in \mathcal {V}} \left (\sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . (k_e \phi _e) \right |_{x = L_e} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . (k_e \phi _e) \right |_{x = 0}\right ) + \\& \sum _{e \in \mathcal {E}}\int _t^{t+\tau } \int _0^{L_e} \left ( \frac {j_e |j_e|^{p-2}}{ \rho _e^{p-1}} -\frac {\partial \phi _e}{\partial x} \right ) k_e\, \mathrm d x \, \mathrm d t \\ =& \sum _{\nu \in \mathcal {V}} \left (\sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . (k_e (\phi _e-\Phi _\nu )) \right |_{x = L_e} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . (k_e (\phi _e-\Phi _\nu )) \right |_{x = 0}\right ) + \\& \sum _{e \in \mathcal {E}}\int _t^{t+\tau } \int _0^{L_e} \left ( \frac {j_e |j_e|^{p-2}}{ \rho _e^{p-1}} -\frac {\partial \phi _e}{\partial x} \right ) k_e\, \mathrm d x \, \mathrm d t \end{align*}

with the nodal mean

\begin{equation*} \Phi _{\nu } :\!= \frac {1}{\text {deg}(\nu )} \left (\sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \phi _e + \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \phi _e \right ), \end{equation*}

and the degree function

\begin{equation*} \text {deg} : \mathcal {V} \longrightarrow \mathbb {N}_0 \quad \nu \mapsto \text {deg}(\nu ) :\!= \left \lvert \left \{e \in \mathcal {E} \, \middle | \, \nu \in e\right \}\right \rvert, \end{equation*}

where for a connected graph it holds true that $\text {deg}\left (\nu \right ) \geq 1$ for all vertices $\nu \in \mathcal {V}$ .

Choosing arbitrary fluxes $k_e$ with compact support inside a single edge immediately yields

\begin{equation*} j_e \left \lvert j_e \right \rvert ^{p - 2} = \rho _e^{p - 1} \frac {\partial \phi _e}{\partial x} \end{equation*}

for each $e \in \mathcal {E}$ . Choosing further arbitrary $k_e$ with support only close to a single vertex $\nu \in \mathcal {V}$ (up to satisfying (9)), we see that further the condition

\begin{equation*} \phi _e = \Phi _\nu \qquad \forall e \in \mathcal {E} \, \text {with} \, \delta ^S(e) = \nu \text { or } \delta ^E(e) = \nu \end{equation*}

holds for the vertex $\nu$ . For dead-end vertices, i.e. deg $\left (\nu \right ) = 1$ , this yields no additional conditions, while for all other vertices this implies continuity of the dual variable.

To relate the dual variable to the energy, we can compute the variation of $L$ with respect to $\left (\rho _{t+\tau }\right )_e$ . As in the metric Wasserstein case, cf. [Reference Ambrosio, Gigli and Savaré2, Reference Santambrogio28], it turns out that the optimality condition becomes

\begin{equation*} \phi _e \left (t+\tau \right ) + E_e'\left (\rho _e\left (t+\tau \right )\right ) = 0. \end{equation*}

In the limit $\tau \downarrow 0$ , we hence obtain the following interface conditions: There exists a nodal function $\Phi : \mathcal {V} \rightarrow \mathbb {R}$ such that

(6) \begin{equation} E_e'(\rho _e)=\Phi _\nu \qquad \forall e \in \mathcal {E} \, \text {with} \, \delta ^S(e) = \nu \text { or } \delta ^E(e) = \nu, \end{equation}

for vertices $\nu \in \mathcal {V}$ . Since $\Phi _\nu$ is not a variable of physical interest, this condition effectively means the continuity of the energy variation across nodes.

To summarize, we obtain that, similar to the case of Wasserstein gradient flows in domains, the gradient flow equations are given by

(7) \begin{align} \frac {\partial \rho _e}{\partial t} + \frac {\partial j_e}{\partial x} & = 0, \end{align}

(8) \begin{align} j_e \left \lvert j_e \right \rvert ^{p - 2} & = - \rho _e^{p-1} \frac {\partial }{\partial x} \left (E'\left (\rho _e\right )\right ), \end{align}

together with the classical Kirchhoff condition for all $\nu \in \mathcal {V}$

(9) \begin{equation} 0 = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} \left . j_e \right |_{x = L_e} - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} \left . j_e \right |_{x = 0}, \end{equation}

and the coupling conditions (6).

5.2. The (ISO3) model for gas networks

In the following, we turn our attention to the (ISO3) model of transport in gas networks, which is derived and embedded in a full model hierarchy in [Reference Domschke, Hiller, Lang, Mehrmann, Morandin and Tischendorf11] (see also Section 2.1). Based on mass density $\rho _e$ and velocity $v_e$ functions defined on edges, the (ISO3) model is given by

\begin{align*} \frac {\partial \rho _e}{\partial t} + \frac {\partial }{\partial x} ( \rho _e v_e) & = 0, \\ \frac {\lambda _e}{2 \mathcal {D}_e} \rho _e v_e \left \lvert v_e \right \rvert & = - \frac {\partial }{\partial x} \left (p_e\left (\rho _e\right )\right ) - g \rho _e \sin \left (\omega _e\right ), \end{align*}

together with the classical Kirchhoff’s condition (C KL) on the vertices of the network.

We start by reformulating the model in terms of the mass density variable $\rho _e$ and the flux variable $j_e$ as:

\begin{align*} \frac {\partial \rho _e}{\partial t} + \frac {\partial j_e}{\partial x} & = 0 \\ j_e \left \lvert j_e \right \rvert & = - \frac {2 \mathcal {D}_e}{\lambda _e} \rho _e \frac {\partial }{\partial x} \left (p_e\left (\rho _e\right )\right ) - \rho _e^2 \frac {2 \mathcal {D}_e g}{\lambda _e} \sin \left (\omega _e\right ) \\ &= -\left ( \rho _e \frac {2 \mathcal {D}_e}{\lambda _e }\frac {\partial }{\partial x} \left (p_e\left (\rho _e\right )\right ) + \rho _e^2\frac {2 \mathcal {D}_e g}{\lambda _e} \sin \left (\omega _e\right )\right ) \end{align*}

Comparing to the general form of the Wasserstein gradient flow (7),(8) for $p = 3$ , we obtain the following expression for the derivative of the driving energy

\begin{equation*} \frac {\partial }{\partial x} E'(\mu )|_e = \frac {2 \mathcal {D}_e}{\lambda _e \rho _e }\frac {\partial }{\partial x} \left (p_e\left (\rho _e\right )\right ) + \frac {2 \mathcal {D}_e g}{\lambda _e} \sin \left (\omega _e\right ) \end{equation*}

on each edge $e \in \mathcal {E}$ , for some energy still to be determined for the (ISO3) model.

As usual, we need to introduce an entropy functional, which – due to the dependence on diameter $\mathcal {D}_e$ and friction coefficient $\lambda _e$ – may vary with the edge. The derivative of the entropy is defined as

\begin{equation*} F_e''(s) :\!= \frac {2 \mathcal {D}_e}{\lambda _e} \frac {p_e'\left (s\right )}s, \end{equation*}

and we define $F_e$ as a second primitive with $F_e \left (0\right ) = 0$ . Then we see that for $\rho _e \gt 0$

\begin{equation*} \frac {\partial }{\partial x} F_e'\left (\rho _e\right ) = \frac {2 \mathcal {D}_e}{\lambda _e \rho _e }\frac {\partial }{\partial x} \left (p_e\left (\rho _e\right )\right ) \end{equation*}

holds. Note that $F_e''(s) \sim \frac {p_e'(s)}s$ is a standard entropy–pressure relation arising in the overdamped limit of fluid models (cf. e.g. [Reference Giesselmann, Lattanzio and Tzavaras16]). In the case of a linear pressure $p_e(s) \sim s$ we arrive at the well-known log-entropy $F_e(s) = s \log s - s$ , but other choices of the pressure such as polynomial laws are possible as well. The essential property is convexity of the entropy, which is equivalent to $p_e$ being monotone.

With the short-hand notation

\begin{equation*} c_e:\!=\frac {2 \mathcal {D}_e g}{\lambda _e} \sin \left (\omega _e\right ), \end{equation*}

we obtain the energy functional

\begin{equation*} E(\mu ) = \sum _ e\int _0^{L_e} \left (F_e\left (\rho _e\right ) + c_e x \rho _e + d_e \rho _e\right ) \, \mathrm d x. \end{equation*}

Note that $d_e$ is an arbitrary integration constant, which however impacts the gradient flow as soon as there is more than one edge, because the mass per edge is not necessarily conserved then. Since a global change by a constant will not affect the result, due to the overall mass conservation on the network, we can however use the condition

(10) \begin{equation} \sum _{e \in \mathcal {E}} d_e = 0. \end{equation}

As the derivation above shows, the gradient flow formulation immediately yields the appropriate interface conditions (6), in addition to the classical Kirchhoff’s law for (ISO3), which are given by

\begin{align*} F_e'(\rho _e) + d_e&=\Phi _\nu \qquad \forall e \in \mathcal {E} \, \text {with} \, \delta ^S(e) = \nu, \\ F_e'(\rho _e)+ c_e L_e + d_e \ &=\Phi _\nu \qquad \forall e \in \mathcal {E} \, \text {with} \, \delta ^E(e) = \nu, \end{align*}

for vertices $\nu \in \mathcal {V}$ .

The interface condition effectively means that the potential $F_e'(\rho _e) + c_e x + d_e$ is continuous across vertices. This may seem quite arbitrary, since a potential is usually specified only up to a constant. However, we have to understand the potential rather as a global quantity on the network, so clearly changes on single edges affect the global structure of the potential. Since there is quite some freedom to determine the constants $d_e$ , we could in turn specify certain interface conditions and determine the associated potential. Again, it is reasonable that different jump conditions on the vertices have the same effect as changing a potential.

5.3. Vanishing diffusion limit

As mentioned above, there are several potentials, and thus interface conditions, that lead to gradient flows. We may want to specify however a potential such that in the limit of vanishing diffusion ( $F_e=0$ ), we obtain a consistent system. The problem of vanishing diffusion has been investigated for linear problems on networks in [Reference Egger and Philippi13] with a focus on incompressibility rather than a gradient structure.

Consistency of the vector field with the vanishing diffusion limit means

(11) \begin{align} d_e&=\Phi _\nu & \forall e \in \mathcal {E} \, \text {with} \, \delta ^S(e) = \nu, \end{align}

(12) \begin{align} c_e L_e + d_e &=\Phi _\nu & \forall e \in \mathcal {E} \, \text {with} \, \delta ^E(e) = \nu, \end{align}

for $\nu \in \mathcal {V}$ , which – together with (10) – can be considered a system of linear equations for the variables $d_e$ , or can even be extended to a linear system for the variables $d_e$ and $\Phi _\nu$ . Hence, it is natural to investigate the (unique) solvability of this system.

The number of unknowns is $\left \lvert \mathcal {E}\right \rvert + \left \lvert \mathcal {V}\right \rvert$ , while the number of equations is $\sum _{\nu \in \mathcal {V}}$ deg $\left (\nu \right ) + 1$ . By the sum formula (also known as the Handshaking-Lemma) [Reference Bondy and Murty6, Theorem 1.1], we have $\sum _{\nu \in \mathcal {V}}$ deg $\left (\nu \right ) + 1 = 2 \left \lvert \mathcal {E}\right \rvert + 1$ . Thus, the number of equations and unknowns coincides if $\left \lvert \mathcal {E}\right \rvert + 1 = \left \lvert \mathcal {V}\right \rvert$ .

It is a simple exercise to see that $\left \lvert \mathcal {E}\right \rvert \geq \left \lvert \mathcal {V}\right \rvert$ implies the existence of a cycle (cf. [Reference Bondy and Murty6, p.15]) in a so-called simple graph (a graph without loops and parallel edges). Hence, in a simple connected graph without a cycle (which seems a suitable condition for gas networks), we have $\left \lvert \mathcal {E} \right \rvert + 1 \leq \left \lvert \mathcal {V} \right \rvert$ , which is however also a lower bound for the number of edges. Thus, in this case, the number of equations and unknowns coincides and we can see that the linear system has full rank.

We finally illustrate the choice of the constants with a simple example.

Example 5.2. Let us consider a very simple tree network consisting of the vertices $\mathcal {V} = \left \{\nu _1, \nu _2, \nu _3, \nu _4\right \}$ , and the edges $\mathcal {E} = \left \{e_1, e_2, e_3\right \}$ with

\begin{align*} & \delta ^S\left (e_1\right ) = \nu _1, \qquad \delta ^S\left (e_2\right ) = \nu _2, \qquad \delta ^S\left (e_3\right ) = \nu _2, \\ & \delta ^E\left (e_1\right ) = \nu _2, \qquad \delta ^E\left (e_2\right ) = \nu _3, \qquad \delta ^E\left (e_3\right ) = \nu _4. \end{align*}

Thus, we obtain the following linear system:

Solving this system yields $d_1=-\frac {2}3 c_1 L_1$ and $d_2 = d_3 = \frac {1}3 c_2 L_2$ .

In the example we see that the constants $d_e$ on edges $e_2$ and $e_3$ only depend on the properties of edge $e_1$ , but not on $c_2, c_3, L_2, L_3$ . Indeed it is easy to see that independence of the edge parameters is always the case for an edge leading to an outgoing boundary of the network.

6.1. Discretizing the constraints

We start by discretizing the interval $[0,L_e]$ associated to every edge into an equidistant grid with $N_x^e$ points. The time interval $[0,T]$ is similarly discretized into an equidistant grid with $N_t$ points. With $\delta _x^e = L_e/N_x^e$ and $\delta _t = T/N_t$ we obtain the grid points

\begin{align*} x_i^e & = (i-1)\delta _x^e \qquad \text {for} \, e \in E, \, i \in \{1, 2, \ldots, N_x^i+1\}, \\ t_k & = (k-1)\delta _t \qquad \text {for} \, k \in \{1, 2, \ldots, N_t+1\}. \end{align*}

Note that in this section, we will slightly deviate from the notation used before, as we will use superscript $e$ and superscript $\nu$ in order to denote the corresponding edge or vertex, respectively, and not to interfere too much with the additional indices coming from the space-time discretization.

To any continuous function $v:[0,L]\times [0,T] \to \mathbb {R}$ we associate a grid function $v_h=(v_{i,k})_{i,k}$ via

\begin{align*} v_{i,k}=v(x_i,t_k) \qquad \text {for} \, i\in \{1,2, \ldots, N_x+1\}, \, k\in \{1, 2, \ldots, N_t+1\}. \end{align*}

Similar notation is used for grid functions associated to the spatial and temporal partitions, respectively. On these functions, we define a scalar product, which we obtain by discretizing the integrals involved using the trapezoidal rule. Given grid functions $u_h=(\rho _{i,k}^u,j_{i,k}^u)$ , $v_h=(\rho _{i,k}^v,j_{i,k}^v)$ , we have

(13) \begin{align} \langle u_h,v_h\rangle = \sum _{k=1}^{N_k+1}\sum _{i=1}^{N_x+1} w_k^tw_j^x\left ( \rho _{i,k}^u\rho _{i,k}^v + j_{i,k}^u j_{i,k}^v\right ),\quad \|v_h\|=\sqrt {\langle v_h,v_h\rangle }, \end{align}

with

\begin{align*} w_i^x =\begin{cases} \frac {\delta _x}{2} &i\in \{1,N_x+1\},\\ \delta _x &\text {otherwise}. \end{cases} \end{align*}

We denote by $w_k^t$ a similar weight function for the composite trapezoidal rule associated with the temporal partition and note that $\sum _{k=1}^{N_t+1}w_k^t=1$ . This renders the linear space of grid functions a Hilbert space.

We will only present the discretization of the constraints for the variant used in (TD BC 1) and (TD BC 2), with obvious modifications for the other cases. We discretize the continuity equation (CE 1) using a centred differencing in space and a backward differencing in time for the interior grid points $i\in \{2, 3, \ldots N_x\},\ k\in \{2, 3, \ldots, N_t+1\}$ , i.e.

(14) \begin{align} \frac {\rho ^e_{i,k}-\rho ^e_{i,k-1}}{\delta _t}+\frac { j^e_{i+1,k}-j^e_{i-1,k}}{2\delta _x}=0. \end{align}

For the boundary we use a one-sided finite difference approximation to approximate $\partial _x j^e$ , i.e. for $i\in \{1,N_x+1\}$ we use

(15) \begin{align} \frac {\rho ^e_{1,k}-\rho ^e_{1,k-1}}{\delta _t}+\frac { j^e_{2,k}-j^e_{1,k}}{\delta _x}&=0,\quad \text {and} \end{align}

(16) \begin{align} \frac {\rho ^e_{N_x+1,k}-\rho ^e_{N_x+1,k-1}}{\delta _t}+\frac { j^e_{N_x+1,k}-j^e_{N_x,k}}{\delta _x}&=0. \end{align}

We also add one half times a finite difference approximation of the Laplace for additional stabilization. The coupling conditions (TD BC 1) and (TD BC 2) at the boundary vertices become

\begin{alignat*} {3} && 0 & = \Bigg ( \left \lvert s^{G,\nu }_k \right \rvert + \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} j^e_{N_x+1,k} \Bigg ) - \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} j^e_{1,k} \qquad \qquad \qquad \quad & \forall \nu \in \partial ^+ \mathcal {V}, \\ \qquad && 0 & = \sum _{\substack {e \in \mathcal {E}: \\ \delta ^E(e) = \nu }} j^e_{N_x+1,k} - \Bigg ( d_k^{G,\nu } + \sum _{\substack {e \in \mathcal {E}: \\ \delta ^S(e) = \nu }} j^e_{1,k} \Bigg ) & \forall \nu \in \partial ^- \mathcal {V}, \end{alignat*}

for all $k\in \{1, 2, \ldots N_t+1\}$ , while the vertex ODEs at the interior vertices (1), become

(17) \begin{align} \frac {\gamma ^\nu _{k}-\gamma ^\nu _{k-1}}{\delta _t} = f^\nu _k, \quad \forall \nu \in \stackrel{\small\circ}{\mathcal {V}}. \end{align}

Since the discretization does not depend on $j_{i, 1}$ , we set $j_{i, 1} = 1$ for $i= 1, 2, \ldots N_x + 1$ . Furthermore, the initial and final conditions of both $\rho ^e$ and $\gamma ^{\nu }$ are taken into account in the following way

(18) \begin{align} \rho ^e_{i, 1} = \rho _0^e \left (x_i\right ) \quad \text {and} \quad \rho ^e_{i, N_t+1} = \rho _T^e \left (x_i\right ) & \qquad \qquad \text {for} \, i \in \{1, 2, \ldots N_x + 1\}, \, e \in \mathcal {E}, \end{align}

(19) \begin{align} \gamma ^{\nu }_1 = \gamma _0^{\nu } \quad \text {and} \quad \gamma ^{\nu }_{N_t + 1} = \gamma _1^{\nu } & \qquad \qquad \text {for} \, \nu \in \stackrel{\small\circ}{\mathcal {V}}. \end{align}

Finally, since the continuity equation is only satisfied approximately, we explicitly add a constraint on the total mass of the graph being a probability density

(20) \begin{align} \sum _{e\in \mathcal {E}}\sum _{i=1}^{N_x+1} w_i^x \left (\rho ^e_{i,k}-\rho ^e_0(x_i)\right ) + \sum _{\nu \in \mathcal {V}} \left (\gamma ^{\nu }_k - \gamma ^{\nu }_0\right )&=0 \qquad \text {for} \, k\in \{2, 3, \ldots N_t+1\}. \end{align}

6.2. Primal-dual algorithm

We now sketch the main ingredients of the primal-dual algorithm, referring to [Reference Carrillo, Craig, Wang and Wei9, Reference Pietschmann and Schlottbom27] for details. As mentioned above, the key idea is to enforce the constraints (14)–(20) only up to a finite precision, i.e. using the norm (13), instead of (14), we only ask for

\begin{align*} \sum _{k=2}^{N_t+1}w_k^t \Bigg ( w_1^x \left (\frac {\rho ^e_{i,k}-\rho ^e_{i,k-1}}{\delta _t}+\frac { j^e_{i+1,k}-j^e_{i-1,k}}{2\delta _x}\right )^2\Bigg ) &\leq \delta _1^2\Bigg . . \end{align*}

These weakened constraints are quadratic and can be written in the form

\begin{align*} Au \in \mathcal {C}_\delta = \left \{x \, \left .\right | \, \|x_j -b_j\|_2 \leq \delta _j, \, j = 1, \ldots, 8\right \}, \end{align*}

where the vector $u$ contains the coefficient of the grid functions $(\rho _h,m_h)$ . We have a total of $8$ constraints since we include, in addition to the continuity equation on the edges, the initial- and final data for $\rho ^e$ and for $\gamma ^{\nu }$ , the coupling conditions at the vertices, vertex dynamics and an additional mass constraint. Note that the weights $w_i^x$ and $w_k^t$ are included in the definition of $A_j$ and $b_j$ , respectively, and the vectors $x_j$ are slices of the vector $x$ corresponding to the number of rows in $A_j$ .

We define the matrix $A$ by vertically concatenating the matrices $A_j$ for $j = 1, \ldots, 8$ . We note that $A$ is the matrix of a linear map from the Hilbert space of grid functions, with the inner product being defined in (13), to a Euclidean space.

We then aim to solve the fully discretized, problem

\begin{align*} \inf _{(\rho _h,m_h,\gamma _h, f_h)} \left (\sum _{e \in \mathcal {E}}\sum _{k=1}^{N_t+1}\sum _{i=1}^{N_x+1} w^x_i w_k^t h(\rho ^e_{i,k},j^e_{i,k}) + \sum _{\nu \in {\mathcal {V}}}\sum _{k=1}^{N_t+1}h(\gamma ^{\nu }_{k},f^{\nu }_k) \right )+ \mathfrak {i}_{\mathcal {C}_\delta }(A u), \end{align*}

where $\mathfrak {i}$ denotes the convex indicator function and with $h$ as defined in (4) and $p=2$ . We now have an unconstrained optimization problem involving a convex, yet non-differentiable, functional and thus we can employ a primal-dual algorithm in the variant given in algorithm1.

Therein, two proximal operators have to be computed: The one with respect to $h$ yields, in each point in the space-time grid, a third order polynomial whose solution can be obtained by explicit formulas, while for the indicator functions we recover the projection onto the respective set, see again [Reference Carrillo, Craig, Wang and Wei9, Reference Pietschmann and Schlottbom27] for details.

Algorithm 1: Primal-dual algorithm for the solution of the optimal transport problems.

In all examples below and unless explicitly stated differently, we make the following parameter choices: $N_x = 150$ , $N_t = 75$ , $T=1$ , $L_e=1$ for all edges $e \in \mathcal {E}$ .

6.3. Example 1: Geodesics and branching

One important observation made in [Reference Erbar, Forkert, Maas and Mugnolo14] is the fact that the logarithmic entropy on metric graphs is not geodesically convex, due to the fact that geodesics may branch. This is illustrated by an explicit example, which we reproduce here numerically, both for the case with and without vertex dynamics, and for the graph topology depicted in Figure 1. As initial and final data we chose

\begin{equation*} \rho _e^0(x):\!= \begin{cases} \mathbb{1}_{[0,0.4]}(x), &e =e_1, \\ 0, & e=e_2, \, e_3,\end{cases} \quad \text { and }\quad \rho _e^1(x) :\!= \begin{cases}0, & e= e_1, \\ \frac {1}{2} \mathbb{1}_{[1-0.4,1]} (x), & e=e_2,\,e_3, \end{cases} \end{equation*}

and, in the case when we allow for a vertex dynamic,

\begin{align*} \gamma _{\nu }^0 = \gamma _{\nu }^1 = 0 \qquad \qquad \forall \nu \in \mathcal {V}. \end{align*}

The results are depicted in Figures 2 and 3 for the case without and with vertex dynamics, respectively. One can observe that with explicit vertex dynamics, the transport over the vertex is slowed down, delaying the whole transport.

Figure 1. Sketch of the graph used in the first example for branching geodesics. Here, no in- or outflux via the boundary is assumed (i.e. $\partial ^+ \mathcal {V} = \partial ^- \mathcal {V} = \emptyset$ ).

Figure 2. Branching geodesic without vertex dynamic: Snapshots of the dynamics of the densities $\rho _e$ at different times.

Figure 3. Branching geodesic with vertex dynamic: Snapshots of the dynamics of the densities $\rho _e$ and $\gamma _{\nu }$ at different times.

6.4. Example 2: Time-dependent in- and outflow

In this example, we study the influence of time-dependent boundary conditions on the solutions of the resulting transport problem with the corresponding vertex dynamics given in definition 3.1, combined with (TD BC 1) and (TD BC 2).

Therefore, we fix the graph depicted in figure 4 and set

\begin{equation*} \partial ^+ \mathcal {V} = \{\nu _1\},\quad \partial ^- \mathcal {V} = \{\nu _3,\,\nu _4\},\,\text { and }\,\stackrel{\small\circ}{\mathcal {V}} = \{\nu _2\}. \end{equation*}

Figure 4. Sketch of the graph used in the second example. We set $\partial ^+ \mathcal {V} = \{\nu _1\}$ , $\partial ^- \mathcal {V} = \{\nu _3,\,\nu _4\}$ and $\stackrel{\small\circ}{\mathcal {V}} = \{\nu _2\}$ .

Figure 5. Snapshots of the dynamics of the densities $\rho _e$ and $\gamma _{\nu }$ with symmetric boundary conditions (InOutSym) at different times.

Figure 6. Snapshots of the dynamics of the densities $\rho _e$ and $\gamma _{\nu }$ with asymmetric boundary conditions (InOutAsym) at different times.

For initial and final data we chose

\begin{align*} \rho _{e_1}^0 = \rho _{e_2}^0 =\rho _{e_3}^0 = \rho _{e_4}^0 = 0.225\quad \text { as well as }\quad \gamma _{\nu _2}^0 = 0.1,\,\gamma ^1_{\nu _3} =0, \end{align*}

and

\begin{align*} \rho _{e_1}^1 = \rho _{e_2}^1 =\rho _{e_3}^1 = 0,\, \rho _{e_4}^1(x) = c_1e^{-\frac {(x-1/2)^2}{(0.2)^2}}\quad \text { as well as }\quad \gamma _{\nu _2}^1 = 0.4,\,\gamma ^1_{\nu _3} =0, \end{align*}

where $c_1$ is chosen such that the total mass is normalized to $1$ . For in- and out-flow we consider two cases:

(InOutSym) \begin{align} s^G_{\nu _1}(t) = t\quad \text { and } \quad d^G_{\nu _3}(t) = d^G_{\nu _4}(t) = \frac {1}{2} t, \end{align}

as well as

(InOutAsym) \begin{align} s^G_{\nu _1}(t) = t,\quad d^G_{\nu _3}(t) = 0 \quad \text { and } \quad d^G_{\nu _4}(t) = t.\end{align}

Note that the values for $\gamma _{\nu _3}$ are only needed in the asymmetric case. The respective results, evaluated at different time steps, are depicted in Figure 5 (for (InOutSym)) and 6 (for (InOutAsym)), respectively. The results show that when only one vertex is an outflow vertex, the mass necessary to obtain the final configuration is predominantly transported via the other edge.