AM  Vol.8 No.8 , August 2017
Domain Decomposition of an Optimal Control Problem for Semi-Linear Elliptic Equations on Metric Graphs with Application to Gas Networks
Abstract
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.

1. Introduction

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 [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 diameter, a 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. c 2 = p ρ (for constant

entropy). 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

ρ t + x ( ρ v ) = 0 t ( ρ v ) + x ( p + ρ v 2 ) = λ 2 D ρ v | v | . (1)

In the particular case, where we have a constant speed of sound c = p ρ , for small velocities | v | c , we arrive at the semi-linear model

ρ t + x ( ρ v ) = 0 t ( ρ v ) + p x = λ 2 D ρ v | v | . (2)

1.2. Network Modeling

Let G = ( V , E ) denote the graph of the gas network with vertices (nodes) V = { n 1 , n 2 , , n | V | } an edges E = { e 1 , e 2 , , e | E | } . Node indices are denoted j J , | J | = | V | , while edges are labelled i I , | I | = | E | . For the sake of uniqueness, we associate to each edge a direction. Accordingly, we introduce the edge-node incidence matrix

d i j = { 1, if e i + q i ) g i ( q i ) ) ρ i d + κ i e i ) e i d x = k J i I k d i k β i x p i ( n k ) e i ( n k ) + i I 0 l i { α i p i e i + β i x p i x e i + g i ( e i + q i ) p i e i + ( g i ( e i + q i ) g i ( q i ) ) ρ i e i + κ i e i 2 } d x . (57)

From now on we assume that all λ k = : λ , μ k : = μ 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

k J i I k d i k e i ( n k ) ( λ β i x e i ( n k ) + μ β i x p i ( n k ) ) = i I 0 l i λ { α i ( e i ) 2 + β i ( x e i ) 2 + ( g i ( e i + q i ) g i ( q i ) ) e i } + μ { α i e i p i + β i x e i x p i + g i ( e i + q i ) e i p i + ( g i ( e i + q i ) g i ( q i ) ) ρ i e i + κ i e i 2 } d x

k J i I k d i k p i ( n k ) ( λ β i x p i ( n k ) μ β i x e i ( n k ) ) = i I 0 l i λ { α i ( p i ) 2 + β i ( x p i ) 2 + g i ( e i + q i ) p i 2 + ( g i ( e i + q i ) g i ( q i ) ) ρ i p i + κ i e i p i } μ { α i e i p i + β i x e i x p i + g i ( e i + q i ) e i p i } d x

We add the latter equations and obtain

k J i I k d i k e i ( n k ) ( λ β i x e i ( n k ) + μ β i x p i ( n k ) ) + k J i I k d i k p i ( n k ) ( λ β i x p i ( n k ) μ β i x e i ( n k ) ) = i I 0 l i λ { α i e i 2 + α i p i 2 + β i ( x e i ) 2 + β i ( x p i ) 2 } d x + i I 0 l i κ i ( μ e i 2 + λ p i e i ) d x + 1 ν k J M i I k ( μ ( p i ( n k ) ) 2 λ d i k e i ( n k ) p i ( n k ) ) + λ i I 0 l i { ( g i ( e i + q i ) g i ( q i ) ) e i + g i ( e i + q i ) p i 2 + ( g i ( e i + y i ) g i ( q i ) ) ρ i p i } d x + μ i I 0 l i { g i ( e i + q i ) e i p i + ( g i ( e i + q i ) g i ( q i ) ) ρ i e i ( g i ( e i + q i ) g i ( q i ) ) e i p i } d x i I 0 l i λ { α i ( e i 2 + p i 2 ) + β i ( x e i ) 2 + β i ( x p i ) 2 } d x + i I 0 l i κ i ( μ e i 2 + λ p i e i ) d x + 1 ν k J M i I k ( μ ( p i ( n k ) ) 2 λ d i k e i ( n k ) p i ( n k ) ) + i I 0 l i ( g i ( e i + q i ) g i ( q i ) ) ρ i ( λ p i + μ e i ) + μ ( g i ( e i + q i ) g i ( q i + θ e i ) ) e i p i d x = I + I I + I I I . (58)

We are going to estimate the third integral. For that matter we assume that g i ( s ) admits a Lipschitz constant L i .

I I I i I 0 l i L i | ρ i | | e i | ( λ | p i | + μ | e i | ) + μ L i ( 1 θ ) | p i | ( e i ) 2 d x i I L i { ( ( μ + λ 2 ) ρ i L ( 0 , l i ) + μ p i L ( 0 , l i ) ) 0 l i e i 2 d x + λ 2 ρ i L ( 0 , l i ) 0 l i p i 2 d x } (59)

The second term contains quadratic expressions an mixed terms. The mixed terms need to be absorbed in the quadratics ones

i I 0 l i κ i ( μ e i 2 + λ p i e i ) d x i I 0 l i κ i ( ( μ λ δ ) e i 2 λ δ p i 2 ) d x 1 ν k J M i I k ( μ ( p i ( n k ) ) 2 λ d i k e i ( n k ) p i ( n k ) ) 1 ν k J M i I k ( ( μ λ δ ) ( p i ( n k ) ) 2 λ δ e i ( n k ) 2 ) (60)

We now combine (58), (59), (60) in order to obtain

k J i I k d i k e i ( n k ) ( λ β i x e i ( n k ) + μ β i x p i ( n k ) ) + k J i I k d i k p i ( n k ) ( λ β i x p i ( n k ) μ β i x e i ( n k ) ) i I 0 l i λ { α i ( e i 2 + p i 2 ) + β i ( x e i ) 2 + β i ( x p i ) 2 } d x + i I 0 l i κ i ( ( μ λ δ ) e i 2 λ δ p i 2 ) d x + 1 ν k J M i I k ( ( μ λ δ ) ( p i ( n k ) ) 2 λ δ e i ( n k ) 2 ) i I L i { ( ( μ + λ 2 ) ρ i L ( 0 , l i ) + μ p i L ( 0 , l i ) ) 0 l i e i 2 d x + λ 2 ρ i L ( 0 , l i ) 0 l i p i 2 d x } (61)

We now group the corresponding quadratic expressions.

k J i I k d i k e i ( n k ) ( λ β i x e i ( n k ) + μ β i x p i ( n k ) ) + k J i I k d i k p i ( n k ) ( λ β i x p i ( n k ) μ β i x e i ( n k ) ) i I ( λ α i + κ i ( μ λ δ ) L i ( ( μ + λ 2 ) ρ i L ( 0 , l i ) + μ p i L ( 0 , l i ) ) λ δ c i 2 ( n ) ) 0 l i e i 2 d x + ( λ α i λ δ L i λ 2 ρ i L ( 0 , l i ) ) 0 l i p i 2 d x + i I 0 l i { λ ( 1 δ c 1 ( n ) ) ( x e i ) 2 + λ ( x p i ) 2 } d x + 1 ν k J M i I k ( ( μ λ δ ) ( p i ( n k ) ) 2 ) , (62)

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

#Math_252# (63)

can be made positive. Moreover, if μ λ δ 0 , we need to absorb the corresponding boundary term again using the estimate [13] . It is obvious that b i can be positive iff ( α i δ L i 1 2 ρ i L ( 0, l i ) ) and λ are positive. The only

parameters that we can select in order achieve positive, respectively, non- negative coefficient are the two parameters κ i , ν i coming from the cost function and the parameters λ 0, μ 0 provided for the algorithm. Moreover, the coefficient α i > 0 becomes relevant. We recall the meaning of α i : it is 1 Δ t ! So it becomes obvious from (63) that the norm of the reference solution to the adjoint equation ρ i and the Lipschitz-constant L i , reflecting the stiffness of the nonlinear term come into play. We thus need small Δ t > 0 to compensate the remaining terms. The question to be discussed below then is as to whether the maximum-norm of the solution ρ i of the adjoint equation which, in turn, involves α i is small against α i . Only in this case, we can choose λ > 0 in order to have b i > 0 . If, on the other side, λ = 0 , we have to compensate L i p i , in this case the adjoint error, by choosing κ i sufficiently large and μ > 0 in order to have a i > 0 . The appearance of the adjoint error, in case μ > 0 , necessitates an a posteriori error estimate. We discuss the following cases λ = 0 , μ > 0 , λ > 0 , μ = 0 and λ > 0 , μ > 0 :

[Case 1.] λ = 0 , μ > 0 :

a i = μ ( κ i L i ( ρ i L ( 0 , l i ) + p i L ( 0 , l i ) ) ) , b i : = 0 (64)

In this case

k J i I k d i k e i ( n k ) ( λ β i x e i ( n k ) + μ β i x p i ( n k ) ) . + k J i I k d i k p i ( n k ) β i ( λ x p i ( n k ) μ x e i ( n k ) ) i I μ ( κ i L i ( ρ i L ( 0 , l i ) + p i L ( 0 , l i ) ) ) 0 l i e i 2 d x + μ ν k J M i I k ( p i ( n k ) ) 2 . (65)

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.

[Case 2.] λ > 0 , μ = 0 :

a i : = λ ( α i κ i 1 δ L i ( 1 2 ρ i L ( 0 , l i ) ) δ c i 2 ( n ) ) b i : = λ ( α i δ L i 1 2 ρ i L ( 0 , l i ) ) . (66)

In this case we obtain

k J i I k d i k e i ( n k ) ( λ β i x e i ( n k ) + μ β i x p i ( n k ) )

+ k J i I k d i k p i ( n k ) β i ( λ x p i ( n k ) μ x e i ( n k ) ) (67)

λ i I { ( α i κ i 1 δ L i 1 2 ρ i L ( 0, l i ) δ c i 2 ( n ) ) 0 l i e i 2 d x

+ ( α i δ L i 1 2 ρ i L ( 0, l i ) δ c i 2 ( n ) ) 0 l i p i 2 d x (68)

+ i I 0 l i ( 1 δ c 1 ( n ) ) ( x e i ) 2 + ( 1 δ c 1 ( n ) ) ( x p i ) 2 d x }

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 c 1 ( n ) , c i 2 ( 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 κ i ( q i q i 0 ) and as q i -dependent coefficient g i ( q i ) . For a given q i , 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 q i .

[Case 3.] λ > 0 , μ > 0 : In this case a i , b i 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.

Theorem 6. Under the positivity assumptions in Cases 1, 2, 3, the iterations converge and the solutions q l = ( q i l ) i I 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, q i l , p i l converge to q i 0 , p i 0 in the energy sense. In Case 1, convergence takes place in the L 2 -sense.

Example 7. We consider the following numerical example. We take for α = 10 , κ = 10 , ν = 1 and f i = α x and Neumann controls at all simple nodes. Clearly, the exact solution of the linear problem, i.e. γ g ( q ) , with γ = 0 and g ( q ) = | q | q , is q i = x , i = 1 , 2 , 3 , 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 γ = 0.1 using the MATLAB routine bvp4 with tolerance e . 6 . 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 α i = 1000 , κ i = 100 , μ = 1 , λ = 0 , γ = 0.1 with relaxation parameter ϵ = 0.1 . Figure 13 reveals the fact that due to the optimality of the functions q i ( x ) = x , 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.

Acknowledgements

The author acknowledges support by the DFG TRR154 Mathematische Modellierung, Simulation und Optimierung am Beispiel von Gasnetzwerken (TPA05).

Cite this paper
Leugering, G. (2017) Domain Decomposition of an Optimal Control Problem for Semi-Linear Elliptic Equations on Metric Graphs with Application to Gas Networks. Applied Mathematics, 8, 1074-1099. doi: 10.4236/am.2017.88082.
References

[1]   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

[2]   LeVeque, R.J. (1992) Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge.

[3]   LeVeque, R.J. (1992) Numerical Methods for Conservation Laws. Birkh?user-Verlag, Basel.

[4]   Smoller, J. (1994) Shock Waves and Reaction—Diffusion Equations (Volume 258 of Grundlehren der Mathematischen Wissenschaften). Springer-Verlag, Berlin.

[5]   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

[6]   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

[7]   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.
https://opus4.kobv.de/opus4-trr154/frontdoor/index/index/docId/121

[8]   Hinze, M. and Volkwein, S. (2002) MIP-Based Instantaneous Control of Mixed- Integer PDE-Constrained Gas Transport Problems. Nonlinear Anal., 50, 1-26.

[9]   Hinze, M. and Volkwein, S. (2001) Instantaneous Control of Vibrating String Networks. In: Online Optimization of Large Scale Systems, Springer-Verlag, Berlin, 229-249.

[10]   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.

[11]   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

[12]   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

[13]   Kato, T. (1976) Perturbation Theory for Linear Operators. 2nd Edition, Springer- Verlag, Berlin-New York, Grundlehren der Mathematischen Wissenschaften, Band 132.

 
 
Top