Stabilized bi-grid projection methods in finite elements for the 2D incompressible Navier-Stokes equations

1.   Introduction

Modeling, simulation and optimization of critical infrastructure systems such as traffic, electricity, water or natural gas networks play an increasingly important role in our society. Many of these problems involve aspects of dynamic energy transportation and distribution in networks. The optimization with respect to efficiency, robustness, or environmental performance requires the use of high-resolution dynamical models in form of time-dependent differential equations. As a particular case, we consider transportation and distribution of natural gas in a network of pipelines, where the high-resolution models refer to the one-dimensional Euler gas equations coupled with further dynamics representing active or passive network elements such as compressors, resistors and control valves.

Detailed models for gas flow in pipeline networks are well established, see e.g., [4,9,17]. Similar maturity has been achieved for time-simulation methods, see e.g., [1,16,29], the analysis of stationary states [22], the time-continuous optimal control of compressors using adjoint techniques [26,37], and feedback stabilization based on classical solutions [15]. If the discrete nature of control valves being open or closed is to be taken into account, then the optimization needs to deal with mixed integer and continuous type variables simultaneously. Similar decisions are often to be taken into account also in other critical infrastructure systems.

In industrial practice the treatment of switched systems is often carried out by first applying a full space-time discretization of the systems and then using a mixed-integer nonlinear programming tool that incorporates the switches as extra variables to be optimized, see e.g., [7,5,47]. Another approach is to restrict the optimization to the stationary case where special purpose techniques can be successfully applied [44]. A temporal expansion of these techniques on a full network currently seems to be out of reach, see, e.g., [23,36]. We proceed in yet another way via a model switching approach, see [24]. Since networks already provide a natural spatial partition by the edges representing pipes, it seems natural to minimize a global model error by selecting either one of several models in a dynamic model hierarchy [16,18] or a stationary model hierarchy [37] for each pipe as a function of time. To do this, it is important to identify regions in the time expanded network problem where stationary models still provide a reasonable approximation in the sense that the global error remains small. A particular difficulty that the model switching problem for gas networks has in common with most of the other mentioned applications is the high-resolution model being a partial differential equation (typically a system of balance laws). While the computation of optimal switching for ordinary differential equations and differential-algebraic equations is theoretically and numerically well studied, see [11,19,27,28,33,38,51,52,53,54], the corresponding theory involving certain types of partial differential equations is still under development [45,46].

Our main contribution in this paper is to provide a theoretical and numerical framework for solving the model switching problem using the example of gas networks. We show that the problem can be casted in the sense of switching among a family of abstract evolution equations on an appropriate Banach space. This allows us to use adjoint based gradient representations for switching time and mode sequence variations recently developed in [45] to characterize locally optimal solutions for the model switching problem by introducing the notion of first order stationarity. This, in turn, motivates a two stage gradient descent approach conceptually introduced and analyzed in [2] for the optimization of switching sequences in the context of ordinary differential equations for numerical solutions of the model switching problem. We provide results for a proof-of-concept implementation for a gas network comprising 340 km of pipes on a 30 min time horizon.

Our approach relies on a semigroup property for networked transport systems. For related results without nodal control, see [6,30,41], and boundary conditions of a delay-differential-type are considered in [48,49]. In [20,21], classical solutions for a system separating a semigroup equation from a linear nodal control condition are derived. The recent work [31] considers mild solutions if the nodal control is semilinear. In contrast, we derive a semigroup formulation for the entire system, allowing for semilinear dynamics both on the edges and the boundary control, based on a characteristic decomposition of the system.

The paper is organized a follows. In Section 2 we provide a more detailed description of the common gas network models. In particular, we briefly discuss a semilinear simplification of the Euler gas equations as the most detailed model on a pipe and consider the corresponding stationary solutions. Moreover, we introduce the model switching problem. In Section 3, we show that the equations on a network coupling these models for the pipes, along with appropriate coupling conditions on the nodes allowing further network elements such as valves and compressors, is well-posed in the sense of being equivalently represented by a nonlinear perturbation of a strongly continuous semigroup in a suitable Sobolev space. In Section 4 we apply the well-posedness result to define an appropriate system of adjoint equations and to derive a stationarity concept based on switching-time gradients and mode-insertion gradients as a first order optimality condition for model optimality. Further, we present details of a conceptual algorithm alternating between switching-time optimization and mode-insertion in order to compute a solution of the model switching problem in the sense of our stationarity concept. In Section 5 we present a numerical study. In Section 6 we discuss applications and directions of future work.

