OPTIMAL MODEL SWITCHING FOR GAS FLOW IN PIPE NETWORKS

. We consider model adaptivity for gas ﬂow in pipeline networks. For each instant in time and for each pipe in the network a model for the gas ﬂow 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 diﬀerential 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.


(Communicated by Axel Klar)
Abstract. 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. [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 mixedinteger 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.1.
A model hierarchy for single pipes. The motion of a compressible nonviscous 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 denotes the density in kg m 3 , v the velocity in m s , P the pressure of the gas in kg m s 2 , g the gravitational constant and h the slope of the pipe. Furthermore θ = λ 2D , where λ 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 and the flux q = 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 s is the specific gas constant. For an ideal gas one has Z(P ) ≡ 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 = P/ obtained under the assumption of a constant gas compressibility factor Z(P, T 0 ) ≡Z. This model neglects inertia effects and assumes small velocities |v| c. For natural gas, indeed one typically has |v| ≤ 10 m s and c ≈ 340 m s , see e.g., [17]. This PDE exhibits two simple characteristics with speeds λ 1 = −c and λ 2 = c.
For stationary situations, we consider the model ∂ x q = 0, obtained from (2) with ∂ t = ∂ t q = 0. Here, the flux q is constant in space and time q(t, x) =q and (t, x) =¯ (x) with¯ being a solution of the momentum equation in (3), which is a Bernoulli-equation. A solution of the momentum equation can therefore be obtained algebraicallȳ 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 ∈ N we consider a network of pipes that we model by a metric graph G = (V, E) with nodes V = (v 1 , . . . , v m ) and edges E = (e 1 , . . . , e n ) ⊆ V × V . For each edge e ∈ 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) = e(2) for any e ∈ E and thus self-loops are not allowed. On the other hand, if v ∈ V is any node, then we define the set of ingoing edges by δ The number |δv| then is called the degree of node v ∈ V .
With each edge e j ∈ 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 ∈ V and e j ∈ δv For each node v ∈ 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 j weighted by given factors α ∈ (0, ∞) m×2 coincide for all incident edges e ∈ δv and can be expressed as The nodal balance equation for a given outflow function q v : [0, T ] → R is similar to a classical Kirchhoff condition for the fluxes q j and can be written as The choice of α corresponds to the nodal types in the network. Prototypically, we will consider the following cases: Junctions: nodes v such that q v ≡ 0 and α k x(v,e k ) = 1 for all e k ∈ δv. Boundary nodes: nodes v such that α k x(v,e k ) = 1 for all e k ∈ δv, but q v ≡ 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 ≡ 0 and |δ + v| = |δ − v| = 1. A description established via the characteristic diagram based on measured specific changes in adiabatic enthalpy H ad of the compression process yields the model where κ is a compressor specific constant,Z is the gas compressibility factor that is assumed to be constant and H ad is within flow dependent and compressor specific bounds obtained from the characteristic diagram. In consistency with the pipe models, we assume that H ad is given by a known referenceH ad . Then we get with a compressor specific factor This yields α k 1 = 1 and α l 0 =ᾱ. The compression ratio of centrifugal and turbo compressors mostly lies within the rangeᾱ ∈ [1, 2], withᾱ ≈ 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 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 ∈ [0, T ] and x ∈ [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 : R 2 → R 2 is globally Lipschitz-continuous. The same transformation applied to (3), together with given boundary values To define an appropriate hyperbolic system to approximate (6), for ε > 0 and any initial values w 0 ∈ L 1 ([0, L], 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: Proof. The global existence and uniqueness ofw follows from standard ODE-theory. For (t, x) ∈ Ω, the broad solution w can be written as Its existence and the continuity on Ω is shown in [8, Theorem 3.1, Theorem 3.4]. Given a fixed t 0 ∈ (0, T ), we obviously can choose ε small enough such that [t 0 , T ] × [0, L] ⊆ Ω. We then have for ε 0 and all (t, 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 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 ∈ V is given by the parameters α and q v as in Section 2.2. On each edge e j of length L j we have an initial gas density j 0 and gas flow q j 0 . For an increasing sequence of switching for a fixedε > 0 withε 1 and denote the right-hand sides in (8) and (2), respectively. Let L j ] be given initial states. Moreover, at any time point t ∈ [τ 0 , τ N +1 ], density and flux additionally satisfy the node coupling conditions for each v ∈ V , explained in Section 2.2. Denote by z d = ( 1 d , q 1 d , . . . , n d , q n d ) the reference solution to (9), (11) for the choice N = 1, µ = 1 and τ = (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 = ( 1 , q 1 , . . . , n , q n ) then define the cost functional with γ 1 , . . . , γ 4 ≥ 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 The challenge now is to choose the sequences µ and τ 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 s.t. z solves (9), (11) for the switching sequence (µ, τ ).
3. Semigroup formulation on networks. In this section, we will set up an abstract formulation of the PDE-system in (9) for τ = (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 ∈ {1, . . . , n}, we now consider the initial boundary value problem for z j = ( j , q j ) on pipe e j ∈ E given by and the coupling conditions for each node v ∈ V . Here, the operators

Lipschitz-continuous and
Moreover, define the operator 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 probleṁ Now we can state the following result: is the infinitesimal generator of a strongly continuous semigroup on Z.
Proof. Note that diag (A, B) is a densely defined, closed operator on Z with a nonempty resolvent set (for instance 0 ∈ ρ(diag(A, B))). Following [43, Chapter 4, Theorem 1.3], to prove the claimed semigroup property, it thus suffices to show that the homogeneous systeṁ Since 1 , q 1 , . . . , n , q n , q v1 , . . . , q vm are absolutely continuous, so is v for each v ∈ V . For each edge e j ∈ E, j = 1, . . . , n, construct the functions j , q j : [0, T ]×[0, 1] → R as follows: for t ∈ (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 The above considerations show that ( 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 ∈ V we have Observe that j and q j are continuous at the boundaries x = 0 and x = L j and that Obviously, we have E(t) = 0 if and only if j (t, .) = q j (t, .) = 0 for all j = 1, . . . , n, E(t) ≥ 0 for all t ≥ 0 and E(0) = 0. Furthermore, But then E ≡ 0, hence z ≡ 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¯ ,q > 0 and define the modified right-hand sidẽ The functionf is continuously differentiable and globally Lipschitz-continuous as a function mapping L 2 ([0, L], R) onto itself for any L > 0. If we replace f by thẽ f in (9), we then have a system fitting the assumptions made above for each fixed k = 1, . . . , N and thus we have a unique solution. If we choose¯ sufficiently small andq sufficiently large, then the numerical simulation with realistic data shows that both the boundaries¯ andq will in fact never be reached by 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 min µ,τ J(µ, τ, z) where A µ k is a strongly continuous semigroup, f µ k is a semilinear perturbation and g µ k ,µ k+1 = id is a transition map that can be chosen trivially here for each k = 1, . . . , N . The ordering cone 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 (µ, τ ) → z(µ, τ ) that can be substituted into J to define the reduced cost functional Φ(µ, τ ) = J(µ, τ, z(µ, τ )). We can now apply [45,Theorem 8] to show that Φ is continuously differentiable with respect to the switching time τ k . Recalling the shortened notation introduced in (13), we find that, while J 2 can be differentiated directly with respect to (µ, τ ), the term J 1 depends on the solution z = z(µ, τ ). Again by [45,Theorem 8], we can state gradient formulae for the derivatives of J 1 based on the adjoint equation Here, f µ z : Z → L(Z) and l z : Z → Z * denote the derivatives of the friction term f µ and l as in (13) with respect to z. Substituting the definitions made in (18) yields that in fact γ ≡ 0 and p = (p j ) j=1,...,n again can be partitioned into edgewise defined functions p j = (p j 1 , p j 2 ) for j = 1, . . . , n satisfying p j where are the derivatives of f 0 and f 1 in (10), as well as the coupling conditions α j x(v,e k ) p j 1 (t, L k x(v, e k )) = α k x(v,e l ) p k 1 (t, L l x(v, e l )), e k , e l ∈ δv, (25) 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 µ is fixed, then we can conclude that a switching sequence τ = (τ k ) k=0,...,N +1 is a KKT-point of the minimization problem (14), if the following holds: For n ∈ {1, . . . , N } set a(τ, n) = min{m ∈ {0, . . . , n} | τ m = τ n } and b(τ, n) = max{m ∈ {n, . . . , N + 1} | τ m = τ n }, then n j=a(τ,n) Similarly, by [45,Theorem 10], the sensitivity of the cost function with respect to introducing a new mode µ on an infinitesimal time interval at the point τ can be represented by the mode insertion gradient given by where µ k is the original mode at time τ . In summary, a switching sequence (µ, τ ) is called stationary, if τ is a KKT-point for µ fixed and if ∂Φ ∂µ (τ ) ≥ 0 for all modes µ and all times τ ∈ [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 ODEdynamical systems [2]. The main idea is a two phase approach as follows: the algorithm is initialized with a switching sequence (µ 0 , τ 0 ), for instance µ 0 = 1 and τ 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 T (0, T ) is used. In a second phase, a new mode µ is inserted at a time point τ where ∂Φ ∂µ (τ ) < 0. If no such point exists, the switching sequence is stationary in the above sense, otherwise the algorithm continues with the first phase again.
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ᾱ, compare (5), and at initial time the network is assumed to be stationary with constant densities 0 on the outer circle (thus α 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 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, lettingz 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 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 , . . . , z k n ) with spatial step size ∆x of the solution at time point t k , it computes the discretized values z k+1 at time t k+1 = t k + ∆t by the explicit finite-volume-scheme for j = 1, . . . , n, where ε ∼ 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  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 timedependent 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.