MODELING, SIMULATION AND OPTIMIZATION OF GAS NETWORKS WITH COMPRESSORS

We consider gas flow in pipeline networks governed by the isothermal Euler equations and introduce a new modeling of compressors in gas networks. Compressor units are modeled as pipe–to–pipe intersections with additional algebraic coupling conditions for the compressor behavior. We prove existence and uniqueness of solutions with respect to these conditions and use the results for numerical simulation and optimization of gas networks.


Introduction.
Recently, there has been intense research in physical phenomena posed on networks with applications for example in traffic, gas and water flow [6,17,15,20,21,2,13,12]. This publication is concerned with the modeling, simulation and optimization of gas networks. In particular the optimization of gas transport in networks is an important industrial problem and has been under investigation for several years, see [30,4,16,37,34,31] and the publications of the Pipeline Simulation Interest Group [29]. Models for transient flow of different mathematical complexity have been proposed and numerical methods for simulation and optimization have been applied, see [16,27]. Our starting point for describing transient gas flow is the isothermal Euler equations as e.g. in [2,3,31,16,9]. Using the properties of solutions to Riemann problems for those equations, coupling conditions for intersections of multiple pipes have been derived and simulated [3,28,10,16].
The contribution of this publication will be in the modeling, simulation and optimization of additional compressors units: in realistic gas networks we usually face a pressure drop along the pipes due to friction effects. To satisfy the demand of consumers, it is necessary to compensate these losses, which is usually done by compressors. The main objective is to drive the gas through the network and to run compressors cost efficiently. Several approaches to this problem are already known: most approaches deal with simplifications of the isothermal Euler equations before modeling and optimizing the compressor run-times, e.g., see [27,40,14,38] for linear and mixed-integer models based on steady-state optimization or linearization approaches and [16,34] for a simplified, nonlinear approach. On the contrary, we propose a modeling of compressors (later its optimization) based on the full isothermal Euler equations without linearization or steady-state assumptions. The idea is to discuss a compressor as a special kind of pipe-to-pipe intersections. This allows us to use the notion of 1 2 -Riemann problems [21] to define solutions for a network problem including compressors stations. Additional conditions have to be imposed to obtain unique solutions, as also seen in other network problems [17,21,3]. For gas dynamics we obtain these conditions by applying a similar reasoning as in the engineering community [28,10,37,31,29].
The outline of the paper is as follows: we recall the isothermal Euler equations and its known coupling conditions for gas networks in Section 2. Next, we present the modeling of compressors and show existence of solutions to a simple network configuration including a controlled compressor. Finally, we present numerical results on simulation and optimization of a sample network.

The isothermal Euler equations as model for transient network flow.
The isothermal Euler equations are derived as simplification of the isentropic equations which in turn are derived from the Euler equations under the assumption that the energy equation is redundant. In the isothermal Euler equations the temperature is constant and we have the pressure law where p is the pressure, z is the natural gas compressibility factor, R the universal gas constant, T the absolute gas temperature, M g is the gas molecular weight and ρ the density. To keep the presentation simple but still valid for real-world gas networks [16,29] we additionally impose the following: there is negligible wall expansion or contraction under pressure loads; i.e. pipes have constant cross-sectional area. All pipes have the same diameter D and the constant a 2 = ZRT /M g , i.e., the sound speed, is the same for all pipes in the network. In all pipes we assume steady state friction [27,39]. The friction factor f g is calculated using Chen's equation [5]: where N Re is the Reynolds number N Re = ρuD/µ, µ the gas dynamic viscosity and ǫ the pipeline roughness, which are again assumed to be the same for all pipes. Finally, the equations which govern the flow inside a pipe j of the network read Here, u j denotes the velocity of a gas and ρ j u j is the flux in pipe j. The first equation is the conservation of mass and the second states the conservation of moment. We use the following notation Before we turn to the modeling of compressors we recall the major results of the discussion of pipe-to-pipe intersections presented in [2,3], see also the approaches in the engineering community [28,10,16].