2.   Gas network modeling

2.1.   A model hierarchy for single pipes

The motion of a compressible non-viscous gas in long high-pressure pipelines is described by the one-dimensional isothermal Euler equations. They are given by a system of nonlinear hyperbolic partial differential equations (PDEs) and consist of the continuity equation and the balance of moments (see, e.g., [9,34,35,50])

where $\varrho$ denotes the density in $\frac{{{\rm{kg}}}}{{{{\rm{m}}^3}}}$, $v$ the velocity in $\frac{{\rm{m}}}{{\rm{s}}}$, $P$ the pressure of the gas in $\frac{{{\rm{kg}}}}{{{\rm{m}}\;{{\rm{s}}^2}}}$, $g$ the gravitational constant and $h'$ the slope of the pipe. Furthermore $\theta = \frac{\lambda}{2D}$, where $\lambda$ is the friction coefficient of the pipe, and $D$ is the diameter of the pipe. The conserved, respectively balanced, quantities of the system are the density $\varrho$ and the flux $q = \varrho v$. Pressure and density are related by the constitutive law for a real gas

where $Z(P, T_0)$ is the gas compressibility factor at constant temperature $T_0$ and $R_{\text{s}}$ is the specific gas constant. For an ideal gas one has $Z(P) \equiv 1$. This system of equations yields a hierarchy of simplified models for pipe dynamics [9,17,25,42].

For instationary situations, we consider the semilinear PDE model

with a constant speed of sound $c = \sqrt{P / \varrho}$ obtained under the assumption of a constant gas compressibility factor $Z(P, T_0) \equiv \bar{Z}$. This model neglects inertia effects and assumes small velocities $|v| \ll c$. For natural gas, indeed one typically has $|v| \le 10\;\frac{{\rm{m}}}{{\rm{s}}}$ and $c \approx 340\;\frac{{\rm{m}}}{{\rm{s}}}$, see e.g., [17]. This PDE exhibits two simple characteristics with speeds $\lambda_1 = -c$ and $\lambda_2 = c$.

For stationary situations, we consider the model

obtained from (2) with $\partial_t \varrho = \partial_t q = 0$. Here, the flux $q$ is constant in space and time $q(t, x) = \bar{q}$ and $\varrho(t, x) = \bar{\varrho}(x)$ with $\bar{\varrho}$ being a solution of the momentum equation in (3), which is a Bernoulli-equation. A solution of the momentum equation can therefore be obtained algebraically

Analogous considerations apply also for the case of non-isothermal flow, see e.g. [9,17]. More generally, similar model hierarchies also exist for other infrastructure systems such as water distribution networks [24].

2.2.   Networks with pipes, valves, and compressors

For $m, n \in \mathbb{N}$ we consider a network of pipes that we model by a metric graph $G = (V, E)$ with nodes $V = (v_1, \ldots, v_m)$ and edges $E = (e_1, \ldots, e_n) \subseteq V \times V$. For each edge $e \in E$, call $e(1)$ the left node and $e(2)$ the right node of $e$. We demand the incident nodes of every edge to be different, so $e(1) \neq e(2)$ for any $e \in E$ and thus self-loops are not allowed. On the other hand, if $v \in V$ is any node, then we define

The number $|\delta v|$ then is called the degree of node $v \in V$.

With each edge $e_j \in E$ of such a network, we associate a pipe model from the hierarchy in Section 2.1 and a given pipe length $L_j >0$. Furthermore, depending on the role of each node in the network, we impose appropriate coupling conditions for the gas density and flow at the boundary of pipes corresponding to edges being incident to that, see [13]. To this end, we define for $v \in V$ and $e_j \in \delta v$

For each node $v \in V$, we then impose a transmission condition for the density and a balance equation for the fluxes at the node. The transmission condition states that the density variables $\varrho^j$ weighted by given factors $\alpha \in (0, \infty)^{m \times 2}$ coincide for all incident edges $e \in \delta v$ and can be expressed as

