1.1. Modeling of Gas Flow in a Single Pipe
The Euler equations are given by a system of nonlinear hyperbolic partial differential equations (PDEs) which represent the motion of a compressible non-viscous fluid or a gas. They consist of the continuity equation, the balance of moments and the energy equation. The full set of equations is given by (see     ). Let denote the density, the velocity of the gas and the pressure. We further denote the friction coefficient of the pipe, the diameter, the cross section area. The variables of the system are , the flux
. We also denote the the speed of sound, i.e. (for constant
entropy). For natural gas we have 340 m/s. In particular, in the subsonic case ( ), the one which we consider in the sequel, two boundary conditions have to be imposed, one on the left and one on the right end of the pipe. We consider here the isothermal case only. Thus, for horizontal pipes
In the particular case, where we have a constant speed of sound , for small velocities , we arrive at the semi-linear model
1.2. Network Modeling
Let denote the graph of the gas network with vertices (nodes) an edges . Node indices are denoted , while edges are labelled . For the sake of uniqueness, we associate to each edge a direction. Accordingly, we introduce the edge-node incidence matrix
From now on we assume that all are independent of the node. We multiply (54) and (56) by λ and (55), (57) by μ and add in the following way λ(54)+ μ(27), λ(55)- μ(56). This leads to
We add the latter equations and obtain
We are going to estimate the third integral. For that matter we assume that admits a Lipschitz constant .
The second term contains quadratic expressions an mixed terms. The mixed terms need to be absorbed in the quadratics ones
We now combine (58), (59), (60) in order to obtain
We now group the corresponding quadratic expressions.
where we have used the boundary estimate due to Kato  . We now need to discuss under which configuration of the parameters the coefficients in front of the quadratic terms
can be made positive. Moreover, if , we need to absorb the corresponding boundary term again using the estimate  . It is obvious that can be positive iff and are positive. The only
parameters that we can select in order achieve positive, respectively, non- negative coefficient are the two parameters coming from the cost function and the parameters provided for the algorithm. Moreover, the coefficient becomes relevant. We recall the meaning of : it is ! So it becomes obvious from (63) that the norm of the reference solution to the adjoint equation and the Lipschitz-constant , reflecting the stiffness of the nonlinear term come into play. We thus need small to compensate the remaining terms. The question to be discussed below then is as to whether the maximum-norm of the solution of the adjoint equation which, in turn, involves is small against . Only in this case, we can choose in order to have . If, on the other side, , we have to compensate , in this case the adjoint error, by choosing sufficiently large and in order to have . The appearance of the adjoint error, in case , necessitates an a posteriori error estimate. We discuss the following cases , and :
[Case 1.] :
In this case
As mentioned above, this case involves both the reference adjoint and the adjoint error. Moreover, in this case the convergence of the iteration is determined solely by the choice of the cost parameters in that has to be sufficiently large, while plays no role.
[Case 2.] :
In this case we obtain
In this case has to absorb all negative terms and, in fact, the penalty parameter acts in the reverse sense than in Case 1. One needs to balance , and in the coefficients in order to make the two coefficients of all quadratic terms positive. We are then left with question as to whether the norm is small against for suitably large . Now, the adjoint of the full network problem has as data the right hand side and as -dependent coefficient . For a given , the corresponding network equation is linear and by Lax-Milgram’s lemma, the solution satisfies an estimate against the data, which, in turn, depend on the original solution .
[Case 3.] : In this case in (63) need to be positive. This can be achieved in general if large and small and is large compared to . A more explicit analysis can be done, but is skipped for the sake of space.
Theorem 6. Under the positivity assumptions in Cases 1, 2, 3, the iterations converge and the solutions of the iterative process (41), describing the local optimality systems on the individual edges, converge to the solution of the optimality system (40). In Case 2,3, converge to in the energy sense. In Case 1, convergence takes place in the -sense.
Example 7. We consider the following numerical example. We take for , , and and Neumann controls at all simple nodes. Clearly, the exact solution of the linear problem, i.e. , with and , is , where the adjoints have Dirichlet traces equal to 1 with Neumann traces being 0 by construction. This, however, is achieved only for very large penalty . We compute the solution with nonlinearity using the MATLAB routine bvp4 with tolerance . As for the domain decomposition, we take 15 steps. The nonlinearity is taken into account using an inner fixed-point loop, where we take 15 iterations. The parameters are chosen as 0.1, respectively. The results are shown in Figures 5-10. In Figure 11 and Figure 12, we display the numerical results for the same setup, but now with , , , , with relaxation parameter . Figure 13 reveals the fact that due to the optimality of the functions , the adjoint, as being forced to have zero Neumann data, is zero in almost the entire interval and is nontrivial in the last part only. Clearly, the three reference solutions and adjoints are plotted on top of the ddm-solutions. No difference is visible.
Figure 5. Forward solutions.
Figure 6. Adjoint solutions.
Figure 7. Optimal tripod.
Figure 8. Errors.
Figure 9. Errors of relaxed iteration.
Figure 10. Nodal errors (state,adjoint).
Figure 11. Optimal tripod.
Figure 12. Optimal state.
Figure 13. Optimal adjoint.
The author acknowledges support by the DFG TRR154 Mathematische Modellierung, Simulation und Optimierung am Beispiel von Gasnetzwerken (TPA05).
 Brouwer, J., Gasser, I. and Herty, M. (2011) Gas Pipeline Models Revisited: Model Hierarchies, Nonisothermal Models, and Simulations of Networks. Multiscale Modeling & Simulation, 9, 601-623. https://doi.org/10.1137/100813580
 Gugat, M., Leugering, G., Martin, A., Schmidt, M., Sirvent, M. and Wintergerst, D. (2016) Towards Simulation Based Mixedinteger Optimization with Differential Equations. Submitted. https://opus4.kobv.de/opus4-trr154/frontdoor/index/index/docId/63
 Gugat, M., Leugering, G., Martin, A., Schmidt, M., Sirvent, M. and Wintergerst, D. (2017) MIP-Based Instantaneous Control of Mixed-Integer PDE-Constrained Gas Transport Problems. Submitted. https://opus4.kobv.de/opus4-trr154/frontdoor/index/index/docId/140
 Hante, F., Leugering, G., Martin, A., Schewe, L. and Schmidt, M. (2017) Challenges in Optimal Control Problems for Gas and Fluid Flow in Networks of Pipes and Canals: From Modeling to Industrial Applications. In: ISIAM-Proceedings, Springer- Verlag, Berlin.
 Lagnese, J.E and Leugering, G. (2003) Domain Decomposition Methods in Optimal Control of Partial Differential Equations (Volume 148 of International Series of Numerical Mathematics). Birkh?user Verlag, Basel.
 Roubicek, T. (2013) Nonlinear Partial Differential Equations with Applications (Volume 153 of International Series of Numerical Mathematics). 2nd Edition, Birkh?user-Basel, Basel. https://doi.org/10.1007/978-3-0348-0513-1
 Kogut, P.I. and Leugering, G. (2011) Optimal Control Problems for Partial Differential Equations on Reticulated Domains, Systems and Control: Foundations and Applications. Springer, New York. https://doi.org/10.1007/978-0-8176-8149-4