OPTIMIZATION OF GAS NETWORKS 83
As in [2] we introduce three assumptions (A1)-(A3) on possible solutions to a network problem to be defined below.
A3 The pressure at a pipe-to-pipe intersection is a constant.
Remark 1. The condition (A3) is frequently used in the engineering literature [28] as well as in the simplified models mentioned before [16,27,38]. It relies mostly on the assumption that the velocity of the gas under consideration is small compared with the applied pressure which is the case for most commercial gas pipelines, see [31]. However, a realistic modeling of a pipe intersection would also have to include further physical effects. Among them, the geometry of the fitting is believed to have the strongest effects on the physical behavior near the fitting, see [8] for a possible mathematical model including geometry effects. Additionally, the material properties as well as possibly cross-section areas of connected pipes influence the coupling conditions. Such effects have been taken into account in the engineering literature by introducing 'minor losses' at fittings. These losses are modeled as pressure losses at the pipe intersection and have been listed in tables depending on the geometry, type of gas and material properties in [10] for example. A more mathematical discussion of different coupling conditions can be found in [2] and in [9,7]. In the latter references the authors replace (A3) by the following assumption: the dynamic pressure is equal among all connected pipes, i.e., ρ j u 2 j + p(ρ j ) = const at pipe fittings. We refer to [9] for more details and a mathematical comparison with the assumption (A3).
Finally, a more detailed picture of the dynamics including the geometry and properties of the pipe walls can be obtained by introducing a 2d model of the pipe intersection.
Due to assumption (6) we model a network of pipes as a directed, finite graph (J , V) where in addition we may connect edges tending to infinity. Each edge j ∈ J corresponds to a pipe. Each pipe j is modeled by an interval [x a j , x b j ]. For edges ingoing or outgoing to the network we have x a j = −∞ or x b j = +∞, respectively. Each vertex v ∈ V corresponds to an intersection of pipes. For a fixed vertex v ∈ V we denote by δ − v (δ + v ) the set of all indices of edges j ∈ J ingoing (outgoing) to the vertex v. On each edge j ∈ J we assume that the dynamics is governed by the isothermal Euler equations (3) for all x ∈ [x a j , x b j ] and t ∈ [0, T ] and supplemented with initial data U 0 j . We give a brief summary of the coupling conditions and existence results at a vertex v ∈ V which are already known.
First, at pipe-to-pipe intersections we couple systems of the type (3) by boundary conditions. To this end we introduce "intermediate" states at the vertex, similar to [17,21,22,2]. We have one intermediate state for each connecting pipe and those states have to satisfy the coupling conditions at the vertex. To be more precise: Definition 1. Consider a single vertex v with ingoing pipes j ∈ δ − v and outgoing pipes j ∈ δ + v . Assume constant initial data U 0 j given on each pipe. A family of functions (U j ) j∈δ − v ∪δ + v is called a solution at the vertex, provided that U j is a weak entropic solution on the pipe j and U j satisfies (A1)-(A2) and (for U j sufficiently regular) The first condition resembles Kirchoff's law and is referred to as Rankine-Hugoniot condition at the vertex. It imposes the conservation of mass at the intersection. The second condition resembles the assumption (A3), i.e., (7), and has been extensively discussed in [2,28,10]. For results on existence and uniqueness for pipe-to-pipe intersections in the sense of Definition 1 we refer to [2,3], in particular to Section 3 of [2], where the existence in the case of vertices with degree less than four has been proven.
Second, we introduce 1 2 -Riemann problems which are the key to construct solutions at pipe intersections and at compressors later on: a 1 2 -Riemann problem at a fixed junction v ∈ V is a Riemann problem for a single pipe j as if the pipe is extended to (−∞, ∞) with initial data U 0 j ,Ū j ordered depending whether the pipe is incoming or outgoing to the junction, i.e., Last, we recall some mathematical properties of Riemann problems for the isothermal Euler equations. Both characteristic families are genuine nonlinear for ρ j > 0. The eigenvalues are λ 1,2 j (U ) = q/ρ∓a. In a Riemann problem for piecewise constant initial data U 0 j (x) with discontinuity at x = 0 we refer to U 0 j (0−) and U 0 j (0+) as left and right states, respectively. We follow the standard theory [11] and notation: we consider the wave curves in the (ρ, q)−plane. They are connected smoothly iñ U up to the second derivative and sketched in Figure 1. A parametrization of the i-shock and i-rarefaction wave curve [26] are given as follows: a point U can be connected to U l by a 1-(Lax-)shock, (resp. 2-(Lax-)shock) if and only if there exists ξ such that with shock speeds given by A parametrization of the 1-rarefaction (resp. 2-rarefaction) curve is given by We refer to the composition of the i-(Lax-)shock curve and the i-rarefaction wave curve through a given (left) stateŨ shortly as i−wave curve. The 1−wave curve is concave and intersects {(ρ, 0) : ρ > 0} exactly once. The 2−wave curve is convex and also intersects {(ρ, 0) : ρ > 0} exactly once. For a given stateŨ we consider the 1-wave curve through the left stateŨ in the first quadrant of the (ρ, q)−plane. The demand function ρ → d(ρ;Ũ ) is then defined as follows: starting at (ρ, q) = (0, 0) we follow the 1-wave curve as long as it is increasing and until we either reach the stateŨ or the turning point of the 1-wave curve. Once we reached one of these points (ρ * , q * ) we extend the function to ρ > ρ * by the constant q * , c.f. [21,25] in the case of scalar conservation laws.
A state U = (ρ, q) is called subsonic, iff |q/ρ| < a and supersonic iff |q/ρ| > a. The line q = aρ in the (ρ, q)−plane is the sonic line. Using the notion of supply and demand functions we obtain the following results.
Lemma 1 (Proposition 3 [2]). We consider an ingoing pipe j ∈ δ − v and a fixed and given state U 0 j =: U l which satisfies (A1) and (A2), i.e., U 0 j is such that ρ 0 j > 0 and q 0 j ≥ 0. Further, we consider an arbitrary real non-negative number q * with Then, there exists a unique stateŪ j =: U r satisfying (A1), (A2) and q r = q * such that the 1 2 -Riemann problem (9a,9b) admits a solution which is either constant or consists of waves with negative speed only. Lemma 2. We consider an outgoing pipe j ∈ δ + v and a fixed given state U 0 j =: U r satisfying (A1) and (A2). Further, we consider a arbitrary real non-negative number q * .
Then, there exists at least one stateŪ j =: U l satisfying (A1), (A2) and q l = q * and such that the 1 2 -Riemann problem (9a, 9c) admits either the constant solution U l ≡ U r or consists of waves of positive speed.
The assertion of the lemma is a consequence of the following facts: For a fixed state U r with q r ≥ 0 the part of the reversed 2-wave curve through U r which belongs to the first quadrant of the (ρ, q)−plane is monotone increasing. Furthermore, if we connect any state U l with q l ≥ 0 to the given state U r (with q r ≥ 0) by a wave of the second family only, then this wave has non-negative speed.
3. Modeling of compressors in gas networks. In real-world gas networks the pressure along a pipe decreases due to strong friction effects [28,37]. The purpose of a compressor station is to increase the pressure of a gas flow q from p in to p out . We consider centrifugal compressors which require under adiabatic conditions the theoretical power of to fulfill the desired pressure increase. Here, c 1 is a compressor dependent constant, γ is the isentropic coefficient of the gas, T in is the temperature and z = z(p, T ) is again the compressibility factor (see American Gas Association or [37]). Usually, the compressor is consuming fuel directly from the gas flow when operating and hereby reducing the flow q. Moreover, a typical compressor station consists of several compressor units accompanied by bypass valves. For a concise mathematical treatment we simplify this situation as follows: we model the compressor station by a single, idealized compressor with characteristic (14). Further, since the fuel consumption is usually less than 0.5% of the gas throughput, we neglect flow reduction due to operating compressors. Last, we consider a constant compressibility factor z(ρ, T ). These assumptions are as in [34,27]. Finally, we obtain the idealized compressor power necessary to increase the pressure from p in to p out and a flow q as for some constant c 0 > 0. In the context of our network model we treat an idealized compressor as a node in the network and therefore as a particular case of a 1-to-1 pipe intersections. Hence, we split the set of all vertices as V = V N ∪ V C , where V N are pipe-to-pipe intersections and V C are nodes having degree two and representing the compressor stations. We assume Definition (1) holds for all v ∈ V N and define for v ∈ V C : are called a solution at vertex v ∈ V C , provided that U i satisfies (A1)-(A2) and U i is a weak entropic solution on the pipe and (for U i sufficiently regular) holds for all t > 0.
Some remarks are in order.