The nodal balance equation for a given outflow function $q^v\colon [0, T] \to \mathbb{R}$ is similar to a classical Kirchhoff condition for the fluxes $q^j$ and can be written as

The choice of $\alpha$ corresponds to the nodal types in the network. Prototypically, we will consider the following cases:

Junctions: nodes $v$ such that $q^v \equiv 0$ and $\alpha^k_{x(v, e_k)} = 1$ for all $e_k \in \delta v$.

Boundary nodes: nodes $v$ such that $\alpha^k_{x(v, e_k)} = 1$ for all $e_k \in \delta v$, but $q^v \not\equiv 0$.

We also refer to $v$ as an entry node, if $q^v <0$, or an exit node, if $q^v>0$.

Compressors: nodes $v$ with $q^v \equiv 0$ and $|\delta^+ v| = |\delta^- v| = 1$. A description established via the characteristic diagram based on measured specific changes in adiabatic enthalpy $H_{\text{ad}}$ of the compression process yields the model

where $\kappa$ is a compressor specific constant, $\bar{Z}$ is the gas compressibility factor that is assumed to be constant and $H_{\text{ad}}$ is within flow dependent and compressor specific bounds obtained from the characteristic diagram. In consistency with the pipe models, we assume that $H_{\text{ad}}$ is given by a known reference $\bar H_{\text{ad}}$. Then we get

with a compressor specific factor

This yields $\alpha^k_1 = 1$ and $\alpha^l_0 = \bar{\alpha}$. The compression ratio of centrifugal and turbo compressors mostly lies within the range $\bar{\alpha} \in [1, 2]$, with $\bar{\alpha} \approx 1.3$ being a typical value, see [13], [39,Chapter 4] and [12,Chapter 5.1]. For realistic long distance gas networks, compressors are set along the pipes in a distance of around $100$-$300\ \text{km}$, see [12,Chapter 5.3].

We note that there are other conceivable ways to model gas pipe junctions, including geometry conditions, leading to possibly nonlinear coupling conditions, see e.g., [3,9,40]. Furthermore, gas networks typically involve additional network components such as valves, resistors, gas coolers, etc. For appropriate coupling conditions we refer to [44].

2.3.   The model switching problem

We now discuss switching between the two models (2) and (3) in order to efficiently resolve the dynamics of the gas flow in a network. The idea is that, with the exception of locally high fluctuation, in realistic scenarios, the solution to (2) is on big parts of the network close to the stationary model (3). In these regions we thus can freeze the solution with an acceptable loss in accuracy to save computational effort. By comparison with the solution fully simulated with (2) we then can set up a cost functional measuring the deviation of the partially frozen solution in some appropriate norm. Adding a performance function measuring the time steps where the costly model (2) is calculated, we can set up the optimization problem of weighting the accuracy against the computational effort. Solving this problem enables us to identify a time-dependent model selection for the simulation of gas dynamics on networks that can be used in further examinations to get a cheaply solvable system.

In order to apply adjoint-based gradient methods to solve this problem efficiently, we consider (3) as a limit of appropriately perturbed instationary semilinear dynamics. To this end, we note that, after a linear transformation, system (2) can be written for $t \in [0, T]$ and $x \in [0, L]$, $L, T>0$, as partially diagonalized system

with the characteristic variable $w$, see Section 3 for details. Let us assume for the moment that $g\colon \mathbb{R} ^2 \to \mathbb{R} ^2$ is globally Lipschitz-continuous. The same transformation applied to (3), together with given boundary values $w^1\in C([0, T], \mathbb{R} ^2)$ yields the ODE-system

To define an appropriate hyperbolic system to approximate (6), for $\varepsilon>0$ and any initial values $w^0 \in L^1([0, L], \mathbb{R} ^2)$ we set up the auxiliary system

Referring to [8,Chapter 3.4] for the definition of a broad solution to system (7), we get the following result:

Lemma 2.1. For any $\varepsilon>0$ there is a unique solution $\bar{w} \in C([0, T] \times [0, L], \mathbb{R} ^2)$ to (6) and a unique broad solution $w\colon [0, T] \times [0, L] \to \mathbb{R} ^2$ to (7). The restriction of $w$ to the set

is a continuous function. Furthermore, for any $t_0 \in (0, T)$ the function $w$ converges uniformly on $[t_0, T]\times[0, L]$ for $\varepsilon \searrow 0$ to $\bar{w}$, i.e.

Proof. The global existence and uniqueness of $\bar{w}$ follows from standard ODE-theory. For $(t, x) \in \Omega$, the broad solution $w$ can be written as

Its existence and the continuity on $\Omega$ is shown in [8,Theorem 3.1,Theorem 3.4]. Given a fixed $t_0 \in (0, T)$, we obviously can choose $\varepsilon$ small enough such that $[t_0, T] \times [0, L] \subseteq \Omega$. We then have

for $\varepsilon \searrow 0$ and all $(t, x) \in [t_0, T] \times [0, L]$. The convergence is uniform, since both $w^1_2$ and $g_2(w)$ are continuous and defined on a compact set, thus uniformly continuous. The same argumentation holds for $w_1$.

Retransforming (7) back to the original variables $\varrho$ and $q$, we get

By Lemma 2.1, this system yields an approximation to model (3).

Now let $G = (V, E)$ be a network, where the type of each node $v \in V$ is given by the parameters $\alpha$ and $q^v$ as in Section 2.2. On each edge $e_j$ of length $L_j$ we have an initial gas density $\varrho^j_0$ and gas flow $q^j_0$. For an increasing sequence of switching times $\tau = (\tau_k)_{k = 0, \ldots, N+1} \subseteq [0, \infty)$ and a finite sequence $\mu = (\mu_k)_{k = 1, \ldots, N}$, where $\mu_k \in \{0, 1\}^n$ denotes for each edge $e_j \in E$, whether model (2) ($\mu_k(j) = 1$) or (8) ($\mu_k(j) = 0$) is used on the time interval $(\tau_{k-1}, \tau_k)$, consider the PDE-system

where $x \in [0, L_j]$ for $j = 1, \ldots, n$ and $t \in (\tau_{k-1}, \tau_k)$ for $k = 1, \ldots, N+1$. Here,

for a fixed $\bar{\varepsilon}>0$ with $\bar{\varepsilon} \ll 1$ and

denote the right-hand sides in (8) and (2), respectively. Let $\varrho^j(0, x) = \varrho^j_0(x)$, $q^j(0, x) = q^j_0(x)$ for $x \in [0, L_j]$ be given initial states. Moreover, at any time point $t \in [\tau_0, \tau_{N+1}]$, density and flux additionally satisfy the node coupling conditions

for each $v \in V$, explained in Section 2.2.

Denote by $z_d = (\varrho^1_d, q^1_d, \ldots, \varrho^n_d, q^n_d)^\top$ the reference solution to (9), (11) for the choice $N = 1$, $\mu = 1$ and $\tau = (0, T)$, which corresponds to the fine model (2) being fully solved on the complete network and the existence of which will be proven in Section 3. For any other $z = (\varrho^1, q^1, \ldots, \varrho^n, q^n)^\top$ then define the cost functional

with $\gamma_1, \ldots, \gamma_4 \geq 0$, where the first term measures the deviation of $z$ from $z_d$, the second term penalizes using the fine model (2) and the third term penalizes the number of switching time points. Note that, since longer pipes mean more computational effort when using the fine model, the lengths $L_j$ of the pipes enter into the cost as well. For later reference, we set

then $J = J_1+J_2$.

The challenge now is to choose the sequences $\mu$ and $\tau$ such that, with $z$ being the corresponding solution to (9), (11), the cost functional $J$ is minimized. Hence our objective is to solve the minimization problem

3.   Semigroup formulation on networks

In this section, we will set up an abstract formulation of the PDE-system in (9) for $\tau = (0, T)$ and any fixed choice of modes per edge and prove that the linear part generates a $C^0$-semigroup. By nonlinear perturbation theory and induction, this yields a well-posedness result of (9) for any finite switching sequence.

For each $j \in \{1, \ldots, n\}$, we now consider the initial boundary value problem for $z^j = (\varrho^j, q^j)^\top$ on pipe $e_j \in E$ given by

and the coupling conditions

