Optimal model switching for gas flow in pipe networks

  • Received: 01 January 2018 Revised: 01 August 2018
  • Primary: 35L50, 76N25, 93C30; Secondary: 35R02, 49J20

  • We consider model adaptivity for gas flow in pipeline networks. For each instant in time and for each pipe in the network a model for the gas flow is to be selected from a hierarchy of models in order to maximize a performance index that balances model accuracy and computational cost for a simulation of the entire network. This combinatorial problem involving partial differential equations is posed as an optimal switching control problem for abstract semilinear evolutions. We provide a theoretical and numerical framework for solving this problem using a two stage gradient descent approach based on switching time and mode insertion gradients. A numerical study demonstrates the practicability of the approach.

    Citation: Fabian Rüffler, Volker Mehrmann, Falk M. Hante. Optimal model switching for gas flow in pipe networks[J]. Networks and Heterogeneous Media, 2018, 13(4): 641-661. doi: 10.3934/nhm.2018029

    Related Papers:

    [1] Fabian Rüffler, Volker Mehrmann, Falk M. Hante . Optimal model switching for gas flow in pipe networks. Networks and Heterogeneous Media, 2018, 13(4): 641-661. doi: 10.3934/nhm.2018029
    [2] Michael Herty . Modeling, simulation and optimization of gas networks with compressors. Networks and Heterogeneous Media, 2007, 2(1): 81-97. doi: 10.3934/nhm.2007.2.81
    [3] Michael Herty, Veronika Sachers . Adjoint calculus for optimization of gas networks. Networks and Heterogeneous Media, 2007, 2(4): 733-750. doi: 10.3934/nhm.2007.2.733
    [4] Mapundi K. Banda, Michael Herty, Axel Klar . Gas flow in pipeline networks. Networks and Heterogeneous Media, 2006, 1(1): 41-56. doi: 10.3934/nhm.2006.1.41
    [5] Markus Dick, Martin Gugat, Günter Leugering . Classical solutions and feedback stabilization for the gas flow in a sequence of pipes. Networks and Heterogeneous Media, 2010, 5(4): 691-709. doi: 10.3934/nhm.2010.5.691
    [6] Martin Gugat, Falk M. Hante, Markus Hirsch-Dick, Günter Leugering . Stationary states in gas networks. Networks and Heterogeneous Media, 2015, 10(2): 295-320. doi: 10.3934/nhm.2015.10.295
    [7] Yue Qiu, Sara Grundel, Martin Stoll, Peter Benner . Efficient numerical methods for gas network modeling and simulation. Networks and Heterogeneous Media, 2020, 15(4): 653-679. doi: 10.3934/nhm.2020018
    [8] Yogiraj Mantri, Michael Herty, Sebastian Noelle . Well-balanced scheme for gas-flow in pipeline networks. Networks and Heterogeneous Media, 2019, 14(4): 659-676. doi: 10.3934/nhm.2019026
    [9] Christian Budde, Marjeta Kramar Fijavž . Bi-Continuous semigroups for flows on infinite networks. Networks and Heterogeneous Media, 2021, 16(4): 553-567. doi: 10.3934/nhm.2021017
    [10] Mapundi K. Banda, Michael Herty, Axel Klar . Coupling conditions for gas networks governed by the isothermal Euler equations. Networks and Heterogeneous Media, 2006, 1(2): 295-314. doi: 10.3934/nhm.2006.1.295
  • We consider model adaptivity for gas flow in pipeline networks. For each instant in time and for each pipe in the network a model for the gas flow is to be selected from a hierarchy of models in order to maximize a performance index that balances model accuracy and computational cost for a simulation of the entire network. This combinatorial problem involving partial differential equations is posed as an optimal switching control problem for abstract semilinear evolutions. We provide a theoretical and numerical framework for solving this problem using a two stage gradient descent approach based on switching time and mode insertion gradients. A numerical study demonstrates the practicability of the approach.



    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.

    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])

    $ tϱ+x(ϱv)=0,t(ϱv)+x(P+ϱv2)=θϱv|v|gϱh,
    $
    (1)

    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

    $ P = R_{\text{s}} \varrho T_0 Z(P, T_0), $

    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

    $ tϱ+xq=0,tq+c2xϱ=θq|q|ϱghϱ,
    $
    (2)

    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

    $ xq=0,c2xϱ=θq|q|ϱghϱ,
    $
    (3)

    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

    $ \bar{\varrho}(x) = {ˉϱ(0)22θˉq|ˉq|c2x,if h=0ˉϱ(0)2exp(2ghx)2θˉq|ˉq|c2exp(2ghx)12gh,if h0.
    $
    (4)

    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].

    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 set of ingoing edges byδ+v={eE|e(2)=v},the set of outgoing edges byδv={eE|e(1)=v},the set of incident edges byδv=δvδ+v.
    $

    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$

    $ x(v, e_j) = {0, if ejδv,1, if ejδ+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

    $ \alpha^k_{x(v, e_k)}\varrho^k(t, L_k x(v, e_k)) = \alpha^l_{x(v, e_l)} \varrho^l(t, L_l x(v, e_l)), \ \ \ \ \forall \, e_k, e_l \in \delta v, \, t \in [0, T]. $

    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

    $ \sum\limits_{e_j \in \delta^+v} q^j(t, L_j) - \sum\limits_{e_j \in \delta^-v} q^j(t, 0) = q^v(t), \ \ \ \ \ t \in [0, T]. $

    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

    $ H_{\text{ad}} = \bar{Z}T_0R_s \frac{\kappa}{\kappa-1}\left( \left(\frac{\varrho^l(0, t)}{\varrho^k(L_k, t)}\right)^{\tfrac{\kappa-1}{\kappa}}-1 \right), \ \ \ \ \ e_k \in \delta^+ v, e_l \in \delta^- v, \, t \in [0, T], $

    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

    $ \varrho^l(0, t) = \bar{\alpha} \varrho^k(L_k, t),\ \ \ \ e_k \in \delta^+ v, e_l \in \delta^- v, \, t \in [0, T] $

    with a compressor specific factor

    $ \bar{\alpha} = \left(1+ \frac{(\kappa-1)\bar H_{\text{ad}}}{\kappa \bar{Z} T_0 R_s} \right)^{\tfrac{\kappa}{\kappa-1}}. $ (5)

    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].

    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

    $ w_t(t, x) +Dw_x(t, x) = g(w(t, x)), \ \ \ D = [c00c]
    , \ \ \ (t, x) \in (0, T) \times (0, L), $

    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

    $ Dwx(t,x)=g(w(t,x)),(t,x)(0,T)×(0,L),w1(t,1)=w11(t),t[0,T],w2(t,0)=w12(t),t[0,T].
    $
    (6)

    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

    $ εwt(t,x)+Dwx(t,x)=εD1g(w(t,x)),(t,x)(0,T)×(0,L),w(0,x)=w0(x),x(0,L),w1(t,1)=w11(t),t[0,T],w2(t,0)=w12(t),t[0,T].
    $
    (7)

    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

    $ \Omega = \{(t, x) \in [0, T]\times[0, L] \, | \, L-\tfrac{c}{\varepsilon}t < x < \tfrac{c}{\varepsilon}t \} $

    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.

    $ \lim\limits_{\varepsilon \searrow 0} \sup\limits_{(t, x) \in [t_0, T]\times[0, L]} |w(t, x)-\bar{w}(t, x)| = 0. $

    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

    $w1(t,x)=w11(tεc(Lx))+ttεc(Lx)1cg1(w(s,x+cε(ts)))ds,w2(t,x)=w12(tεcx)+ttεcx1cg2(w(s,cε(ts)))dy.
    $

    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

    $w2(t,x)=w12(tεcx)+ttεcx1cg2(w(s,cε(ts)))dyw12(t)+x01cg2(w(s,y))dy=ˉw2(t,x)
    $

    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

    $ ϱt+1εqx=θc2q|q|ϱghc2ϱ,qt+c2εϱx=0.
    $
    (8)

    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

    $ [tϱj(t,x)tqj(t,x)]
    +\frac{1}{\varepsilon_{\mu_k(j)}} [01c20]
    [xϱj(t,x)xqj(t,x)]
    = f^{\mu_k(j)}\left([ϱj(t,x)qj(t,x)]
    \right) $
    (9)

    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,

    $ \varepsilon_{\mu_k(j)} = {ˉε, if μk(j)=0,1, if μk(j)=1
    $

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

    $ f^0(\varrho, q) = [θc2q|q|ϱghc2ϱ0]
    \ \ \ \text{and} \ \ \ f^1(\varrho, q) = [0θq|q|ϱghϱ]
    $
    (10)

    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

    $ αix(v,ei)ϱi(t,Lix(v,ei))=αjx(v,ej)ϱj(t,Ljx(v,ej)),    ei,ejδv,ejδ+vqj(t,Lj)ejδvqj(t,0)=qv(t)
    $
    (11)

    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

    $ J(μ,τ,z)=nj=1T0Lj0γ1(ϱj(t,x)ϱjd(t,x))2+γ2(qj(t,x)qjd(t,x))2dxdt=+γ3Nk=1nj=11Ljτk+1τk(εμk(j)ˉε)2dt,=+γ4N
    $
    (12)

    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

    $ l(t,z)=nj=1Lj0γ1(ϱj(t,x)ϱjd(t,x))2+γ2(qj(t,x)qjd(t,x))2dx,J1=12T0l(t,z)dt,J2=γ3Nk=1nj=11Ljτk+1τk(εμk(j)ˉε)2dt+γ4N,
    $
    (13)

    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

    $ min(μ,τ)   J(μ,τ,z)s.t.   z solves (9),(11) for the switching sequence (μ,τ).
    $
    (14)

    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

    $zjt(t,x)+Ajzjx(t,x)=fj(zj(t,x)),t[0,T],x[0,Lj],zj(0,x)=zj0(x),x[0,Lj],
    $
    (15)

    and the coupling conditions

    $ αkx(v,ek)zk1(t,Lkx(v,ek))=αlx(v,el)zl1(t,Llx(v,el))    ek,elδv,
    $
    (16)
    $ ejδ+vzj2(t,Lj)ejδvzj2(t,0)=qv0(t)
    $
    (17)

    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

    $ A^j = \frac{1}{\varepsilon_j} [01c2j0]
    \ \ \ \ \ \text{ for some } \varepsilon_j, c_j > 0 $

    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

    $ Z = \left[ \underset{j=1}{\overset{n}{\mathop{\otimes }}}\,{{L}^{2}}([0,{{L}_{j}}],{{\mathbb{R}}^{2}}) \right] \otimes L^2([0, \infty), \mathbb{R} ^m), $

    the vectors

    $z=((z1),,(zn),qv1,,qvm),z0=((z10),,(zn0),qv10,,qvm0)
    $

    the operators

    $ A=[A1An]x,     B=1mx     andf(z)=(f1(z1),,fn(zn),0,,0).
    $
    (18)

    Moreover, define the operator

    $ \text{diag}(A, B) = [A00B]
    $

    on the domain

    $ D([A00B])={z=((z1),,(zn),qv1,,qvm)Z|z is absolutely continuous,αkx(v,ek)zk1(Lkx(v,ek))=αlx(v,el)zl1(Llx(v,el))    vV,ek,elδv,ejδ+vzj2(Lj)ejδvzj2(0)=qv(0)    vV}
    $
    (19)

    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

    $ ˙z(t)+[A00B]z(t)=f(z),   t[0,T],z(0)=z0.
    $
    (20)

    Now we can state the following result:

    Theorem 3.1. The operator

    $ \left(D\left([A00B]
    \right), [A00B]
    \right) $

    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

    $ ˙z(t)+[A00B]z(t)=0,    t[0,T],z(0)=z0.
    $
    (21)

    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

    $ z_0 = (\varrho^1_0, q^1_0, \ldots, \varrho^n_0, q^n_0, q^{v_1}_0, \ldots, q^{v_m}_0)^\top \in D(\text{diag}(A, B)) $

    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

    $ q^{v_k}(t, s) = q^{v_k}_0(s+t). $

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

    $ αv=(ejδ+vαj1cj+ejδvαj0cj)>0
    $

    and

    $ ϱv(t)=α1v[ejδ+v(cjϱj0(Ljcjεjt)+qj0(Ljcjεjt))+ejδv(cjϱj0(cjεjt)qj0(cjεjt))+qv(t,0)].
    $

    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

    $ ϱj(t,0)=αj0ϱej(1)(t),qj(t,0)=αj0cjϱej(1)(t)+cjϱj0(cjεjt)+qj0(cjεjt),ϱj(t,Lj)=αj1ϱej(2)(t),qj(t,Lj)=αj1cjϱej(2)(t)cjϱj0(Ljcjεjt)+qj0(Ljcjεjt).
    $

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

    $ \lim\limits_{t \searrow 0} \varrho^j(t, x) = \varrho^j_0(x)\ \ \ \ \ \ \text{and} \ \ \ \ \ \ \lim\limits_{t \searrow 0} q^j(t, x) = q^j_0(x) \ \ \ \ \ \ \text{ for } x \in \{0, L_j\}. $

    We skip these calculations here for brevity. Next, set

    $ [ϱjqj]
    (t, x) = \frac{1}{2}[1εjcjcjεj1]
    [ϱjqj]
    \left(\frac{c_j}{\varepsilon_j} t, x\right) +\frac{1}{2}[1εjcjcjεj1]
    [ϱjqj]
    \left(-\frac{c_j}{\varepsilon_j} t, x \right) $

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

    $ [ϱjqj]
    (t, x) = {[ϱj0qj0](xt), if xt[0,Lj],[ϱjqj](εj(tx)cj,0), if xt<0,[ϱjqj](εj(t(Ljx))cj,Lj), if xt>Lj.
    $

    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

    $ejδ+vqj(t,Lj)ejδvqj(t,0)=ejδ+v[αj1cjϱej(2)(t)cjϱj0(Ljcjεjt)+qj0(Ljcjεjt)]=ejδv[αj0cjϱej(1)(t)+cjϱj0(cjεjt)+qj0(cjεjt)]=(ejδ+vαj1cj+ejδvαj0cj)ϱv(t)=ejδ+v[cjϱj0(Ljcjεjt)qj0(Ljcjεjt)]ejδv[cjϱj0(cjεjt)+qj0(cjεjt)]=αvϱv(t)(αvϱv(t)qv(t))=qv(t)
    $

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

    $Ajx[ϱjqj](t,x)=12Aj[1εjcjcjεj1]x[ϱjqj](cjεjt,x)+12Aj[1εjcjcjεj1]x[ϱjqj](cjεjt,x)=εj2cj[cjεj1cj2ε2jcjεj]t[ϱjqj](cjεjt,x)+εj2cj[cjεj1cj2ε2jcjεj]t[ϱjqj](cjεjt,x)=t[ϱjqj](t,x)
    $

    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

    $ E(t) = \frac{1}{2}\sum\limits_{e_j \in E} \int_0^{L_j} (\varrho^j(t, x))^2 + \frac{1}{c_j^2}(q^j(t, x))^2 dx. $

    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,

    $ddtE(t)=ejELj0ϱj(t,x)tϱj(t,x)+1c2jqj(t,x)tqj(t,x)dx=ejE1εjLj0ϱj(t,x)xqj(t,x)+qj(t,x)xϱj(t,x)dxP.I.=ejE1εj[qj(t,Lj)ϱj(t,Lj)qj(t,0)ϱj(t,0)]=vV[ejδ+v1εjqj(t,Lj)ϱj(t,Lj)ejδv1εjqj(t,0)ϱj(t,0)]=vVϱv(t)[ejδ+vαj1εjqj(t,Lj)ejδvαj0εjqj(t,0)]=0,
    $

    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

    $ \tilde{f}(z) = [0θmin{z2|z2|,ˉq}max{z1,ˉϱ}ghz1]
    . $

    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).

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

    $ minμ,τJ(μ,τ,z)s.t.˙z(t)=Aμkz(t)+fμk(t,z(t)),k{1,,N},t(τk1,τk),z(τk)=gμk,μk+1(z(τk)),k{1,,N},z(τ0)=z0,τT(0,T).
    $
    (22)

    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

    $ \mathcal{T}(0, T) = \{ \tau = (\tau_1, \ldots, \tau_N) \in \mathbb{R} ^N \mid 0 = \tau_0 \leq \tau_1 \leq \ldots \leq \tau_N \leq \tau_{N+1} = T \}. $

    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

    $ [pt(t)γt(t)]+[Aμk00B][p(t)γ(t)]=[[fμkz(z(t))]p(t)+lz(z(t))0],t(τk1,τk),[p(T)γ(T)]=0,k=1,,N+1.
    $
    (23)

    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

    $ [pj1pj2]t+1εμk(j)[0c2j10][pj1pj2]x=(fμk(j)z(zj))[pj1pj2]+γ1[ϱjϱjdqjqjd][pj1pj2](T,x)=0,
    $
    (24)

    where

    $ f^0_z(\varrho, q) = [θc2qj|qj|ϱj|ϱj|gh2θc2|qj|ϱj00]
    \ \ \ \ \ \text{and} \ \ \ \ \ f^1_z(\varrho, q) = [00θqj|qj|ϱj|ϱj|gh2θ|qj|ϱj]
    $

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

    $ αjx(v,ek)pj1(t,Lkx(v,ek))=αkx(v,el)pk1(t,Llx(v,el)),      ek,elδv,
    $
    (25)
    $ ejδ+vpj2(t,Lj)ejδvpj2(t,0)=0.
    $
    (26)

    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

    $ Φτk=ejE[Lj0pj(τk,x)[(Aμk)j(Aμk1)j]zj(τk,x)dx+γ3[(εμk(j)ˉε)2(εμk1(j)ˉε)2]].
    $
    (27)

    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

    $ nj=a(τ,n)Φτj(τ)0     and     b(τ,n)j=nΦτj(τ)0.
    $
    (28)

    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

    $ Φμ(τ)=ejE[Lj0pj(τ,x)[(Aμk)j(Aμ)j]zj(τ,x)dx+γ3[(εμk(j)ˉε)2(εμˉε)2]],
    $
    (29)

    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].

    Algorithm 1 Two-phase gradient descent for stationary switching sequences
    Require: Initial switching sequence $(\mu^0, \tau^0)$ with $N$ modes, Armijo-parameters $\beta \in (0, 1)$ and $\gamma \in (0, 1)$
    1: Set $k=0$, solve (20) for $z$ and (24) for $p$.
    2: Calculate the switching time gradient $\frac{\partial \Phi}{\partial \tau^k}=\left(\frac{\partial \Phi}{\partial \tau^k_n}\right)_{n=1, \ldots, N-1}$ in (27).
    3: while $\tau^k$ does not satisfy (28) do
    4:     Find a step size $s_k = \max \{ \beta^l \, | \, l=0, 1, 2, \ldots \}$ such that
         $ \Phi\left(P\left(\tau^k-s^k\frac{\partial \Phi}{\partial \tau^k}\right)\right) \leq \Phi(\tau^k) -\gamma \frac{\partial\Phi}{\partial \tau^k}^\top \left[\tau^k-P\left(\tau^k-s^k\frac{\partial \Phi}{\partial \tau^k}\right)\right] $
    5:     Set $\tau^k \gets P\left(\tau^k-s^k\frac{\partial \Phi}{\partial \tau^k}\right)$.
    6:     Solve (9) for $z$ and (24) for $p$.
    7:     Calculate the switching time gradient $\frac{\partial \Phi}{\partial \tau^k}=\left(\frac{\partial \Phi}{\partial \tau^k_n}\right)_{n=1, \ldots, N-1}$ in (27).
    8: end while
    9: if the mode insertion gradient $\frac{\partial \Phi}{\partial \mu'}(\tau') \geq 0$ in (29) $\ \ \forall \, \mu' \in \{0, 1\}, \, \tau' \in [0, T]$ then
    10:    return $(\mu, \tau)=(\mu^k, \tau^k)$
    11: end if
    12: Find mode $\mu'$ and $\tau' \in [0, T]$ with $\frac{\partial \Phi}{\partial \mu'}(\tau') < 0$ in (29).
    13: Find $n \in \{0, \ldots, N\}$ such that $\tau' \in [\tilde{\tau}^k_n, \tilde{\tau}^k_{n+1})$.
    14: Set $\mu^{k+1} \gets (\mu^k_1, \ldots, \mu^k_n, \mu', \mu^k_{n+1}, \ldots, \mu^k_N)$
    15: Set $\tau^{k+1} \gets (\tau^k_1, \ldots, \tau^k_n, \tau', \tau', \tau^k_{n+1}, \ldots, \tau^k_N)$
    16: Set $k \gets k+1$, $N \gets N+1$ and go to 2.

     | Show Table
    DownLoad: CSV

    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.

    Figure 1.  A gas network with a supply node $N_1$ and two costumer nodes $N_2$ and $N_3$.

    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.

    Figure 2.  Snapshot of the fully simulated solution showing density (solid, blue) and flux (dashed, red, scaled by $0.05$). On the outer pipes 1 to 5 we see a lot of fluctuation due to the oscillatory boundary flows. The pipes 6 to 9 of the inner circle, however, remain nearly constant.
    Figure 3.  (A): resulting optimized switching sequence showing, for each time step from $t_0 = 0\ \text{s}$ to $T = 1800\ \text{s}$ and each edge $e_1, \ldots, e_{10}$, if the solution is calculated with the fine model (white) or frozen (black). (B), (C): filtered results with two different filters. (D): $L^2$-error relative to maximum values of the solution $\bar{z}$ corresponding to freezing edges $6$ to $10$ completely.

    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

    $zk+12j+12=12(zkj+1+zkj)Δt2ΔxA(zkj+1zkj),zk+1j=zkjΔtΔxA(zk+12j+12zk+12j12)+ε(zkj+12zkj+zkj1)
    $

    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

    $ zv=2|δv|(eδ+v(ze)kn1+[1001]eδv(ze)k2)
    $

    and then using the ghost cell values

    $ (ze)kn+1=[1001](zv(ze)kn1)         for all eδ+v,(ze)k0=zv[1001](ze)k2         for all eδv.
    $

    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

    $zv=21+ˉα((ze+)kn1+[1001](ze)k2),(ze+)kn+1=[1001](zv(ze+)kn1),(ze)k0=[ˉα001]zv[1001](ze)k2.
    $

    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.

    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.

    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.

    [1] M. A. Adewumi and J. Zhou, Simulation of Transient Flow in Natural Gas Pipelines, 27th Annual Meeting of PSIG (Pipeline Simulation Interest Group), Albuquerque, NM, 1995, URL https://www.onepetro.org/conference-paper/PSIG-9508.
    [2] Gradient descent approach to optimal mode scheduling in hybrid dynamical systems. Journal of Optimization Theory and Applications (2008) 136: 167-186.
    [3] Coupling conditions for gas networks governed by the isothermal Euler equations. Networks and Heterogeneous Media (2006) 1: 295-314.
    [4] Gas flow in pipeline networks. Networks and Heterogeneous Media (2006) 1: 41-56.
    [5] MPEC problem formulations and solution strategies with chemical engineering applications. Computers & Chemical Engineering (2008) 32: 2903-2913.
    [6] Flows in networks with delay in the vertices. Mathematische Nachrichten (2012) 285: 1603-1615.
    [7] L. T. Biegler, Nonlinear Programming: Concepts, Algorithms, and Applications to Chemical Processes, vol. 10 of MOS-SIAM Series on Optimization, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2010. doi: 10.1137/1.9780898719383
    [8] A. Bressan, Hyperbolic Systems of Conservation Laws: The One-dimensional Cauchy Problem, vol. 20, Oxford University Press on Demand, 2000.
    [9] Gas pipeline models revisited: model hierarchies, nonisothermal models, and simulations of networks. SIAM Journal on Multiscale Modeling and Simulation (2011) 9: 601-623.
    [10] J. C. Butcher, Numerical Methods for Ordinary Differential Equations, John Wiley & Sons, 2016. doi: 10.1002/9781119121534
    [11] Optimal control of a class of hybrid systems. IEEE Transactions on Automatic Control (2001) 46: 398-415.
    [12] G. Cerbe, Grundlagen der Gastechnik, Hanser, 2016. doi: 10.3139/9783446449664
    [13] Cascading of fluctuations in interdependent energy infrastructures: Gas-grid coupling. Applied Energy (2015) 160: 541-551.
    [14] P. J. Davis and P. Rabinowitz, Methods of Numerical Integration, Courier Corporation, 2007.
    [15] Classical solutions and feedback stabilization for the gas flow in a sequence of pipes. Networks and Heterogeneous Media (2010) 5: 691-709.
    [16] Adaptive refinement strategies for the simulation of gas flow in networks using a model hierarchy. Electronic Transactions Numerical Analysis (2018) 48: 97-113.
    [17] P. Domschke, B. Hiller, J. Lang and C. Tischendorf, Modellierung von Gasnetzwerken: Eine Übersicht, Technical report, Technische Universität Darmstadt, 2017, URL https://opus4.kobv.de/opus4-trr154/frontdoor/index/index/docId/191.
    [18] Adjoint-based error control for the simulation and optimization of gas and water supply networks. Journal of Applied Mathematics and Computing (2015) 259: 1003-1018.
    [19] Transition-time optimization for switched-mode dynamical systems. IEEE Transactions on Automatic Control (2006) 51: 110-115.
    [20] Maximal controllability for boundary control problems. Applied Mathematics & Optimization (2010) 62: 205-227.
    [21] Vertex control of flows in networks. Networks and Heterogeneous Media (2008) 3: 709-722.
    [22] Stationary states in gas networks. Networks and Heterogeneous Media (2015) 10: 295-320.
    [23] M. Hahn, S. Leyffer and V. M. Zavala, Mixed-Integer PDE-Constrained Optimal Control of Gas Networks, Mathematics and Computer Science, URL https://www.mcs.anl.gov/papers/P7095-0817.pdf.
    [24] F. M. Hante, G. Leugering, A. Martin, L. Schewe and M. Schmidt, Challenges in Optimal Control Problems for Gas and Fluid Flow in Networks of Pipes and Canals: From Modeling to Industrial Applications, in Industrial Mathematics and Complex Systems: Emerging Mathematical Models, Methods and Algorithms (eds. P. Manchanda, R. Lozi and A. H. Siddiqi), Springer Singapore, Singapore, 2017, 77-122. doi: 10.1007/978-981-10-3758-0_5
    [25] Modeling and simulation of a gas distribution pipeline network. Applied Mathematical Modelling (2009) 33: 1584-1600.
    [26] Adjoint calculus for optimization of gas networks. Networks and Heterogeneous Media (2007) 2: 733-750.
    [27] Optimal switching between autonomous subsystems. Journal of the Franklin Institute (2014) 351: 2675-2690.
    [28] Second-order switching time optimization for nonlinear time-varying dynamic systems. IEEE Transactions on Automatic Control (2011) 56: 1953-1957.
    [29] Transient analysis of isothermal gas flow in pipeline networks. Chemical Engineering Journal (2000) 76: 169-177.
    [30] Spectral properties and asymptotic periodicity of flows in networks. Mathematische Zeitschrift (2005) 249: 139-162.
    [31] Mild solution and constrained local controllability of semilinear boundary control systems. Journal of Dynamical and Control Systems (2017) 23: 735-751.
    [32] C. B. Laney, Computational Gasdynamics, Cambridge University Press, 1998. doi: 10.1017/CBO9780511605604
    [33] Control parametrization enhancing technique for optimal discrete-valued control problems. Automatica (1999) 35: 1401-1407.
    [34] R. J. Le, Veque, Numerical Methods for Conservation Laws, Birkhäuser, 1992. doi: 10.1007/978-3-0348-8629-1
    [35] R. J. Le, Veque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, 2002. doi: 10.1017/CBO9780511791253
    [36] A mixed integer approach for time-dependent gas network optimization. Optimization Methods and Software (2010) 25: 625-644.
    [37] V. Mehrmann, M. Schmidt and J. Stolwijk, Model and Discretization Error Adaptivity within Stationary Gas Transport Optimization, to appear, Vietnam Journal of Mathematics, URL https://arXiv.org/abs/1712.02745, Preprint 11-2017, Institute of Mathematics, TU Berlin, 2017.
    [38] Hybrid systems of differential-algebraic equations - Analysis and numerical solution. Journal of Process Control (2009) 19: 1218-1228.
    [39] E. S. Menon, Gas pipeline Hydraulics, CRC Press, 2005.
    [40] Pipe networks: Coupling constants in a junction for the isentropic Euler equations. Energy Procedia (2015) 64: 140-149.
    [41] D. Mugnolo, Semigroup Methods for Evolution Equations on Networks, Springer, 2014. doi: 10.1007/978-3-319-04621-1
    [42] Simulation of transient gas flows in networks. International Journal for Numerical Methods in Fluids (1984) 4: 13-24.
    [43] A. Pazy, Semigroups of Linear Operators and Applications to Partial Differential Equations, vol. 44, Springer-Verlag, New York, 1983. doi: 10.1007/978-1-4612-5561-1
    [44] Validation of nominations in gas network optimization: Models, methods, and solutions. Optimization Methods and Software (2015) 30: 15-53.
    [45] Optimal switching for hybrid semilinear evolutions. Nonlinear Analysis and Hybrid Systems (2016) 22: 215-227.
    [46] Optimality Conditions for Switching Operator Differential Equations. Proceedings in Applied Mathematics and Mechanics (2017) 17: 777-778.
    [47] S. Sager, Reformulations and Algorithms for the Optimization of Switching Decisions in Nonlinear Optimal Control, Journal of Process Control, 19 (2009), 1238-1247, URL https://mathopt.de/PUBLICATIONS/Sager2009b.pdf.
    [48] E. Sikolya, Semigroups for Flows in Networks, PhD thesis, Eberhard-Karls-Universität Tübingen, 2004.
    [49] Flows in networks with dynamic ramification nodes. Journal of Evolution Equations (2005) 5: 441-463.
    [50] J. Smoller, Shock Waves and Reaction-Diffusion Equations, vol. 258 of Grundlehren der mathematischen Wissenschaften, Springer, 1983. doi: 10.1007/978-1-4612-0873-0
    [51] Switched-mode systems: Gradient-descent algorithms with Armijo step sizes. Discrete Event Dynamic Systems: Theory and Applications (2015) 25: 571-599.
    [52] Optimal control of switched autonomous systems. Proceedings of the 41st IEEE Conference on Decision and Control (2002) 4: 4401-4406.
    [53] Optimal control of switched systems based on parameterization of the switching instants. IEEE Transactions on Automatic Control (2004) 49: 2-16.
    [54] Optimal control of hybrid switched systems: A brief survey. Discrete Event Dynamic Systems: Theory and Applications (2015) 25: 345-364.
  • This article has been cited by:

    1. Simone Göttlich, Falk M. Hante, Andreas Potschka, Lars Schewe, Penalty alternating direction methods for mixed-integer optimal control with combinatorial constraints, 2021, 188, 0025-5610, 599, 10.1007/s10107-021-01656-9
    2. Sara Grundel, Michael Herty, Hyperbolic discretization of simplified Euler equation via Riemann invariants, 2022, 106, 0307904X, 60, 10.1016/j.apm.2022.01.006
    3. Zlatinka Dimitrova, Flows of Substances in Networks and Network Channels: Selected Results and Applications, 2022, 24, 1099-4300, 1485, 10.3390/e24101485
    4. Falk M. Hante, 2020, Chapter 7, 978-981-15-0927-8, 157, 10.1007/978-981-15-0928-5_7
    5. Andreas Rupp, Markus Gahn, Guido Kanschat, Partial differential equations on hypergraphs and networks of surfaces: Derivation and hybrid discretizations, 2022, 56, 2822-7840, 505, 10.1051/m2an/2022011
  • Reader Comments
  • © 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5247) PDF downloads(328) Cited by(5)

Figures and Tables

Figures(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog