Domain Decomposition of an Optimal Control Problem for Semi-Linear Elliptic Equations on Metric Graphs with Application to Gas Networks

We consider optimal control problems for the flow of gas in a pipe network. The equations of motions are taken to be represented by a semi-linear model derived from the fully nonlinear isothermal Euler gas equations. We formulate an optimal control problem on a given network and introduce a time discretization thereof. We then study the well-posedness of the corresponding time-discrete optimal control problem. In order to further reduce the complexity, we consider an instantaneous control strategy. The main part of the paper is concerned with a non-overlapping domain decomposition of the semi-linear elliptic optimal control problem on the graph into local problems on a small part of the network, ultimately on a single edge.


Introduction 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 [1] [2] [3] [4]).Let ρ denote the density, v the velocity of the gas and p the pressure.We further denote λ the friction coefficient of the pipe, D the .For natural gas we have 340 m/s.In particular, in the subsonic case ( v c < ), 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 p c ρ = , for small velocities v c  , we arrive at the semi-linear model ( )

Network Modeling
Let ( ) .For the sake of uniqueness, we associate to each edge a direction.Accordingly, we introduce the edge-node incidence matrix 1, if node is the left node of the edge , 1, if node is the right node of the edge , 0, else.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 ( ) p n coincide for all edges incident at node j n , i.e.

{ }
: 1, , 0 . We express the transmission conditions at the nodes in the following way.We introduce the edge degree : The nodal balance equation for the fluxes can be written as in instant of the From the considerations above we conclude the following system of semilinear hyperbolic equations on the metric graph G : ,0 x t T a q x t q x t c q x t p x t i x t T p x t D a h p n t q n t u t i j t T p n t p n t i k j t T d q n t j t T p x p ,0 0 , , 0 , 0, , .
To the best knowledge of the authors, for problem (5), no published result seems to be available.

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 [5] [6] [7] 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 ( ) the pairs ( ) , p q .In (7), the quantities , , , p q p q are given constants that determine the feasible pressures and flows in the pipe i , while , i i u u 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.

Time Discretization
We now consider the time discretization of (5) such that [ ] Accordingly, we denote 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 1 n t + .This strategy has is known as rolling horizon approach, the simplest case of the moving horizon paradigm, see e.g.[8] [9].Thus, for each . , , satisfies 8 at time level 1.
It is now convenient to discard the actual time level 1 n + and redefine the states at the former time as input data.To this end, we introduce 1 : 1 : and rewrite (8) as .
We now differentiate the first equation in (11) with respect to x and insert the result into the second equation.After renaming , we consider the semi-linear elliptic problem on the graph G with Neumann controls at simple nodes.
where we set ( ) ( ) . We then consider in the rest of the paper the following optimal control problem: . , , satisfies 12 .
i S i i i j p q u i j I p q u I p q x u s t p q u 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 m edges.We arrange the graph such that the central node is located at for all right ends of the edges and i x =  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 , 0 0, 0, 2, , .

Domain Decomposition
We provide an iterative non-overlapping domain decomposition that can be interpreted as an Uzawa method (Alg3, in the sense of Glowinski Then ( ) ( )  .With this notation, the general concept is easily established.We set for any Applying k  to both sides of (16), we obtain , .
Clearly, if the transmission conditions (17), (18) hold at the multiple node k n , 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 l 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 , We use the facts 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 g e q g q + + − , 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 1 l q + .We consider the following integration by parts formula after multiplying by a test function φ .
( ) ( ) e g e q g q x α φ β φ φ We now take the test function to be equal to e e e e g e q g q e x α β We use the boundary condition at the interfaces in the form , .
This identity is used in the identity (24), evaluated for the error: ( ) We obtain g e e e e g e q g q e x α β We assume and define the bilinear form We define the corresponding quadratic form applied to l e ( ) ( ) a e e e e e e x 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 l to zero.Then we obtain
Clearly, for 0 i α > , this shows that the ( )  -error strongly tends to zero.Then also the traces tend to zero as l → ∞ and, therefore, the iteration converges.
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 1 γ = 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 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

Domain Decomposition for Optimal Control Problems
We pose the following optimal control problem with Neumann boundary controls: S N 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 l → ∞ , the system (41) tends to (40).Now, (41) decomposes the fully connected problem (40) to a problem on a single edge i ∈  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.
• 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, 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 ( ) i i ρ ∈ .This provides a guess for the controls ( ) S j j u ∈ .Therefore, we can establish the states ( ) i i q ∈ .The states ( ) i i q ∈ , in turn, are inserted into the adjoint problem and that system is then solved for ( ) i i ρ ∈ , which closes the cycle.With this method, we keep the optimization in an outer loop and solve both the forward system for the states ( ) i i q ∈ and the adjoint system for ( ) i i ρ ∈ , individually.For given ( ) i i q ∈ , 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 ( ) i i q ∈ and ( ) i i ρ ∈ .The effect of parallelization can, therefore, be used for the solves in the inner loop, while the outer loop is sequential.
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 i ρ for each i ∈  then solve for i q which is introduced then in the local adjoint equations.This is then followed by the solve for i ρ and the update of the boundary data , ik ik g h 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).

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 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 [11], +   is strictly monotone and continuous and, hence, the semi-linear problem admits a unique solution q ∈  .Clearly, for regular right hand sides f , the solution is in ( ) D  .