for each node $v \in V$. Here, $L_j>0$, $z^j_0\colon [0, L_j] \to \mathbb{R} ^2$, $f^j \in C^1(\mathbb{R} ^2, \mathbb{R} ^2)$ is globally Lipschitz-continuous and

for each $j \in \{1, \ldots, n\}$. Moreover, $\alpha \in (0, \infty)^{m \times 2}$ and $q^v_0\colon [0, \infty) \to \mathbb{R}$ is a given outflow for each $v \in V$. Introduce the space

the vectors

the operators

Moreover, define the operator

on the domain

of all vectors of absolutely continuous functions that satisfy the boundary conditions (16) and (17). Note that, though the inflow functions $q^v$ is only evaluated at the origin, it also is shifted in time due to the transport-type evolution represented by operator $B$ - this together enforces the originally time-dependent coupling condition (17). We then can reformulate system (15)-(17) as the abstract initial-value problem

Now we can state the following result:

Theorem 3.1. The operator

is the infinitesimal generator of a strongly continuous semigroup on $Z$.

Proof. Note that $\text{diag}(A, B)$ is a densely defined, closed operator on $Z$ with a nonempty resolvent set (for instance $0 \in \rho(\text{diag}(A, B))$). Following [43,Chapter 4,Theorem 1.3], to prove the claimed semigroup property, it thus suffices to show that the homogeneous system

has a unique solution $z \in C^1([0, T], Z)$ for any $T>0$ and every choice $z_0 \in D(\text{diag}(A, B))$. So let

be given and first assume $T \leq \bar{T} : = \min\left\{ \frac{\varepsilon_jL_j}{2c_j} \, | \, j = 1, \ldots, n \right\}>0$. We recognise in the operator $B$ a transport equation with velocity $-1$ for the variables $q^{v_1}, \ldots, q^{v_m}$ with the well-known unique solution

Next, for each note $v \in V$ we can set

and

Since $\varrho^1, q^1, \ldots, \varrho^n, q^n, q^{v_1}, \ldots, q^{v_m}$ are absolutely continuous, so is $\varrho^{v}$ for each $v \in V$. For each edge $e_j \in E$, $j = 1, \ldots, n$, construct the functions $\varrho^j, q^j\colon [0, T] \times [0, 1] \to \mathbb{R}$ as follows: for $t \in (0, T]$ set

Substituting the coupling conditions stated in (19), one can indeed show that, with these definitions,

We skip these calculations here for brevity. Next, set

for $(t, x) \in (0, T]\times(0, L_j)$, where

The above considerations show that $(\varrho^j_*, q^j_*)$ is continuous and piecewise absolutely continuous, thus absolutely continuous everywhere. To see that these functions in fact solve system (21), first note that the initial condition and the coupling condition (16) are satisfied by construction. Moreover, for each $v \in V$ we have

Observe that $\varrho^j$ and $q^j$ are continuous at the boundaries $x = 0$ and $x = L_j$ and that

for a.e. $(t, x) \in (0, T]\times(0, L_j)$, thus the above construction indeed yields a solution to (21) for $T \leq \bar{T}$ with $z(t) \in D$ for $t \in [0, T]$. For the case $T > \bar{T}$ note that the same steps can be repeated successively to expand the solution for the times $[\bar{T}, 2\bar{T}], [2\bar{T}, 3\bar{T}]$ and so on. To prove uniqueness of the solution, consider the case $z_0 = 0$ of initial data and introduce the energy

Obviously, we have $E(t) = 0$ if and only if $\varrho^j(t, .) = q^j(t, .) = 0$ for all $j = 1, \ldots, n$, $E(t) \geq 0$ for all $t\geq 0$ and $E(0) = 0$. Furthermore,

since $\varrho^v \equiv 0$ for all $v \in V$. But then $E \equiv 0$, hence $z \equiv 0$. This concludes the proof for the semigroup property.

We now can apply standard arguments from semigroup theory for further discussion of (20): by [43,Chapter 6,Theorem 1.2 and 1.5], if the inhomogeneity $f$ is continuously differentiable and globally Lipschitz-continuous, (20) has a unique classical solution $z$. Though these assumptions do not apply to the friction term $f$ stated in (2), we can use the following technical consideration to still get a unique solution: choose any $\bar{\varrho}, \bar{q}>0$ and define the modified right-hand side