Remark 2.
In the case P = 0 we recover solutions in the sense of Definition (1) and in particular hypothesis (A3): p(ρ 1 (x b v , t)) = p(ρ 2 (x a v , t)). The definition easily extends to the case of time-dependent power changes by replacing P by a function t → P(t).
If we want to include flux reduction due to fuel consumption we have to modify (16a) as follows: where c 3 describes the percentage of flux used by the compressor.
Next, we prove existence and uniqueness under restrictions on the initial data. By the same reasoning as in [17,2] we only expect solutions in the case of a single pipe intersection.
Recall the following properties of the 1-wave curve in the first quadrant of the (ρ, q)−plane. Given a left state U l with q l ≥ 0. Then, for any U on the 1-rarefaction wave curve through U l , the slope of the 1-rarefaction wave curve at U is equal to λ 1 (U ). Hence, if and only if U l is subsonic, then we have λ 1 (U l ) = q l /ρ l − a < 0. Therefore, for any subsonic state U l , the demand d(ρ l ; U l ) exceeds q l . This is summarized in the following lemma. As a consequence of the previous lemma and Lemma 1 we obtain: given a subsonic left state U l . Then, all subsonic states U r belonging to the 1-wave curve through U l can be connected to U l by a wave of negative speed. Using Proposition 3 [2] and Lemma 2 we conclude: given a subsonic right state U r , then all left states which can be connected to U r by a wave of positive speed belongs to the reversed 2-wave curve through the right state U r . Now, consider a vertex v ∈ V C and subsonic initial data U 0 1 ≡ U 0 2 . Obviously, the constant solution U i (x, t) = U 0 1 for i = 1, 2 is a solution in the sense of definition (2), if P = 0. The next result is a perturbation result of this (uncontrolled) situation and yields existence in the case P ≥ 0 but sufficiently small. Theorem 1. Consider a vertex v ∈ V C and subsonic initial data U 0 1 ≡ U 0 2 satisfying (A1) and (A2) and q 0 1 = 0. Then, for any P ≥ 0 sufficiently small, there exists a solution U i (x, t) i = 1, 2 in the sense of Definition 2 which is subsonic for all t > 0. The solution is unique in this class.
We give an example of the statesŪ 1 andŪ 2 in the phase space (ρ, q) in Figure 2: the initial state is U 0 1 = U 0 2 = (1, 1) and we continuously increase P starting from P = 0. The further parameters are γ = 5/3 (monoatomic gas), c 0 = 1 and a = 5. We observe that a working compressor influences the states before and after the compressor unit. Figure 2 shows the different statesŪ 1 andŪ 2 in the left part and the dependence ofρ 1 ,ρ 2 andq on the compressor power P in the right part of the figure. Finally, we discuss the general case U 0 1 = U 0 2 . Here, we proceed as follows: given subsonic states U 0 1 and U 0 2 which satisfy (A1) and with u 0 1 > 0 and u 0 2 > 0. Let P ≥ 0. Then for U 0 1 − U 0 2 sufficiently small, there exists an intermediate state U m = (ρ m , q m ) as intersection of the 1-wave curve through the left state U 0 1 and the reversed 2-wave curve through the right state U 0 2 , due to general existence theory for Riemann problems of genuine nonlinear conservation laws [11]. If U m is subsonic and satisfies u m > 0, then, we proceed similar as in the proof of Theorem 1: we consider the function where the functions (ρ 1 (ξ), q 1 (ξ)) and (ρ 2 (τ ), q 2 (τ )) the parametrization of the 1wave curve through left state U 0 1 , respectively the reversed 2-wave curve through the right state U 0 2 is. Due to the existence of U m , there is a unique ξ m and a unique τ m , such that The specific form of the functions ρ 1 (ξ) and q 1 (ξ) depend on whether U m is connected to U 0 1 by a 1-rarefaction wave or by a 1-shock wave. Similarly, this holds true for ρ 2 (τ ) and q 2 (τ ). As in the proof of the previous theorem, we compute the determinant of DF at the point (ξ m , τ m ) to apply the inverse function theorem: Hence, detDF (ξ m , τ m ) only vanishes, if the vectors and d dτ ρ 2 (τ m ) d dτ q 2 (τ m ) are linearly dependent. But this is impossible, due to the shape of the wave curves in the ρ − q−plane and since the system (3) is strictly hyperbolic in U m . Therefore, detDF (ξ m , τ m ) = 0. We conclude now exactly as in the previous proof: there exists a neighborhood N (τ m , ξ m ) such that F is bijective and this yields statesŪ 1 andŪ 2 on the 1-wave curve through U 0 1 and reversed 2-wave curve through U 0 2 , respectively, and such that (16) is satisfied. Due Lemma 1 and 2 we obtain: the solutions to the 1 2 -Riemann problems (9a,9b) (resp. (9a,9c) ) consists of waves of non-positive (resp. non-negative) speed and satisfy U 1 (x b 1 , t) =Ū 1 (resp. U 2 (x a 2 , t) =Ū 2 ). Hence, we obtain a solution in the sense of Definition 2. We summarize this findings in the following theorem.
Theorem 2. Given subsonic states U 0 1 and U 0 2 satisfying (A1) and u 0 1 , u 0 2 > 0. Let U 0 1 − U 0 2 be sufficiently small, such that the point of intersection U m of the 1-wave curve through the left state U 0 1 and the reversed 2-wave curve through the right state additionally satisfies 0 < q m /ρ m < a.
Then, for any P ≥ 0 sufficiently small, there exists a solution {U i } i=1,2 in the sense of Definition 2.

Some remarks are in order.
Remark 3. In Lemma 2 and in the case P = 0 we obtainŪ 1 =Ū 2 = U m .
The condition 0 < u m ≡ q m /ρ m is necessary in the case P = 0. Otherwise the proof of Theorem 1 fails, since the inverse function theorem does not apply.
As simplification, the compressor has been treated as a node in the network and therefore no friction effects apply at the compressor station.
Even so the assumption on U i to be subsonic for all times t > 0 restricts the set of possible solutions, it is satisfied for realistic gas network simulations, see [16,37,36]. 4. Optimization problems and numerical results. Typically, gas suppliers have to guarantee specific pressures at particular nodes in the network, e.g., customers: for some given functions φ v (t) and a set of demand nodes V D . Therefore, the task is to minimize the compressor power such that still certain pressures are delivered at demand nodes in the network. Usually, the fuel consumption of compressor stations is used as an objective function for the optimization problem for the compressor power P, see [27,16]. The fuel consumption per unit time for a single idealized compressor can be modeled in a simplified situation by [34] with a positive constant c 2 .
Next, we present numerical results for some sample networks and for the optimization problem (24). As in Remark (3) we use the a diameter of D = 1 [m] for each pipe in the network. For the numerical results we use a higher-order finite volume scheme [23,26,33,18]. The source term is integrated exactly after splitting the system. In all computations we use a space discretization of 400 gridpoints in space and a time-discretization according to the CFL condition. Example 1. We show the influence of a compressor on the dynamic on the pipes. Consider two connected pipes with a compressor located at x = 1 with constant initial data U 0 1 = [1, 1] and U 0 2 = [1, 2] and sound speed a = 5. We neglect friction effects in this example and set c 0 = 1 and γ = 5/3. We simulation solutions in the sense of Definition 2. We present results on the evolution of the densities ρ 1,2 and the fluxes q 1,2 in the case P = 0 in Figure 3 and in the case P = 3/4 in Figure 4. In both figures the area {(t, x) : x < 1} belongs to pipe one and {(t, x) : x > 1} to pipe two. We observe the different shapes of the solution depending on the state of the compressor: the additional discontinuity in the density in Figure 4 due to the working compressor as well the additional kink in the evolution of the flux. As expected, the flux is in all cases continuous at x = 1, due to (8) and (16), respectively. The solution in the uncontrolled case in Figure 3 coincides with the solution to classical Riemann problem for the isothermal Euler equations with initial data U l := U 0 1 and U r := U 0 2 .
Example 2. We treat a simple optimization problem under realistic gas flow conditions. We assume a network of two connected pipes with a compressor located at x 0 = 200 [cm]. The sound speed is a = 377 [m/s]. We present results including friction inside the pipes where the friction factor is given by equation (2). The other parameters are γ = 5/3 and c 0 = 1. We consider a steady state as an initial condition and use the following boundary conditions. We fix a pressure of p(x, t) = 65 [bar]  at the inlet x = 0, an inflow of q(x, t) = 4 · 10 5 [m 3 /h] and the same flow rate at x = 300[cm]. In the uncontrolled case (P = 0), we observe a constant pressure decrease due to friction and a final pressure of about p(x = 300, t) = 35 [bar], see Figure 5. Then, we optimize for a constant compressor power P(t) = P, such that we deliver a pressure ofp = 42 [bar] at T = 1.5[s] at x = 300[cm], i.e., we solve a problem (24) with φ(t) = χ T (t) 65 [bar] and c 2 = c 0 = 1. We present the optimized pressure evolution in Figure 5. We observe the discontinuity in the pressure at x = 200[cm] due to the running compressor unit. Moreover, the compressor effects both, the incoming and the outgoing pressure. The increase at x = 200[cm] needs approximately 0.5[s] to influence the pressure at the boundary. In Figure 6 we present parts of the history of the optimization: we depict the pressure evolution p(x, T ) at the desired time T for several choices of P occurring during the optimization. In this case, a bisection method yields a final compressor power of P = 3.7125 within 12 iterations such that p(x = 300, T ) −p ≤ 10 −3 .   Example 3. In the last example we investigate the network of Figure 7 with two compressors. Initially, we assume a constant state in the network. Then, we consider a situation of a sinodial pressure change on both inlet pipes. If both compressors are turned off, this results in a pressure drop at the outlet (B) as depicted in Figure 8. We reformulate the optimization problem to determine the compressor load:  Figure 9. Contour lines of the pressure in the top row of pipes in the non controlled (top) and optimized situation (bottom). The pipe-to-pipe and pipe-to-compressors nodes are at x = 1/3 and x = 2/3, respectively.
with positive weights α 1,2 . Depending on the objective other approaches should be used to treat the inequality constraints, for example, log-barrier methods or SQP approaches. For our sample problem the previous formulation combined with a direct-search method is sufficient. The further parameters of the model problem are α 1 = 100, α 2 = 1000, a = 5, c 2 = c 0 = 1, γ = 5/3 and φ v (t) = 25. The directsearch method terminated after 30 evaluations of the objective function with an optimal compressor energy of P opt 1 = 0.3 and P opt 2 = 5 · 10 −4 for the compressor in the top and bottom row, respectively. The corresponding pressure evolution at the demand node is also shown in Figure 8 and we observe that the constraint (22) is satisfied. The evolution of the pressure inside the pipes in the top row is depicted in Figure 8 for the two situations P = (0, 0) (uncontrolled) and P = P opt .
The pipe-to-pipe and pipe-to-compressor nodes are located at x = 1/3 and x = 2/3. In the uncontrolled case we see that the sinodial pressure of the top inlet is transported through the pipe system together with the pressure change arriving at time t ≈ 7 from the bottom inlet. The pressure is continuous through x = 2/3 since the compressor is not working, c.f. (16). On the other hand, we observe a discontinuous pressure in the controlled case and a far more complex wave pattern due to the generated and interacting waves arising from both compressor stations in the controlled case. The effects are a superposition of waves generated by the compressors (e.g. x = 2/3, t = 0) and the inflow profiles.
We conclude with a remark on the optimization methods applied here.

Remark 4.
In both examples a gradient-free method has been used to determine the optimal compressor powers. In the first case, the optimization problem has only one (time-independent) unknown and therefore the bisection method can be used to efficiently solve the problem. The method can be found in standard textbooks, e.g. [35]. In the other example we used the nonlinear, gradient free Nelder-Mead simplex method, see [24]. This method extends the simplex algorithm to nonlinear minimization problems. The used algorithms are limited to low-dimensional optimization problems and they are in general very inefficient when applied to higher dimensional problems. If, for example, we are interested a time-dependent compressor control, then a more efficient treatment of the optimization problems using gradient information on the cost functional is necessary. This can be achieved for example by deriving an adjoint and sensitivity calculus for (24). Having gradient information at hand quasi-Newton method can be used to efficiently compute the optimal compressor loads, see again [35] for a variety of such methods. The derivation and the properties of the adjoint equations will be subject to a forthcoming publication.

5.
Conclusion. We presented a model for compressors in gas networks governed by the isothermal Euler equations. We prove existence of solutions to compressor controlled gas networks and implemented the derived conditions in a numerical scheme for optimizing the compressor fuel consumption subject to pressure demands at terminal nodes of the network. In the numerical examples, we showed, how the compressor influences the pressure and flow dynamics in the pipe. Future work will concentrate on the efficient treatment of the derived optimal control problem.