Smoothness of the Control-to-State-Map
Let ( ) ˆt q u be the solution of ( 12) with u replaced with û tu + and let q be the solution of (12).We denote by : e q q = − the difference of these solutions.
We obtain .
e x g e q x g q x x i e

n e n i j k e n i n e n tu i n d e n k
α β Dividing by t and letting t tend to zero in (12) implies with ( )( ) ( )( ) .

g q x e x x i e n e n i j k e n i n e n u i n d e n k
For the solution q of (12), applying the standard Lax-Milgram Lemma, e′ 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 [12] in order to establish the first order optimality conditions (40).
Theorem 4.Under the assumption (33), for * f ∈  , there exists a unique solution q ∈  of ( 12).In addition, the mapping from u into q 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).

A Priori Error Estimates for the Optimality System
We denote the errors : e q q = − and :  .These errors solve the system equations: , .
We prove the following Lemma 5.The solutions , i i e p for i ∈  of (48) satisfies the estimate More precisely, for , Remark 5.1.As a result, for small data , ij ij g h , we have small solutions.
Proof of Lemma 5: We multiply the equations in (48) by 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 i e , i p by i e and i p , e g e q g q e x d e n e n e e g e q g q e x α β p g e q p g e q g q e p x d p n p n p p g e q p g e q g q p e p x Now, we reverse the roles and obtain e g e q g q p x d e n p n e p e p g e q g q e p x α β p p g e q p g e q g q d e e x d p n e n p e p e g e q p e g e q g q e e x α β ρ κ 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 p n e e g e q g q e e p e p g e q e p g e q g q e e x λβ µβ p n e n p p g e q p g e q g q p e p e p e p g e q e p x λβ µβ We add the latter equations and obtain ( g e y g q p x g e q e p g e q g q e g e q g q e p x e p e p x e p e x p n d e n p n g e q g q ρ µ ρ p e g e q g q e e p x I II III We are going to estimate the third integral.For that matter we assume that 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 [13].We now need to discuss under which configuration of the parameters the coefficients in front of the quadratic terms µ > , necessitates an a posteriori error estimate.We discuss the following cases 0, 0 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 i κ has to be sufficiently large, while i α plays no role.
In this case we obtain   In this case i α has to absorb all negative terms and, in fact, the penalty parameter i κ acts in the reverse sense than in Case 1.One needs to balance κ , δ and n in the coefficients ( ) ( ) , i c n c n in order to make the two coefficients of all quadratic terms positive.We are then left with question as to whether the norm i ρ is small against i α for suitably large i α .Now, the adjoint i ρ of the full network problem has as data the right hand side ( ) q q κ − − and as i q -dependent coefficient ( ) For a given i q , the corresponding network equation is linear and by Lax-Milgram's lemma, the solution i ρ satisfies an estimate against the data, which, in turn, depend on the original solution i q .
[Case 3.] 0, 0 λ µ > > : In this case , i i a b in (63) need to be positive.This can be achieved in general if i α large and i κ small and µ is large compared to λ .A more explicit analysis can be done, but is skipped for the sake of space.( ) , with 0 γ = and ( ) 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

=.
the cross section area.The variables of the system are ρ , the flux q a v ρ We also denote c the the speed of sound, i.e.

1 x
=  and controlled Neumann conditions at i x =  .We obtain, accordingly 1 10 tol e = − .The system with domain decomposition is solved with classical finite differences of second order.Figure1shows 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.

Figure 1 .
Figure 1.The tripod with disciniuity at the central node.

Figure 3 .
Figure 3. Error history at the central node.
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.
) describes a self-adjoint positive definite operator in the Hilbert space  :

Theorem 6 .
Under the positivity assumptions in Cases 1, 2, 3, the iterations converge and the solutions process (41), describing the local optimality systems on the individual edges, converge to the solution of the optimality system (40).In Case 2,3, , the energy sense.In Case 1, convergence takes place in the 2 L -sense.Example 7. We consider the following numerical example.We take for 10 controls at all simple nodes.Clearly, the exact solution of the linear problem, i.e.
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 0.1 γ = using the MATLAB routine bvp4 with tolerance .6 e − .

Figures 5 - 10 .
Figures 5-10.In Figure 11 and Figure 12, we display the numerical results for the same setup, but now with 1000 i α =
the graph of the gas network with vertices (nodes) To this end, it is sufficient to consider three cases: a.) the edge i connects a controlled M k ∈ at which the domain decomposition is active, c.) the edge i connects two multiple nodes , M j k ∈ .Case a.):