The function $\tilde{f}$ is continuously differentiable and globally Lipschitz-continuous as a function mapping $L^2([0, L], \mathbb{R} )$ onto itself for any $L>0$. If we replace $f$ by the $\tilde{f}$ in (9), we then have a system fitting the assumptions made above for each fixed $k = 1, \ldots, N$ and thus we have a unique solution. If we choose $\bar{\varrho}$ sufficiently small and $\bar{q}$ sufficiently large, then the numerical simulation with realistic data shows that both the boundaries $\bar{\varrho}$ and $\bar{q}$ will in fact never be reached by $\varrho$ and $q$, respectively. In this case, however, the solution coincides with the unique solution to the original problem (9). The same argumentation holds for the right-hand side in (8). In summary, this proves the necessary results on well-posedness of the system (20) that we need to examine the optimization problem (14).

4.   Optimality conditions and a solution method

By the results developed in Section 3, we find that the optimization problem (14) is of the form

where $A^{\mu_k}$ is a strongly continuous semigroup, $f^{\mu_k}$ is a semilinear perturbation and $g^{\mu_k, \mu_{k+1}} = \text{id}$ is a transition map that can be chosen trivially here for each $k = 1, \ldots, N$. The ordering cone $\mathcal{T}(0, T)$ is for fixed $T>0$ defined by

In general, problem (22) does not have a unique global minimum. Instead, we will define a concept of stationarity fitted to the hybrid nature of the problem and set up an algorithm capable of finding such stationary points.

First, however, we apply on system (22) some results from the preliminary work in [45]: by Theorem 3.1 and the subsequent remarks, [45,Lemma 2] can be applied, yielding a control-to-state-map $(\mu, \tau) \mapsto z(\mu, \tau)$ that can be substituted into $J$ to define the reduced cost functional $\Phi(\mu, \tau) = J(\mu, \tau, z(\mu, \tau))$. We can now apply [45,Theorem 8] to show that $\Phi$ is continuously differentiable with respect to the switching time $\tau_k$. Recalling the shortened notation introduced in (13), we find that, while $J_2$ can be differentiated directly with respect to $(\mu, \tau)$, the term $J_1$ depends on the solution $z = z(\mu, \tau)$. Again by [45,Theorem 8], we can state gradient formulae for the derivatives of $J_1$ based on the adjoint equation

Here, $f^\mu_z\colon Z \to L(Z)$ and $l_z\colon Z \to Z^*$ denote the derivatives of the friction term $f^\mu$ and $l$ as in (13) with respect to $z$. Substituting the definitions made in (18) yields that in fact $\gamma \equiv 0$ and $p = (p^j)_{j = 1, \ldots, n}$ again can be partitioned into edgewise defined functions $p^j = (p^j_1, p^j_2)^\top$ for $j = 1, \ldots, n$ satisfying

where

are the derivatives of $f^0$ and $f^1$ in (10), as well as the coupling conditions

By [45,Lemma 6], system (23) (and thus (24), (25), (26) as well) have a unique mild solution and, applying [45,Theorem 8], the switching time gradient is given by

If the mode sequence $\mu$ is fixed, then we can conclude that a switching sequence $\tau = (\tau_k)_{k = 0, \ldots, N+1}$ is a KKT-point of the minimization problem (14), if the following holds: For $n \in \{1, \ldots, N\}$ set $a(\tau, n) = \min \{ m \in \{0, \ldots, n\} \, | \, \tau_m = \tau_n \}$ and $b(\tau, n) = \max \{ m \in \{n, \ldots, N+1\} \, | \, \tau_m = \tau_n \}$, then

Similarly, by [45,Theorem 10], the sensitivity of the cost function with respect to introducing a new mode $\mu'$ on an infinitesimal time interval at the point $\tau'$ can be represented by the mode insertion gradient given by

