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
In contrast to the classical notion of discrete graphs, the graphs considered here are known as metric graphs, in the sense, that the edges are continuous curves. In fact, we consider here straight edges, along which differential equations hold. The pressure variables coincide for all edges incident at node , i.e. . We express the transmission conditions at the nodes in the following way. We introduce the edge degree . We distinguish now between multiple nodes , where , which we denote , whereas for simple nodes , for which , we write . The set of simple nodes decomposes then into those simple nodes, where Dirichlet conditions hold and Neumann nodes . Then the continuity conditions read as follows
The nodal balance equation for the fluxes can be written as in instant of the classical Kirchhoff-type condition
From the considerations above we conclude the following system of semi- linear hyperbolic equations on the metric graph :
To the best knowledge of the authors, for problem (5), no published result seems to be available.
2. Optimal Control Problems and Outline
We are now in the position to formulate optimal control problems on the level of (5). There are currently two different approaches towards optimizing and/or control the flow of gas flow through pipe networks. The first one aims at optimizing decision variables such as on-off-states for valves and compressors or zero-full-supply and demand variables for input and exit nodes, respectively. Valves and compressors can be modelled as transmission conditions at a serial node. We refer to    and refrain in the sequel from discussing issues of valves and compressors. The combined discrete and continuous optimization will be the subject of a forthcoming publication. We now describe the general format for an optimal control problem associated with the semi-linear model equations.
In (6), is a penalty parameter and a continuous function on the pairs . In (7), the quantities are given constants that determine the feasible pressures and flows in the pipe , while describe control constraints. In the continuous-time case the inequalities are considered as being satisfied for all times and everywhere along the pipes. In the sequel, we will not consider control constraints and state-constraints and, moreover, even reduce to a time semi-discretization.
We now consider the time discretization of (5) such that is decomposed into break points with widths . Accordingly, we denote . We consider a mixed implicit-explicit Euler scheme which takes in the friction term in an explicit manner.
We then obtain the optimal control problem on the time-discrete level:
In (9), we consider edgewise given cost functions e.g.
It is clear that (9) involves all time steps in the cost functional. We would like to reduce the complexity of the problem even further. To this aim we consider what has come to be known as instantaneous control. This amounts to reducing the sums in the cost function of (9) to the time-level . This strategy has is known as rolling horizon approach, the simplest case of the moving horizon paradigm, see e.g.   . Thus, for each and given , we consider the problems
It is now convenient to discard the actual time level and redefine the states at the former time as input data. To this end, we introduce , , , and rewrite (8) as
We now differentiate the first equation in (11) with respect to and insert the result into the second equation. After renaming , and introducing , we consider the semi-linear elliptic problem on the graph with Neumann controls at simple nodes.
where we set . We then consider in the rest of the paper the following optimal control problem:
Example 1. Before we embark on the non-overlapping domain decomposition method in the context of the instantaneous control paradigm (13), we look into the situation of a star graph with a central multiple node and edges. We arrange the graph such that the central node is located at for all right ends of the edges and at the left ends of the edges. This is the situation that we consider in our examples. We assumed that the first edge satisfies a homogeneous Dirichlet condition at and controlled Neumann conditions at . We obtain, accordingly
3. Domain Decomposition
We provide an iterative non-overlapping domain decomposition that can be interpreted as an Uzawa method (Alg3, in the sense of Glowinski). See the monograph  for details. The idea for this algorithm originates from a decoupling of the transmission conditions. To this end, we define the flux vector and the pressure vectors at a given node . Given a vector , we define
Then and for . With this notation, the general concept is easily established. We set for any :
Applying to both sides of (16), we obtain
But then (16) reduces to
which, in turn, implies
Clearly, if the transmission conditions (17), (18) hold at the multiple node , then (16) is also fulfilled. Thus, (16) is equivalent to the transmission conditions (17), (18). These new conditions (16) are now relaxed in an iterative scheme as follows. We use as iteration number.
We have the following relations:
This gives rise to the definition of a fixed point mapping. To this end, we need to look into the behavior of the interface in terms of , that is
We use the facts and and show
We now formulate a relaxed version of a fixed point iteration: for
Up to now, the relations concerning the iteration at the interfaces do not involve the state equation explicitly. For the analysis of the convergence of the iterates, we need to specify the equations.
The Non-Overlapping Domain Decomposition
We are interested in the errors between the solutions to the problem (12) and
Thus, we introduce . Then solves a non-linear differential equation with nonlinearity , zero right hand side and homogeneous boundary conditions at the simple nodes. As we noted above, the full transmission conditions are equivalent to (16). Hence, the error satisfies the same iterative Robin-type boundary conditions as . We consider the following integration by parts formula after multiplying by a test function .
We now take the test function to be equal to and obtain:
We use the boundary condition at the interfaces in the form
This identity is used in the identity (24), evaluated for the error:
and define the bilinear form
We define the corresponding quadratic form applied to
which is certainly bounded below by . Then the error iteration is
and, thus, the error does not increase. That it actually decreases to zero is shown next. But first we look at the relaxed version of the iteration (25). We take norms and calculate in order to obtain (for
We iterate in (36) or (37) down from to zero. Then we obtain
Clearly, for , this shows that the -error strongly tends to zero. Then also the traces tend to zero as and, therefore, the iteration converges.
Theorem 2. Under assumption (33), for each the iteration (25) with (26) and (21), (22) converges as . The convergence of the solutions is in the sense of (38).
Example 3. We show a numerical example, where three edges span a tripod. The first edge (see Figure 1) satisfies homogeneous Dirichlet conditions at the exterior node, while for the other two edges satisfy homogeneous Neumann conditions at the exterior nodes. In particular, we take , , , , . The nonlinearity is weighted by a factor and there are 10 fixed point iterations in order to handle the nonlinearity. The system without domain decomposition is solved using the MATLAB routine bvp4c with error tolerance . The system with domain decomposition is solved with classical finite differences of second order. Figure 1 shows the tripod, where we display the original solutions and the ones obtained by the domain decomposition. No difference is visible. Notice the discontinuity of the state at the central node. This is contrast to the classical nodal conditions known in the literature, where the states are continous across the multiple node, while the Neumann traces satisfying the Kirchhoff condition. We display the individual solutions―again without and with domain decomposition in Figure 2. There is no visible difference. Figure 3 shows the nodal errors at the central node. We see the nodal errors regarding the conservation of flows and the two
Figure 1. The tripod with disciniuity at the central node.
Figure 2. The three edges individually.
Figure 3. Error history at the central node.
Figure 4. Iteration history L¥-errors.
continuity conditions of the derivatives at the central node. In Figure 4, we display the relative -errors of the solutions, where the errors are taken with respect to the computed solution without domain decomposition.
4. Domain Decomposition for Optimal Control Problems
We pose the following optimal control problem with Neumann boundary controls:
The corresponding optimality system then reads as follows:
The idea now is to use a domain decomposition similar to the original system on the network. We design a method that allows to interpret the decomposed optimality system (41) as an edge-wise optimality system of an optimal control problem formulated on an individual edge. To this end, we introduce the following local system:
The same arguments that led from (16), (15) to (17), (18) apply to show that, upon convergence as , the system (41) tends to (40). Now, (41) decomposes the fully connected problem (40) to a problem on a single edge with inhomogeneous Robin-type boundary conditions. The question is as to whether the decomposed optimality system (41) is in fact itself an optimality system on that edge. If so, then it is possible to parallelize the optimization problems rather than the forward and backward solves. Let us, therefore, now consider the following optimization problems on a single edge. The idea is to introduce a virtual control that aims at controlling classical inhomogeneous Neumann condition including the iteration history at the interface as inhomogeneity to the Robin-type condition that appears in the decomposition. To this end, it is sufficient to consider three cases: a.) the edge connects a controlled Neumann node with a multiple node at which the domain decomposition is active, b.) the edge connects a controlled Neumann node with multiple node at which the domain decomposition is active, c.) the edge connects two multiple nodes .
・ If we write down the optimality systems for (42), (43) and (44), respectively, and combine the results, we arrive at (41).
・ This shows that within the loop of iterations that restore the transmission conditions at the multiple nodes, we can reformulate the system (41) as the optimality system of an optimal control problem formulated on a single edge, with input data coming from the iteration history that involves all nodes adjacent at the ends of the given edge.
・ This means that we can actually decompose the optimization problem given on the graph into a sequence of local optimization problems given on an individual edge.
・ The resulting optimization problem on the individual edges are strictly convex, thus, admitting a unique global solution.
Remark 4.2. There are at least two ways to use the proposed ddm-approach.
1) In the first approach, we consider (40) and start with a guess for the adjoint variables . This provides a guess for the controls . Therefore, we can establish the states . The states , in turn, are inserted into the adjoint problem and that system is then solved for , which closes the cycle. With this method, we keep the optimization in an outer loop and solve both the forward system for the states and the adjoint system for , individually. For given , the adjoint system is a linear elliptic problem on the graph. To this system the ddm-method above applies and converges. As we have established above, the forward problem admits a convergent ddm-algorithm. This finally means that in the inner loop we can use convergent ddm-iterations for finding and . The effect of parallelization can, therefore, be used for the solves in the inner loop, while the outer loop is sequential.
2) In the second approach, we decompose the coupled system (40) to (41). The resulting decoupled problem is then the optimality system for the virtual optimal control problems (42), (43) or (44), as seen above. In this case, there is no outer loop other than the ddm-iteration which is completely parallel. Still, the local optimality systems have to be solved in a way describes in the first approach. Namely, we provide an initial guess for for each then solve for which is introduced then in the local adjoint equations. This is then followed by the solve for and the update of the boundary data which are used in the communication at the next ddm-iteration. In this, admittedly, more elegant approach, the constrained minimization problem on the entire graph can be decomposed to minimization problems on a single edge. As we will see below, unfortunately, but expectedly, the convergence is no longer global as in the first approach, but rather local. This means that only if we start close to a solution of (40), or if we have a priori estimates and tune the parameters accordingly, we can prove convergence of the unique solutions of (41) to those of (40).
5.1. Wellposedness of the Primal Problem
The semi-linear network problem (12) admits a unique solution. This is true, as the linear part of problem (12) describes a self-adjoint positive definite operator in the Hilbert space :
Indeed, we also define the energy space
and the operator as follows:
It is a matter of applying standard integration by parts to show that, indeed, is symmetric and positive definite in such that can be extended as a self-adjoint operator in . Then it is standard to show that can be extended to a bounded coercive map from into its dual . If we assume (33) and define the Nemitskji operator then is strictly monotone and continuous. Hence, according to  , is strictly monotone and continuous and, hence, the semi-linear problem admits a unique solution . Clearly, for regular right hand sides , the solution is in .
5.2. Smoothness of the Control-to-State-Map
Let be the solution of (12) with replaced with and let be the solution of (12). We denote by the difference of these solutions. We obtain
Dividing by and letting tend to zero in (12) implies with
For the solution of (12), applying the standard Lax-Milgram Lemma, uniquely solves the linear elliptic network problem (47) and, therefore, satisfies standard energy estimates. As the cost function in (39) is convex, according to the classical Weierstrass theorem, problem (39) admits a unique solution. One can then verify the conditions for the Ioffe-Tichomirov Theorem  in order to establish the first order optimality conditions (40).
Theorem 4. Under the assumption (33), for , there exists a unique solution of (12). In addition, the mapping from into is Gateaux differentiable. Moreover, the optimal control problem (39) admits a unique solution. The optimal solution is characterized by the optimality system of first order (40).
5.3. A Priori Error Estimates for the Optimality System
We denote the errors and for and . These errors solve the system equations:
We prove the following
Lemma 5. The solutions for of (48) satisfies the estimate
More precisely, for , we obtain
Remark 5.1. As a result, for small data , we have small solutions.
Proof of Lemma 5: We multiply the equations in (48) by and , respectively, integrate and then use integration by parts. For the sake of brevity, we leave the full arguments to the reader.
We multiply the state equation for the errors , by and , respectively.
Now, we reverse the roles and obtain
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).