where $\mu_k$ is the original mode at time $\tau'$. In summary, a switching sequence $(\mu, \tau)$ is called stationary, if $\tau$ is a KKT-point for $\mu$ fixed and if $\frac{\partial \Phi}{\partial \mu'}(\tau') \geq 0$ for all modes $\mu'$ and all times $\tau' \in [0, T]$. Any global minimum of the problem (14) then is stationary.

In order to compute such stationary switching signals, we consider a conceptual algorithm originally proposed for optimal mode scheduling in hybrid ODE-dynamical systems [2]. The main idea is a two phase approach as follows: the algorithm is initialized with a switching sequence $(\mu^0, \tau^0)$, for instance $\mu^0 = 1$ and $\tau^0 = (0, T)$ where no switching occurs and the system is solved by keeping the mode constant at 1.

In a first phase, the positions of available switching time points are optimized, while conserving their order, by using a projected-gradient method with Armijo step size. To this end a projection $P$ onto the ordering cone $\mathcal{T}(0, T)$ is used. In a second phase, a new mode $\mu'$ is inserted at a time point $\tau'$ where $\frac{\partial \Phi}{\partial \mu'}(\tau') < 0$. If no such point exists, the switching sequence is stationary in the above sense, otherwise the algorithm continues with the first phase again.

A more precise description of the procedure is given in Algorithm 1. Note that Algorithm 1 does not necessarily terminate with the global solution to the minimization problem (22), but with a stationary point in the above sense. A convergence analysis for this algorithm can be found in [2]. A possible implementation of the projection $P$ can be found in [19].

5.   Numerical study

As a proof of concept we consider a gas network outlined in Figure 1b. At the boundary node $N_1$ gas is supplied to the network whereas nodes $N_2$ and $N_3$ are customer nodes where gas can possibly be taken out of the network. This can be seen as a simplified example of a big regional gas pipeline network with a local distribution subnetwork. In the particular scenario we are looking at, there is gas transported from the supply $N_1$ to node $N_2$ to satisfy a given demand while $N_3$ is inactive. All pipes are assumed to be horizontal ($h' = 0$), the compressor is assumed to be running at a constant compression factor $\bar{\alpha}$, compare (5), and at initial time the network is assumed to be stationary with constant densities $\varrho_0$ on the outer circle (thus $\alpha \varrho_0$ on the inner circle) and flux $q_0$ everywhere, moreover the outflows at the nodes $N_1$ and $N_2$ are outlined in Figure 1c. See the table in Figure 1a for specific values for those and other constants.

Due to the almost decoupled inner and outer circle, it can be expected that the numerical solution highly varies on the outer circle connecting $N_1$ and $N_2$ but is near to constant on the inner circle. This is confirmed by the simulation of the full model, see Figure 2 for a snapshot of the simulation at the time $t = 900 \ \text{s}$. We therefore can suspect that the simulation on the inner circle can be widely frozen with only some small losses in the accuracy of the solution. Indeed, letting $\bar{z}$ be the distinguished solution to (9), (11) resulting from freezing the solution completely on the inner edges $6$ to $10$, our simulations show that the $L^2$-errors of the density and the flux relative to the respective maximum values in $z_d$ is less than $1\%$ for $\varrho$ and less than $6.5\%$ for $q$; see Figure 3d. Moreover, compared to a full simulation $z_d$ about half of the computation time in the sense of $J_2$ defined in (13) is saved.

In our implementation the system (9), (11) is solved via splitting of the hyperbolic part and the friction term. For the simulation of the hyperbolic part we use the 2-step-Richtmyer-method with artificial viscosity, see [32,Chapter 18.1]. Given a system matrix $A$ and a discretization $z^k = (z^k_1, \ldots, z^k_n)^\top$ with spatial step size $\Delta x$ of the solution at time point $t_k$, it computes the discretized values $z^{k+1}$ at time $t_{k+1} = t_k+\Delta t$ by the explicit finite-volume-scheme

for $j = 1, \ldots, n$, where $\varepsilon \sim 0.05$ introduces an additional, minor but stabilizing smoothing to the solution. The discretization is chosen in a way such that the cells $z^k_1$ and $z^k_n$ are centered at the boundary points $x = 0$ and $x = 1$, respectively. In order to handle the boundary values and the inner cells simultaneously with the same scheme, we add appropriate ghost cells on each end of the spatial domain. The appropriate ghost cell values $(z^e)^k_0, (z^e)^k_{n+1}$ for each $e \in E$ are derived from the coupling conditions. We mention without further proof that these can be solved explicitly for each node $v \in V$ by first setting the weighted mean value

and then using the ghost cell values

Here, in each time step, we incorporate only those edges where the solution is actively calculated. In the special case if $v \in V$ is a compressor, we only have one ingoing edge $e_+ \in \delta^+v$ and one outgoing edge $e_- \in \delta^-v$ and instead have to set

For the friction term, we add an explicit Runge-Kutta-step using the classical Runge-Kutta-scheme of order 2 or midpoint rule, see [10,Chapter 23]. The same methods are used on the same discretization grid backwards in time to solve the adjoint system stated in (24), (25), (26).

To evaluate the gradient formulae (27) and (29), we use the trapezoidal rule over the spatial grid, see [14,Chapter 2.1]. The expression $[(A^{\mu_k})^j-(A^{\mu_{k-1}})^j] z^j(\tau_k, x)$ occurring in both (27) and (29) represents the difference of the time derivatives of the solution depending on which mode $\mu_k$ is switched to. In the numerical scheme, this is realized by calculating a time step of the solution for each choice of $\mu_k$ and then substracting the forward difference quotients.

The fixed temporal discretization grid for the actual solution is supplemented by a grid of switching time points that may vary due to the superordinated optimization where the data needed for the gradient formulae is calculated. For our study we start the optimization with the fully frozen solution, where in no time step the model (2) is actually calculated, and iteratively insert regions of active numerical solving wherever the gradients indicate a major loss of accuracy. Applying the projected-gradient method with sequential mode insertion described in Algorithm 1 yields the results shown in Figure 3a. We observe that it is in fact almost unnecessary to calculate the fine model on the edges 6 to 10 of the inner circle and that the algorithm indeed approximates the distinguished solution $\bar{z}$ as expected. Note that Algorithm 1 does not remove switching points during its process, because doing so might lead to the iteration being caught in an infinite loop, where the same switching points are added and removed repeatedly. This way, however, there remain scattered short intervals where the mode is switched twice almost instantly, which can be interpreted as the algorithm eliminating the respective intermediate modes. We suggest a post-processing step to remove this scattering, for instance by applying the following filtering rule with $d \in \mathbb{N}$ a fixed parameter: if within $d$ time steps after a switching point in the numerical solution the mode is switched again, then both switching time points are removed. Relative to the initial cost, we get for the optimized solution a cost reduction of $96.4\%$ with $65$ switching points, for $d = 5$ shown in Figure 3b a cost reduction of $96.3\%$ with $19$ switching points and for $d = 10$ shown in Figure 3c a cost reduction of $95.1\%$ with $7$ switching points.

6.   Conclusion

We have presented an application of the theory of switching systems to a model hierarchy for the dynamics of gas in a pipeline network. A semigroup formulation was given for the model on gas networks including time-dependent outflow at each node as well as a linearized model for compressors that allowed us to prove the existence of unique solutions. Using adjoint-based gradient representations for switching abstract evolution systems, we implemented a two stage projected-gradient descent method for the optimization of switching between different levels of accuracy in the hierarchy. As an example, we optimized the simulation of a small supply network by freezing the calculation on edges whenever the numerical error made is small compared to the computational costs.

Our prototypical approach can be applied in a similar fashion to realistic industrial networks. The technique can also be extended to identify a model switching strategy for a reduced model using a range of different parameter configurations such as representing different boundary flow scenarios, compressor settings and valve positions. The resulting reduced simulation model can then be used for a time expanded mixed-integer optimization technique based on full discretization with a minimum of variables or within a bilevel optimisation method which optimises a cost functional on the outer level with optimal efficiency on the lower level as proposed, for example, in [24]. Further directions include a combination of this method in a network of submodel hierarchies such as those used in [37] for the optimisation of stationary models. Moreover, it remains an open question, whether the chosen numerical method lead to a consistent discretization of the optimality system. Well-balanced finite-volume schemes also seem to be an promising alternative approach in this context.

Acknowledgments

The authors thank the Deutsche Forschungsgemeinschaft for their support within the projects A03 and B03 in the Sonderforschungsbereich/-Transregio 154 Mathematical Modelling, Simulation and Optimization using the Example of Gas Networks. Furthermore, the authors are indebted to the anonymous reviewers who considerably improved the exposition of this paper.