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

Günter Leugering

Chair of Applied Mathematics II, Friedrich-Alexander-University Erlangen-Nuremberg, Erlangen, Germany.

DOI: 10.4236/am.2017.88082                

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.

Keywords

Optimal Control, Gas Networks, Euler’s Equation, Hierarchy of Models, Semi-Linear Approximation, Non-Overlapping Domain Decomposition

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 $\rho$ denote the density, $v$ the velocity of the gas and $p$ the pressure. We further denote $\lambda$ the friction coefficient of the pipe, $D$ the diameter, $a$ the cross section area. The variables of the system are $\rho$ , the flux

$q=a\rho v$ . We also denote $c$ the the speed of sound, i.e. $c^2=\frac{\partial p}{\partial \rho }$ (for constant

entropy). For natural gas we have 340 m/s. In particular, in the subsonic case ( $\left|v\right|<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

$\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x}\left(\rho v\right)=0\\ \frac{\partial }{\partial t}\left(\rho v\right)+\frac{\partial }{\partial x}\left(p+\rho v^2\right)=-\frac{\lambda }{2D}\rho v\left|v\right|\mathrm{.}\end{array}$ (1)

In the particular case, where we have a constant speed of sound $c=\sqrt{\frac{p}{\rho }}$ , for small velocities $\left|v\right|\ll c$ , we arrive at the semi-linear model

$\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x}\left(\rho v\right)=0\\ \frac{\partial }{\partial t}\left(\rho v\right)+\frac{\partial p}{\partial x}=-\frac{\lambda }{2D}\rho v\left|v\right|.\end{array}$ (2)

1.2. Network Modeling

Let $G=\left(V,E\right)$ denote the graph of the gas network with vertices (nodes) $V=\left\{n_1,n_2,\cdots ,n_{\left|V\right|}\right\}$ an edges $E=\left\{e_1,e_2,\cdots ,e_{\left|E\right|}\right\}$ . Node indices are denoted $j\in \mathcal{J},\left|\mathcal{J}\right|=\left|V\right|$ , while edges are labelled $i\in \mathcal{I}\mathrm{,}\left|\mathcal{I}\right|=\left|E\right|$ . For the sake of uniqueness, we associate to each edge a direction. Accordingly, we introduce the edge-node incidence matrix

$d_{ij}=\{\begin{array}{l}-\mathrm{1,}\text{if}\text{node}{n}_j\text{is}\text{the}\text{left}\text{node}\text{of}\text{the}\text{edge}{e}_i\mathrm{,}\\ +\mathrm{1,}\text{if}\text{node}{n}_j\text{is}\text{the}\text{right}\text{node}\text{of}\text{the}\text{edge}{e}_i\mathrm{,}\\ \mathrm{0,}\text{else}\mathrm{.}\end{array}$

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_i\left(n_j\right)$ coincide for all edges incident at node $n_j$ , i.e. $i\in {\mathcal{I}}_j\mathrm{<mo>\; :\; </mo>}=\left\{i\in \mathrm{1,}\cdots ,E|d_{ij}\ne 0\right\}$ . We express the transmission conditions at the nodes in the following way. We introduce the edge degree $d_j\mathrm{<mo>\; :\; </mo>}=\left|{\mathcal{I}}_j\right|$ . We distinguish now between multiple nodes $n_j$ , where $d_j>1$ , which we denote ${\mathcal{J}}^M$ , whereas for simple nodes $n_j$ , for which $d_i=1$ , we write ${\mathcal{J}}^S$ . The set of simple nodes decomposes then into those simple nodes, where Dirichlet conditions hold ${\mathcal{J}}_D^S$ and Neumann nodes ${\mathcal{J}}_N^S$ . Then the continuity conditions read as follows

$p_i\left(n_j\mathrm{,}t\right)=p_k\left(n_j\mathrm{,}t\right)\mathrm{,}\forall i\mathrm{,}k\in {\mathcal{I}}_j\mathrm{,}j\in {\mathcal{J}}^M\mathrm{,}t\in \left(\mathrm{0,}T\right)\mathrm{.}$ (3)

The nodal balance equation for the fluxes can be written as in instant of the classical Kirchhoff-type condition

${\displaystyle \underset{i\in {\mathcal{I}}_j}{\sum }}{d}_{ij}{q}_i\left(n_j\mathrm{,}t\right)=\mathrm{0,}j\in {\mathcal{J}}^M\mathrm{,}t\in \left(\mathrm{0,}T\right)\mathrm{.}$ (4)

From the considerations above we conclude the following system of semi- linear hyperbolic equations on the metric graph $G$ :

$\begin{array}{l}{\partial }_t{p}_i\left(x,t\right)+\frac{c_i^2}{a_i}{\partial }_x{q}_i\left(x,t\right)=0,i\in \mathcal{I},\left(x,t\right)\in \left(0,{\mathcal{l}}_i\right)\times \left(0,T\right)\\ {\partial }_t{q}_i\left(x,t\right)+{\partial }_x{p}_i\left(x,t\right)=-\frac{\lambda c_i^2}{2D_i{a}_i^2}\frac{q_i\left(x,t\right)\left|q_i\left(x,t\right)\right|}{p_i\left(x,t\right)},i\in \mathcal{I},\left(x,t\right)\in \left(0,{\mathcal{l}}_i\right)\times \left(0,T\right)\\ h_j\left(p_i\left(n_j,t\right),q_i\left(n_j,t\right)\right)=u_j\left(t\right),i\in {\mathcal{I}}_j,j\in {\mathcal{J}}^S,t\in \left(0,T\right)\\ p_i\left(n_j,t\right)=p_k\left(n_j,t\right),\forall i,k\in {\mathcal{I}}_j,j\in {\mathcal{J}}^M,t\in \left(0,T\right)\\ {\displaystyle \underset{i\in {\mathcal{I}}_j}{\sum }}{d}_{ij}{q}_i\left(n_j,t\right)=0,j\in {\mathcal{J}}^M,t\in \left(0,T\right)\\ p_i\left(x,0\right)=p_{i,0}\left(x\right),q_i\left(x,0\right)=q_{i0}\left(x\right),x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}.\end{array}$ (5)

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

$\begin{array}{l}\underset{\left(p\mathrm{,}q\mathrm{,}u\right)\in \Xi }{min}I\left(p\mathrm{,}q\mathrm{,}u\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{T}{\int }}}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{I}_i\left(p_i\mathrm{,}{q}_i\right)\text{d}x\text{d}t+\frac{\nu }{2}{\displaystyle \underset{j\in {\mathcal{J}}^S}{\sum }}{\displaystyle \underset{0}{\overset{T}{\int }}}{\left|u_j\left(t\right)\right|}^2\text{d}t\\ s\mathrm{.}t\mathrm{.}\left(p\mathrm{,}q\mathrm{,}u\right)\text{satisfies}\left(\text{5}\right)\mathrm{,}\end{array}$ (6)

$\Xi \mathrm{<mo>\; :\; </mo>}=\left\{\left(p\mathrm{,}q\mathrm{,}u\right)\mathrm{<mo>\; :\; </mo>}{\underset{\_}{p}}_i\le p_i\le {\bar{p}}_i\mathrm{,}{\underset{\_}{q}}_i\le q_i\le {\bar{q}}_i\mathrm{,}{\underset{\_}{u}}_i\le u_i\le {\bar{u}}_i\mathrm{,}i\in \mathcal{I}\right\}\mathrm{.}$ (7)

In (6), $\nu >0$ is a penalty parameter and $I_i\left(\cdot \mathrm{,}\cdot \right)$ a continuous function on the pairs $\left(p\mathrm{,}q\right)$ . In (7), the quantities ${\underset{\_}{p}}_i\mathrm{,}{\underset{\_}{q}}_i\mathrm{,}{\bar{p}}_i\mathrm{,}{\bar{q}}_i$ are given constants that determine the feasible pressures and flows in the pipe $i$ , while ${\underset{\_}{u}}_i\mathrm{,}{\bar{u}}_i$ 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 $\left[\mathrm{0,}T\right]$ is decomposed into break points $t_0=0<t_1<\cdots <t_N=T$ with widths $\Delta t_n:=t_{n+1}-t_n,n=0,\cdots ,N-1$ . Accordingly, we denote $p_i\left(x,t_n\right)=:p_{i,n}\left(x\right),q_i\left(x,t_n\right)=:q_{i,n}\left(x\right),n=0,\cdots ,N-1$ . We consider a mixed implicit-explicit Euler scheme which takes $p_i$ in the friction term in an explicit manner.

$\begin{array}{l}\frac{1}{\Delta t}{p}_{i,n+1}\left(x\right)+\frac{c_i^2}{a_i}{\partial }_x{q}_{i,n+1}\left(x\right)=\frac{1}{\Delta t}{p}_{i,n}\left(x\right),x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ \frac{1}{\Delta t}{q}_{i,n+1}\left(x\right)+{\partial }_x{p}_{i,n+1}\left(x\right)\\ =-\frac{\lambda c_i^2}{2D_i{a}_i^2}\frac{q_{i,n+1}\left(x\right)\left|q_{i,n+1}\left(x\right)\right|}{p_{i,n}\left(x\right)}+\frac{1}{\Delta t}{q}_{i,n}\left(x\right),x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ g_j\left(p_{i,n+1}\left(n_j\right),q_{i,n+1}\left(n_j\right)\right)=u_{j,n+1},i\in {\mathcal{I}}_j,j\in {\mathcal{J}}^S\\ p_{i,n+1}\left(n_j\right)=p_{k,n+1}\left(n_j\right),\forall i,k\in {\mathcal{I}}_j,j\in {\mathcal{J}}^M\\ {\displaystyle \underset{i\in {\mathcal{I}}_j}{\sum }}{d}_{ij}{q}_{i,n+1}\left(n_j\right)=0,j\in {\mathcal{J}}^M,\\ p_{i,0}\left(x\right)=p_{i,0}\left(x\right),q_{i,0}\left(x\right)\\ =q_{i0}\left(x\right),x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}.\end{array}$ (8)

We then obtain the optimal control problem on the time-discrete level:

$\begin{array}{l}\underset{\left(p\mathrm{,}q\mathrm{,}u\right)}{min}I\left(p\mathrm{,}q\mathrm{,}u\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{n=1}{\overset{N}{\sum }}}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{I}_i\left(p_{i\mathrm{,}n}\mathrm{,}{q}_{i\mathrm{,}n}\right)\text{d}x+\frac{\nu }{2}{\displaystyle \underset{j\in {\mathcal{J}}^S}{\sum }}{\displaystyle \underset{n\mathrm{<mo>\; =\; </mo>1}}{\overset{N}{\sum }}}{\left|u_j\left(n\right)\right|}^2\\ s\mathrm{.}t\mathrm{.}\left(p\mathrm{,}q\mathrm{,}u\right)\text{satisfies}\left(8\right)\mathrm{.}\end{array}$ (9)

In (9), we consider edgewise given cost functions e.g.

$I_i\left(p_{i,n},q_{i,n}\right)\left(x\right):=\frac{{\kappa }_i}{2}\left\{{\left|p_{i,n}\left(x\right)-p_{i,n}^d\left(x\right)\right|}^2+{\left|q_{i,n}\left(x\right)-q_{i,n}^d\left(x\right)\right|}^2\right\},x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}.$

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 $t_{n+1}$ . 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 $n=1,\cdots ,N-1$ and given $p_{i\mathrm{,}n}\mathrm{,}{q}_{i\mathrm{,}n}$ , we consider the problems

$\begin{array}{l}\underset{\left(p\mathrm{,}q\mathrm{,}u\right)}{min}I\left(p\mathrm{,}q\mathrm{,}u\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{I}_i\left(p_i\mathrm{,}{q}_i\right)\text{d}x+\frac{\nu }{2}{\displaystyle \underset{j\in {\mathcal{J}}^S}{\sum }}{\left|u_j\right|}^2\\ s\mathrm{.}t\mathrm{.}\left(p\mathrm{,}q\mathrm{,}u\right)\text{satisfies}\left(8\right)\text{at}\text{time}\text{level}n+1.\end{array}$ (10)

It is now convenient to discard the actual time level $n+1$ and redefine the states at the former time as input data. To this end, we introduce ${\alpha }_i:=\frac{1}{\Delta t}$ , $f_i^1:=\frac{1}{\Delta t}{p}_{i,n}\left(x\right)$ , $f_i^2:=\frac{1}{\Delta t}{q}_{i,n}\left(x\right)$ , ${\gamma }_i\left(x\right):=\frac{\lambda c_i^2}{2D_i{a}_i^2}\frac{1}{p_{i,n}\left(x\right)}$ and rewrite (8) as

$\begin{array}{l}{\alpha }_i{p}_i\left(x\right)+\frac{c_i^2}{a_i}{\partial }_x{q}_i\left(x\right)=f_i^1,x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ {\alpha }_i{q}_i\left(x\right)+{\partial }_x{p}_i\left(x\right)+{\gamma }_i\left(x\right)q_i\left(x\right)\left|q_i\left(x\right)\right|=f_i^2,x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ g_j\left(p_i\left(n_j\right),q_i\left(n_j\right)\right)=u_j,i\in {\mathcal{I}}_j,j\in {\mathcal{J}}^S\\ p_i\left(n_j\right)=p_k(n_j),\forall i,k\in {\mathcal{I}}_j,j\in {\mathcal{J}}^M\\ {\displaystyle \underset{i\in {\mathcal{I}}_j}{\sum }}{d}_{ij}{q}_i\left(n_j\right)=0,j\in {\mathcal{J}}^M.\end{array}$ (11)

We now differentiate the first equation in (11) with respect to $x$ and insert the result into the second equation. After renaming $f_i:=f_i^2-\frac{1}{{\alpha }_i}{\partial }_x{f}_i^1$ , $\mathcal{I}=\left\{i=1,\cdots ,m\right\}$ and introducing ${\beta }_i=\frac{c_i^2}{a_i{\alpha }_i}$ , we consider the semi-linear elliptic problem on the graph $G$ with Neumann controls at simple nodes.

$\begin{array}{l}{\alpha }_i{q}_i\left(x\right)-{\beta }_i{\partial }_{xx}{q}_i\left(x\right)+g_i\left(x;q_i\left(x\right)\right)=f_i\left(x\right),x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ {\beta }_i{\partial }_x{q}_i\left(n_k\right)={\beta }_j{\partial }_x{q}_j\left(n_k\right),i\ne j\in {\mathcal{I}}_k,k\in {\mathcal{J}}^M\\ q_i\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_D^S\\ {\partial }_x{q}_i\left(n_k\right)=u_k,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_N^S\\ {\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{q}_i\left(n_k\right)=0,k\in {\mathcal{J}}^M,\end{array}$ (12)

where we set $g_i\left(x;s\right):={\gamma }_i\left(x\right)s\left|s\right|$ . We then consider in the rest of the paper the following optimal control problem:

$\begin{array}{l}\underset{\left(p,q,u\right)}{\min }I\left(p,q,u\right):={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{I}_i\left(p_i,q_i\right)\text{d}x+\frac{\nu }{2}{\displaystyle \underset{j\in {\mathcal{J}}^S}{\sum }}{\left|u_j\right|}^2\\ s.t.\left(p,q,u\right)\text{satisfies}\left(12\right).\end{array}$ (13)

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 $x=0$ for all right ends of the edges and $x={\mathcal{l}}_i$ 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 $x={\mathcal{l}}_1$ and controlled Neumann conditions at $x={\mathcal{l}}_i$ . We obtain, accordingly

$\begin{array}{l}{\alpha }_i{q}_i\left(x\right)-{\beta }_i{\partial }_{xx}{q}_i\left(x\right)+{\alpha }_i{\gamma }_i\left(x\right)q_i\left(x\right)\left|q_i\left(x\right)\right|=f_i\left(x\right),x\in \left(0,{\mathcal{l}}_i\right),i=1,\cdots ,m\\ {\beta }_i{\partial }_x{q}_i\left(0\right)={\beta }_j{\partial }_x{q}_j\left(0\right),i\ne j=1,\cdots ,m\\ {\displaystyle \underset{i=1}{\overset{m}{\sum }}}{q}_i\left(0\right)=0,q_1\left({\mathcal{l}}_1\right)=0,{\partial }_x{q}_i\left({\mathcal{l}}_i\right)=u_{i}i=2,\cdots ,m.\end{array}$ (14)

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 [10] for details. The idea for this algorithm originates from a decoupling of the transmission conditions. To this end, we define the flux vector $q^k\mathrm{<mo>\; :\; </mo>}={\left(d_{ik}{q}_i\left(n_k\right)\mathrm{,}i\in {\mathcal{I}}_k\right)}^{\text{T}}$ and the pressure vectors $\partial q^k\mathrm{<mo>\; :\; </mo>}={\left({\beta }_i{\partial }_x{q}_i\left(n_k\right)\mathrm{,}i\in {\mathcal{I}}_k\right)}^{\text{T}}$ at a given node $n_k\mathrm{,}k\in {\mathcal{J}}^M$ . Given a vector $z:=\left(z_i,i\in {\mathcal{I}}_k\right)$ , we define

${\left({\mathcal{S}}^k\left(z\right)\right)}_i\mathrm{<mo>\; :\; </mo>}=\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{z}_j-z_i\mathrm{.}$ (15)

Then ${\left(S^k\right)}^2=I$ and $S^k\left(1\right)=1$ for $\mathrm{1\; <mo>\; :\; </mo>}=\left(\mathrm{1,}\cdots \mathrm{,1}\right)\in {ℝ}^{d_k}$ . With this notation, the general concept is easily established. We set for any $\sigma >0$ :

$q^k+\sigma \partial q^k=\sigma {\mathcal{S}}^k\left(\partial q^k\right)-{\mathcal{S}}^k\left(q^k\right)\mathrm{.}$ (16)

Applying ${\mathcal{S}}^k$ to both sides of (16), we obtain

${\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{q}_i\left(n_k\right)=0.$ (17)

But then (16) reduces to

$\partial q_i\left(n_k\right)=\frac{1}{d_k}{\displaystyle \underset{k\in {\mathcal{I}}_k}{\sum }}\partial q_j\left(n_k\right)\mathrm{,}i\in {\mathcal{I}}_k\mathrm{,}$

which, in turn, implies

${\beta }_i{\partial }_x{q}_i\left(n_k\right)={\beta }_j{\partial }_x{q}_j\left(n_k\right)\mathrm{,}i\ne j\in {\mathcal{I}}_k\mathrm{,}k\in {\mathcal{J}}^M\mathrm{.}$ (18)

Clearly, if the transmission conditions (17), (18) hold at the multiple node $n_k$ , 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.

${\left(p^k\right)}^{l+1}+\sigma \partial {\left(q^k\right)}^{l+1}=\sigma {\mathcal{S}}^k\left({\left(\partial q^k\right)}^l\right)-{\mathcal{S}}^k\left({\left(q^k\right)}^l\right)=\mathrm{<mo>\; :\; </mo>}{\left(g^k\right)}^{l+1}\mathrm{.}$ (19)

We have the following relations:

${\left(g^k\right)}^{l+1}={\mathcal{S}}^k\left(2\sigma \partial {\left(q^k\right)}^l-{\left(g^k\right)}^l\right)\mathrm{.}$ (20)

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 $g^k\mathrm{,}k\in {\mathcal{J}}^M$ , that is

$g\in \mathcal{X}:={\Pi }_{k\in {\mathcal{J}}^M}{\Pi }_{i\in {\mathcal{I}}_k},{\Vert g\Vert }_{\mathcal{X}}^2:={\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\frac{1}{\sigma }{\left|g_{i,k}\right|}^2,$ (21)

$\mathcal{T}\mathrm{<mo>\; :\; </mo>}\mathcal{X}\to \mathcal{X}\mathrm{,}$ (22)

${\left(\mathcal{T}g\right)}_{i\mathrm{,}k}\mathrm{<mo>\; :\; </mo>}={\mathcal{S}}^k{\left(2\sigma \partial \left(q^k\right)-g^k\right)}_i\mathrm{,}k\in {\mathcal{J}}^M\mathrm{,}i\in {\mathcal{I}}_k\mathrm{,}$

${\left(\mathcal{T}g\right)}_k=\left\{{\left(\mathcal{T}\right)}_{i\mathrm{,}k}\mathrm{,}i\in {\mathcal{T}}_k\right\}\mathrm{,}$

$\mathcal{T}g=\left\{{\left(\mathcal{T}g\right)}_k\mathrm{,}k\in {\mathcal{J}}^M\right\}\mathrm{.}$

Now,

${\Vert \mathcal{T}g\Vert }_{\mathcal{X}}^2={\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\frac{1}{\sigma }{\left|{\mathcal{S}}^k{\left(2\mu \partial \left(q^k\right)-g^k\right)}_i\right|}^2\mathrm{.}$ (23)

We use the facts ${\displaystyle {\sum }_{i\in {\mathcal{I}}_k}}{\left({\mathcal{S}}^k{g}^k\right)}_i^2={\displaystyle {\sum }_{i\in {\mathcal{I}}_k}}{\left(g^k\right)}_i^2$ and ${\displaystyle {\sum }_{i\in {\mathcal{I}}_k}}{\left({\mathcal{S}}^k{q}^k\right)}_i{\left({\mathcal{S}}^k{g}^k\right)}_i={\displaystyle {\sum }_{i\in {\mathcal{I}}_k}}{q}_i^k{g}_i^k$ and show

${\Vert \mathcal{T}g\Vert }_{\mathcal{X}}^2={\Vert g\Vert }_{\mathcal{X}}^2-4{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(g_i^k-\sigma \partial q_i^k\right)\partial q_i^k\mathrm{.}$ (24)

We now formulate a relaxed version of a fixed point iteration: for $ϵ\in \left[\mathrm{0,1}\right)$

$g^{l+1}=\left(1-ϵ\right)\mathcal{T}\left(g^l\right)+ϵg^l\mathrm{.}$ (25)

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

#Math_108# (26)

Thus, we introduce $e^{l+1}:=q^{l+1}-q$ . Then $e^{l+1}$ solves a non-linear differential equation with nonlinearity $g_i\left(e_i^{l+1}+q_i\right)-g_i\left(q_i\right)$ , 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 $q^{l+1}$ . We consider the following integration by parts formula after multiplying by a test function $\varphi$ .

$0={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{e}_i^{l+1}-{\beta }_i{\partial }_{xx}{e}_i^{l+1}+g_i\left(e_i^{l+1}+q_i\right)-g_i\mathrm{<mo\; stretchy="false">\; (\; </mo>}{q}_i\mathrm{<mo\; stretchy="false">\; )\; </mo>}\right){\varphi }_i\text{d}x$ (27)

$=-{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\beta }_i{\partial }_x{e}_i^{l+1}\left(n_k\right){\varphi }_i\left(n_k\right)$ (28)

#Math_116# (29)

We now take the test function to be equal to $e_i^{l+1}$ and obtain:

$\begin{array}{l}{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\partial }_x{e}_i^{l+1}\left(n_k\right)e_i^{l+1}\left(n_k\right)\\ ={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{e}_i^{l+1}{e}_i^{l+1}+{\beta }_i{\partial }_x{e}_i^{l+1}{\partial }_x{e}_i^{l+1}+\left(g_i\left(e_i^{l+1}+q_i\right)-g_i\left(q_i\right)\right)e_i^{l+1}\right)\text{d}x.\end{array}$ (30)

We use the boundary condition at the interfaces in the form

$d_{ik}{e}_i^{l+1}\left(n_k\right)=g_{ik}^{l+1}-\sigma {\partial }_x{e}_i^{l+1}\left(n_k\right)\mathrm{,}k\in {\mathcal{J}}^M\mathrm{.}$

This identity is used in the identity (24), evaluated for the error:

${\Vert \mathcal{T}g\Vert }_{\mathcal{X}}^2={\Vert g\Vert }_{\mathcal{X}}^2-4{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(g_i^k-\sigma \partial e_i^k\right)\partial e_i^k\mathrm{.}$ (31)

We obtain

#Math_121# (32)

We assume

$\left(g_i\left(x\mathrm{<mo>\; ;\; </mo>}s\right)-g_i\left(x\mathrm{<mo>\; ;\; </mo>}t\right)\right)\left(s-t\right)\ge \mathrm{0,}\forall x\in \left(\mathrm{0,}{\mathcal{l}}_i\right)\mathrm{,}i\in \mathcal{I}$ (33)

and define the bilinear form

$a_i\left({\psi }_i\mathrm{,}\varphi \right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{\psi }_i{\varphi }_i+{\beta }_i{\partial }_x{\psi }_i{\partial }_x{\varphi }_i\right)\text{d}x\mathrm{.}$ (34)

We define the corresponding quadratic form applied to $e^l$

$a_i\left(e_i^l\mathrm{,}{e}_i^l\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{e}_i^l{e}_i^l+{\beta }_i{\partial }_x{e}_i^l{\partial }_x{e}_i^l\right)\text{d}x\mathrm{,}$ (35)

which is certainly bounded below by ${\Vert e_i\Vert }_{L^2\left(\mathrm{0,}{\mathcal{l}}_i\right)}$ . Then the error iteration is

${\Vert g^{l+1}\Vert }_{\mathcal{X}}^2={\Vert \mathcal{T}{g}^l\Vert }_{\mathcal{X}}^2={\Vert g^l\Vert }_{\mathcal{X}}^2-4{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{a}_i\left(e_i^l\mathrm{,}{e}_i^l\right)$ (36)

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 $ϵ\in \left[\mathrm{0,1}\right]$

${\Vert g^{l+1}\Vert }_{\mathcal{X}}^2\le {\Vert g^l\Vert }_{\mathcal{X}}^2-4\left(1-ϵ\right){\displaystyle \underset{i\in \mathcal{I}}{\sum }}{a}_i\left(e_i^l\mathrm{,}{e}_i^l\right)\mathrm{.}$ (37)

We iterate in (36) or (37) down from $l$ to zero. Then we obtain

$\left\{g^l\right\}\text{is}\text{bounded}\mathrm{,}{a}_i\left(e_i^l\mathrm{,}{e}_i^l\right)\to \mathrm{0,}l\to \infty \mathrm{.}$ (38)

Clearly, for ${\alpha }_i>0$ , this shows that the $H^1\left(\mathrm{0,}{\mathcal{l}}_i\right)$ -error strongly tends to zero. Then also the traces tend to zero as $l\to \infty$ and, therefore, the iteration converges.

Theorem 2. Under assumption (33), for each $ϵ\in \left[\mathrm{0,1}\right)$ the iteration (25) with (26) and (21), (22) converges as $l\to \infty$ . 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 $dx=1/1000$ , ${\alpha }_i=10$ , $f_1=1$ , $f_2=0.1$ , $f_3=0.5$ . The nonlinearity is weighted by a factor $\gamma =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 $tol=1e-10$ . 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 $L^{\infty }$ -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:

$\begin{array}{l}\underset{\left(q\mathrm{,}u\right)}{min}I\left(q\mathrm{,}u\right)\mathrm{<mo>\; :\; </mo>}=\frac{\kappa }{2}{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\Vert q_i-q_i^0\Vert }^2+\frac{\nu }{2}{\displaystyle \underset{k\in {\mathcal{J}}_N^S}{\sum }}{\left|u_k\right|}^2\\ \text{subject}\text{to}\\ {\alpha }_i{q}_i\left(x\right)-{\beta }_i{\partial }_{xx}{q}_i\left(x\right)+g_i\left(x\mathrm{<mo>\; ;\; </mo>}{q}_i\left(x\right)\right)=f_i\left(x\right)\mathrm{,}x\in \left(\mathrm{0,}{\mathcal{l}}_i\right)\mathrm{,}i\in \mathcal{I}\\ q_i\left(n_k\right)=\mathrm{0,}i\in {\mathcal{I}}_k\mathrm{,}{n}_k\in {\mathcal{J}}_D^S\\ {\partial }_x{q}_i\left(n_k\right)=u_k\mathrm{,}i\in {\mathcal{I}}_k\mathrm{,}{n}_k\in {\mathcal{J}}_N^S\\ {\beta }_i{\partial }_x{q}_i\left(n_k\right)={\beta }_j{\partial }_x{q}_j\left(n_k\right)\mathrm{,}i\ne j\in {\mathcal{I}}_k\mathrm{,}k\in {\mathcal{J}}^M\mathrm{,}\\ {\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{q}_i\left(n_k\right)=\mathrm{0,}k\in {\mathcal{J}}^M\mathrm{.}\end{array}$ (39)

The corresponding optimality system then reads as follows:

$\begin{array}{l}{\alpha }_i{q}_i-{\beta }_i{\partial }_{xx}{q}_i+g_i\left(q_i\right)=f_i\text{in}\left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ {\alpha }_i{\rho }_i-{\beta }_i{\partial }_{xx}{\rho }_i+{g^{\prime }}_i\left(q_i\right){\rho }_i=-\kappa \left(q_i-q_i^0\right)\text{in}\left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ q_i\left(n_k\right)=0,{\rho }_i\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_D^S\\ {\partial }_x{q}_i\left(n_k\right)=\frac{1}{\nu }{\rho }_i\left(n_k\right),{\partial }_x{\rho }_i\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_N^S\\ {\beta }_i{\partial }_x{q}_i\left(n_k\right)={\beta }_j{\partial }_x{q}_j\left(n_k\right),{\beta }_i{\partial }_x{\rho }_i\left(n_k\right)={\beta }_j{\partial }_x{\rho }_j\left(n_k\right),i\ne j\in {\mathcal{I}}_k,k\in {\mathcal{J}}^M\\ {\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{q}_i\left(n_k\right)=0={\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\rho }_i\left(n_k\right),k\in {\mathcal{J}}^M.\end{array}$ (40)

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:

$\begin{array}{l}{\alpha }_i{q}_i^{l+1}-{\beta }_i{\partial }_{xx}{q}_i^{l+1}=f_i\text{in}\left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ {\alpha }_i{\rho }_i^{l+1}-{\beta }_i{\partial }_{xx}{\rho }_i^{l+1}=-\kappa \left(q_i^{l+1}-q_i^0\right)\text{in}\left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ q_i^{l+1}\left(n_k\right)=0,{\rho }_i^{l+1}\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_D^S\\ {\partial }_x{q}_i^{l+1}\left(n_k\right)=\frac{1}{\nu }{\rho }_i^{l+1}\left(n_k\right),{\partial }_x{\rho }_i^{l+1}\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_N^S\\ d_{ik}{q}_i^{l+1}\left(n_k\right)+{\lambda }_k{\partial }_x{q}_i^{l+1}\left(n_k\right)+{\mu }_k{\partial }_x{\rho }_i^{l+1}\left(n_k\right)\\ ={\lambda }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{q}_k^l\left(n_k\right)-{\beta }_i{\partial }_x{q}_i^l\left(n_k\right)\right)+{\mu }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{\rho }_k^l\left(n_k\right)-{\beta }_i{\partial }_x{\rho }_i^l\left(n_k\right)\right)\\ -\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{q}_j^l\left(n_k\right)-d_{ik}{q}_i^l\left(n_k\right)\right)=g_{ik}^{l+1},\\ d_{ik}{\rho }_i^{l+1}\left(n_k\right)+{\lambda }_k{\partial }_x{\rho }_i^{l+1}\left(n_k\right)-{\mu }_k{\partial }_x{q}_i^{l+1}\left(n_k\right)\\ ={\lambda }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{\rho }_k^l\left(n_k\right)-{\beta }_i{\partial }_x{\rho }_i^l\left(n_k\right)\right)-{\mu }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{q}_j^l\left(n_k\right)-{\beta }_i{\partial }_x{q}_i^l\left(n_k\right)\right)\\ -\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\rho }_k^l\left(n_k\right)-d_{ik}{\rho }_i^l\left(n_k\right)\right)=h_{ik}^{l+1},k\in {\mathcal{J}}^M.\end{array}$ (41)

The same arguments that led from (16), (15) to (17), (18) apply to show that, upon convergence as $l\to \infty$ , the system (41) tends to (40). Now, (41) decomposes the fully connected problem (40) to a problem on a single edge $i\in \mathcal{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. To this end, it is sufficient to consider three cases: a.) the edge $i$ connects a controlled Neumann $j\in {\mathcal{J}}_N^S$ node with a multiple node $k\in {\mathcal{J}}^M$ at which the domain decomposition is active, b.) the edge $i$ connects a controlled Neumann node $j\in {\mathcal{J}}_N^S$ with multiple node $k\in {\mathcal{J}}^M$ at which the domain decomposition is active, c.) the edge $i$ connects two multiple nodes $j\mathrm{,}k\in {\mathcal{J}}^M$ .

Case a.):

#Math_162# (42)

Case b.):

#Math_163# (43)

Case c.):

#Math_164# (44)

Remark 4.1.

・ 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 ${\left({\rho }_i\right)}_{i\in \mathcal{I}}$ . This provides a guess for the controls ${\left(u_j\right)}_{j\in {\mathcal{J}}^S}$ . Therefore, we can establish the states ${\left(q_i\right)}_{i\in \mathcal{I}}$ . The states ${\left(q_i\right)}_{i\in \mathcal{I}}$ , in turn, are inserted into the adjoint problem and that system is then solved for ${\left({\rho }_i\right)}_{i\in \mathcal{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 ${\left(q_i\right)}_{i\in \mathcal{I}}$ and the adjoint system for ${\left({\rho }_i\right)}_{i\in \mathcal{I}}$ , individually. For given ${\left(q_i\right)}_{i\in \mathcal{I}}$ , 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 ${\left(q_i\right)}_{i\in \mathcal{I}}$ and ${\left({\rho }_i\right)}_{i\in \mathcal{I}}$ . 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 ${\rho }_i$ for each $i\in \mathcal{I}$ then solve for $q_i$ which is introduced then in the local adjoint equations. This is then followed by the solve for ${\rho }_i$ and the update of the boundary data $g_{ik}\mathrm{,}{h}_{ik}$ 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. Wellposedness

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 $\mathcal{H}$ :

$\mathcal{H}={\displaystyle {\prod }_{i\in \mathcal{I}}{H}^1\left(\mathrm{0,}{\mathcal{l}}_i\right)}\mathrm{.}$

Indeed, we also define the energy space

$\mathcal{V}\mathrm{<mo>\; :\; </mo>}=\left\{q\in \mathcal{H}|q_i\left(n_k\right)=\mathrm{0,}i\in {\mathcal{I}}_k\mathrm{,}{n}_k\in {\mathcal{J}}_D^S\mathrm{,}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{q}_i\left(n_k\right)=\mathrm{0,}k\in {\mathcal{J}}^M\right\}$

and the operator $\mathcal{A}$ as follows:

$\mathcal{A}q\mathrm{<mo>\; :\; </mo>}={\alpha }_i{q}_i\left(x\right)-{\beta }_i{\partial }_{xx}{q}_i\left(x\right)\text{in}\mathcal{H}$

$\begin{array}{l}D\left(\mathcal{A}\right)\mathrm{<mo>\; :\; </mo>}=\{q\in \mathcal{V}|{\beta }_i{\partial }_x{q}_i\left(n_k\right)={\beta }_j{\partial }_x{q}_j\left(n_k\right)\mathrm{,}i\ne j\in {\mathcal{I}}_k\mathrm{,}k\in {\mathcal{J}}^M\mathrm{,}\\ {\partial }_x{q}_i\left(n_k\right)=\mathrm{0,}i\in {\mathcal{I}}_k\mathrm{,}{n}_k\in {\mathcal{J}}_N^S\}\mathrm{.}\end{array}$ (45)

It is a matter of applying standard integration by parts to show that, indeed, $\mathcal{A}$ is symmetric and positive definite in $\mathcal{H}$ such that $\mathcal{A}$ can be extended as a self-adjoint operator in $\mathcal{H}$ . Then it is standard to show that $\mathcal{A}$ can be extended to a bounded coercive map from $\mathcal{V}$ into its dual ${\mathcal{V}}^{\mathrm{<mo>\; *\; </mo>}}$ . If we assume (33) and define the Nemitskji operator $\mathcal{G}\left(q\right)\left(x\right)\mathrm{<mo>\; :\; </mo>}=g\left(q\left(x\right)\right)$ then $\mathcal{G}$ is strictly monotone and continuous. Hence, according to [11] , $\mathcal{A}+\mathcal{G}$ is strictly monotone and continuous and, hence, the semi-linear problem admits a unique solution $q\in \mathcal{V}$ . Clearly, for regular right hand sides $f$ , the solution is in $D\left(\mathcal{A}\right)$ .

5.2. Smoothness of the Control-to-State-Map

Let ${\hat{q}}_t\left(\hat{u}\right)$ be the solution of (12) with $u$ replaced with $u+t\hat{u}$ and let $q$ be the solution of (12). We denote by $e:=\hat{q}-q$ the difference of these solutions. We obtain

#Math_204# (46)

Dividing by $t$ and letting $t$ tend to zero in (12) implies with $e^{\prime }:=e^{\prime }\left(u\right)\left(\hat{u}\right)=\delta e\left(u\right)\left(\hat{u}\right)$

$\begin{array}{l}{\alpha }_i{e^{\prime }}_i\left(x\right)-{\beta }_i{\partial }_{xx}{e^{\prime }}_i\left(x\right)+{g^{\prime }}_i\left(q_i\right)\left(x\right){e^{\prime }}_i\left(x\right)=0,x\in \left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ {\beta }_i{\partial }_x{e^{\prime }}_i\left(n_k\right)={\beta }_j{\partial }_x{e^{\prime }}_j\left(n_k\right),i\ne j\in {\mathcal{I}}_k,k\in {\mathcal{J}}^M\\ {e^{\prime }}_i\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_D^S\\ {\partial }_x{e^{\prime }}_i\left(n_k\right)={\hat{u}}_k,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_N^S\\ {\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e^{\prime }}_i\left(n_k\right)=0,k\in {\mathcal{J}}^M.\end{array}$ (47)

For the solution $q$ of (12), applying the standard Lax-Milgram Lemma, $e^{\prime }$ 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\in {\mathcal{V}}^{\mathrm{<mo>\; *\; </mo>}}$ , there exists a unique solution $q\in \mathcal{V}$ 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).

5.3. A Priori Error Estimates for the Optimality System

We denote the errors $e_i^l:=q_i^l-q_i$ and $p_i^l:={\rho }_i^l-{\rho }_i$ for $i\in \mathcal{I}$ and $l=0,1,\cdots$ . These errors solve the system equations:

$\begin{array}{l}{\alpha }_i{e}_i^{l+1}-{\beta }_i{\partial }_{xx}{e}_i^{l+1}+g_i\left(e_i^{l+1}+q_i\right)-g_i\left(q_i\right)=0,\text{in}\left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ {\alpha }_i{p}_i^{l+1}-{\beta }_i{\partial }_{xx}{p}_i^{l+1}+{g^{\prime }}_i\left(e_i^{l+1}+q_i\right)p_i^{l+1}+\left({g^{\prime }}_i\left(e_i^{l+1}+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i\\ =-\kappa e_i^{l+1}\text{in}\left(0,{\mathcal{l}}_i\right),i\in \mathcal{I}\\ e_i^{l+1}\left(n_k\right)=0,p_i^{l+1}\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_D^S\\ {\partial }_x{e}_i^{l+1}\left(n_k\right)=\frac{1}{\nu }{p}_i^{l+1}\left(n_k\right),{\partial }_x{p}_i^{l+1}\left(n_k\right)=0,i\in {\mathcal{I}}_k,n_k\in {\mathcal{J}}_N^S\\ d_{ik}{e}_i^{l+1}\left(n_k\right)+{\lambda }_k{\partial }_x{e}_i^{l+1}\left(n_k\right)+{\mu }_k{\partial }_x{p}_i^{l+1}\left(n_k\right)\\ ={\lambda }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{e}_k^l\left(n_k\right)-{\beta }_i{\partial }_x{e}_i^l\left(n_k\right)\right)+{\mu }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{p}_k^l\left(n_k\right)-{\beta }_i{\partial }_x{p}_i^l\left(n_k\right)\right)\\ -\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e}_j^l\left(n_k\right)-d_{ik}{e}_i^l\left(n_k\right)\right)=g_{ik}^{l+1},\\ d_{ik}{p}_i^{l+1}\left(n_k\right)+{\lambda }_k{\partial }_x{p}_i^{l+1}\left(n_k\right)-{\mu }_k{\partial }_x{e}_i^{l+1}\left(n_k\right)\\ ={\lambda }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{p}_k^l\left(n_k\right)-{\beta }_i{\partial }_x{p}_i^l\left(n_k\right)\right)-{\mu }_k\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{\beta }_j{\partial }_x{e}_j^l\left(n_k\right)-{\beta }_i{\partial }_x{e}_i^l\left(n_k\right)\right)\\ -\left(\frac{2}{d_k}{\displaystyle \underset{j\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{p}_k^l\left(n_k\right)-d_{ik}{p}_i^l\left(n_k\right)\right)=h_{ik}^{l+1},k\in {\mathcal{J}}^M.\end{array}$ (48)

We prove the following

Lemma 5. The solutions $e_i\mathrm{,}{p}_i$ for $i\in \mathcal{I}$ of (48) satisfies the estimate

${\Vert e_i\Vert }_{H^1\left(\mathrm{0,}{\mathcal{l}}_i\right)}^2+{\Vert p_i\Vert }_{H^1\left(\mathrm{0,}{\mathcal{l}}_i\right)}^2+\gamma \left(e_i{\left(0\right)}^2+e_i{\left({\mathcal{l}}_i\right)}^2+p_i{\left(0\right)}^2+p_i{\left({\mathcal{l}}_i\right)}^2\right)\le C\left\{g_{ij}^2+h_{ij}^2\right\}\mathrm{.}$ (49)

More precisely, for ${\lambda }_k=\lambda ,\forall k\in {\mathcal{J}}^M$ , we obtain

$\left(e_i{\left({\mathcal{l}}_i\right)}^2+p_i{\left({\mathcal{l}}_i\right)}^2\right)\le \frac{4\nu }{\lambda }{\displaystyle \underset{i=1}{\overset{2}{\sum }}}\left\{g_{ij}^2+h_{ij}^2\right\}.$ (50)

Remark 5.1. As a result, for small data $g_{ij}\mathrm{,}{h}_{ij}$ , we have small solutions.

Proof of Lemma 5: We multiply the equations in (48) by $e_i^{l+1}$ and $p_i^{l+1}$ , respectively, integrate and then use integration by parts. For the sake of brevity, we leave the full arguments to the reader.

5.4. Convergence

$\left(g,h\right)\in \mathcal{X}:={\displaystyle {\prod }_{k\in {\mathcal{J}}^M}{\displaystyle {\prod }_{i\in {\mathcal{I}}_k}{ℝ}^2}},{\Vert \left(g,h\right)\Vert }_{\mathcal{X}}^2:={\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left({\left|g_{i,k}\right|}^2+{\left|h_{i,k}\right|}^2\right),$ (51)

$\mathcal{T}\mathrm{<mo>\; :\; </mo>}\mathcal{X}\to \mathcal{X}\mathrm{,}$ (52)

$\begin{array}{l}\mathcal{T}{\left(g\mathrm{,}h\right)}_{i\mathrm{,}k}\mathrm{<mo>\; :\; </mo>}=({\mathcal{S}}^k{\left(2\left({\lambda }_k{\partial }_x\left(e^k\right)+{\mu }_k{\partial }_x\left(p^k\right)\right)-g^k\right)}_i\mathrm{,}\\ {\mathcal{S}}^k{\left(2\left({\lambda }_k{\partial }_x\left(p^k\right)-{\mu }_k{\partial }_x\left(e^k\right)\right)-h^k\right)}_i),k\in {\mathcal{J}}^M\mathrm{,}i\in {\mathcal{I}}_k\mathrm{,}\end{array}$

$\mathcal{T}{\left(g\mathrm{,}h\right)}_k=\left\{\left(\mathcal{T}\right){\left(g\mathrm{,}h\right)}_{i\mathrm{,}k}\mathrm{,}i\in {\mathcal{I}}_k\right\}\mathrm{,}$

$\mathcal{T}\left(g\mathrm{,}h\right)=\left\{\mathcal{T}{\left(g\mathrm{,}h\right)}_k\mathrm{,}k\in {\mathcal{J}}^M\right\}\mathrm{.}$

Now,

#Math_233# (53)

We multiply the state equation for the errors $e_i$ , $p_i$ by $e_i$ and $p_i$ , respectively.

$\begin{array}{c}0={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{e}_i-{\beta }_i{\partial }_{xx}{e}_i+g_i\left(e_i+q_i\right)-g_i\left(q_i\right)\right)e_i\text{d}x\\ =-{\displaystyle \underset{k\in \mathcal{I}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\beta }_i{\partial }_x{e}_i\left(n_k\right)e_i\left(n_k\right)\\ +{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left\{{\alpha }_i{\left(e_i\right)}^2+{\beta }_i{\left({\partial }_x{e}_i\right)}^2+\left(g_i\left(e_i+q_i\right)-g_i\left(q_i\right)\right)e_i\right\}\text{d}x,\end{array}$ (54)

$\begin{array}{c}0={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{p}_i-{\beta }_i{\partial }_{xx}{p}_i+{g^{\prime }}_i\left(e_i+q_i\right)p_i+\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i+{\kappa }_i{e}_i\right)p_i\text{d}x\\ =-{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\beta }_i{\partial }_x{p}_i\left(n_k\right)p_i\left(n_k\right)+{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\{{\alpha }_i{\left(p_i\right)}^2+{\beta }_i{\left({\partial }_x{p}_i\right)}^2\\ +{g^{\prime }}_i\left(e_i+q_i\right)p_i^2+\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i{p}_i+{\kappa }_i{e}_i{p}_i\}\text{d}x.\end{array}$ (55)

Now, we reverse the roles and obtain

$\begin{array}{c}0={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{e}_i-{\beta }_i{\partial }_{xx}{e}_i+g_i\left(e_i+q_i\right)-g_i\left(q_i\right)\right)p_i\text{d}x\\ =-{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\beta }_i{\partial }_x{e}_i\left(n_k\right)p_i\left(n_k\right)\\ +{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left\{{\alpha }_i{e}_i{p}_i+{\beta }_i{\partial }_x{e}_i{\partial }_x{p}_i+\left(g_i\left(e_i+q_i\right)-g_i\left(q_i\right)\right)e_i{p}_i\right\}\text{d}x,\end{array}$ (56)

$\begin{array}{c}0={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{p}_i-{\beta }_i{\partial }_{xx}{p}_i+{g^{\prime }}_i\left(e_i+q_i\right)p_i+\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_{i}d+{\kappa }_i{e}_i\right)e_i\text{d}x\\ =-{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{\beta }_i{\partial }_x{p}_i\left(n_k\right)e_i\left(n_k\right)+{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\{{\alpha }_i{p}_i{e}_i+{\beta }_i{\partial }_x{p}_i{\partial }_x{e}_i\\ +{g^{\prime }}_i\left(e_i+q_i\right)p_i{e}_i+\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i{e}_i+{\kappa }_i{e}_i^2\}\text{d}x.\end{array}$ (57)

From now on we assume that all ${\lambda }_k=:\lambda ,{\mu }_k:=\mu$ 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

$\begin{array}{l}{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{e}_i\left(n_k\right)+\mu {\beta }_i{\partial }_x{p}_i\left(n_k\right)\right)\\ ={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\lambda \left\{{\alpha }_i{\left(e_i\right)}^2+{\beta }_i{\left({\partial }_x{e}_i\right)}^2+\left(g_i\left(e_i+q_i\right)-g_i\left(q_i\right)\right)e_i\right\}\\ +\mu \left\{{\alpha }_i{e}_i{p}_i+{\beta }_i{\partial }_x{e}_i{\partial }_x{p}_i+{g^{\prime }}_i\left(e_i+q_i\right)e_i{p}_i+\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i{e}_i+{\kappa }_i{e}_i^2\right\}\text{d}x\end{array}$

$\begin{array}{l}{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{p}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{p}_i\left(n_k\right)-\mu {\beta }_i{\partial }_x{e}_i\left(n_k\right)\right)\\ ={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\lambda \left\{{\alpha }_i{\left(p_i\right)}^2+{\beta }_i{\left({\partial }_x{p}_i\right)}^2+{g^{\prime }}_i\left(e_i+q_i\right)p_i^2+\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i{p}_i+{\kappa }_i{e}_i{p}_i\right\}\\ -\mu \left\{{\alpha }_i{e}_i{p}_i+{\beta }_i{\partial }_x{e}_i{\partial }_x{p}_i+g_i\left(e_i+q_i\right)e_i{p}_i\right\}\text{d}x\end{array}$

We add the latter equations and obtain

$\begin{array}{l}{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{e}_i\left(n_k\right)+\mu {\beta }_i{\partial }_x{p}_i\left(n_k\right)\right)\\ +{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{p}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{p}_i\left(n_k\right)-\mu {\beta }_i{\partial }_x{e}_i\left(n_k\right)\right)\\ ={\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\lambda \left\{{\alpha }_i{e}_i^2+{\alpha }_i{p}_i^2+{\beta }_i{\left({\partial }_x{e}_i\right)}^2+{\beta }_i{\left({\partial }_x{p}_i\right)}^2\right\}\text{d}x+{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{\kappa }_i\left(\mu e_i^2+\lambda p_i{e}_i\right)\text{d}x\\ +\frac{1}{\nu }{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(\mu {\left(p_i\left(n_k\right)\right)}^2-\lambda d_{ik}{e}_i\left(n_k\right)p_i\left(n_k\right)\right)\\ +\lambda {\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left\{\left(g_i\left(e_i+q_i\right)-g_i\left(q_i\right)\right)e_i+{g^{\prime }}_i\left(e_i+q_i\right)p_i^2+\left({g^{\prime }}_i\left(e_i+y_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i{p}_i\right\}\text{d}x\\ +\mu {\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left\{{g^{\prime }}_i\left(e_i+q_i\right)e_i{p}_i+\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i{e}_i-\left(g_i\left(e_i+q_i\right)-g_i\left(q_i\right)\right)e_i{p}_i\right\}\text{d}x\\ \ge {\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\lambda \left\{{\alpha }_i\left(e_i^2+p_i^2\right)+{\beta }_i{\left({\partial }_x{e}_i\right)}^2+{\beta }_i{\left({\partial }_x{p}_i\right)}^2\right\}\text{d}x+{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{\kappa }_i\left(\mu e_i^2+\lambda p_i{e}_i\right)\text{d}x\\ +\frac{1}{\nu }{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(\mu {\left(p_i\left(n_k\right)\right)}^2-\lambda d_{ik}{e}_i\left(n_k\right)p_i\left(n_k\right)\right)\\ +{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i\right)\right){\rho }_i\left(\lambda p_i+\mu e_i\right)+\mu \left({g^{\prime }}_i\left(e_i+q_i\right)-{g^{\prime }}_i\left(q_i+\theta e_i\right)\right)e_i{p}_i\text{d}x\\ =I+II+III.\end{array}$ (58)

We are going to estimate the third integral. For that matter we assume that ${g^{\prime }}_i\left(s\right)$ admits a Lipschitz constant $L_i$ .

$\begin{array}{c}III\le {\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{L}_i\left|{\rho }_i\right|\left|e_i\right|\left(\lambda \left|p_i\right|+\mu \left|e_i\right|\right)+\mu L_i\left(1-\theta \right)\left|p_i\right|{\left(e_i\right)}^2\text{d}x\\ \le {\displaystyle \underset{i\in \mathcal{I}}{\sum }}{L}_i\left\{\left(\left(\mu +\frac{\lambda }{2}\right){\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}+\mu {\Vert p_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right){\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{e}_i^2\text{d}x+\frac{\lambda }{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{p}_i^2\text{d}x\right\}\end{array}$ (59)

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

$\begin{array}{l}{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{\kappa }_i\left(\mu e_i^2+\lambda p_i{e}_i\right)\text{d}x\ge {\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{\kappa }_i\left(\left(\mu -\frac{\lambda }{\delta }\right)e_i^2-\lambda \delta p_i^2\right)\text{d}x\\ \frac{1}{\nu }{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(\mu {\left(p_i\left(n_k\right)\right)}^2-\lambda d_{ik}{e}_i\left(n_k\right)p_i\left(n_k\right)\right)\\ \ge \frac{1}{\nu }{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(\left(\mu -\frac{\lambda }{\delta }\right){\left(p_i\left(n_k\right)\right)}^2-\lambda \delta e_i{\left(n_k\right)}^2\right)\end{array}$ (60)

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

$\begin{array}{l}{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{e}_i\left(n_k\right)+\mu {\beta }_i{\partial }_x{p}_i\left(n_k\right)\right)\\ +{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{p}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{p}_i\left(n_k\right)-\mu {\beta }_i{\partial }_x{e}_i\left(n_k\right)\right)\\ \ge {\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\lambda \left\{{\alpha }_i\left(e_i^2+p_i^2\right)+{\beta }_i{\left({\partial }_x{e}_i\right)}^2+{\beta }_i{\left({\partial }_x{p}_i\right)}^2\right\}\text{d}x\\ +{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{\kappa }_i\left(\left(\mu -\frac{\lambda }{\delta }\right)e_i^2-\lambda \delta p_i^2\right)\text{d}x+\frac{1}{\nu }{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(\left(\mu -\frac{\lambda }{\delta }\right){\left(p_i\left(n_k\right)\right)}^2-\lambda \delta e_i{\left(n_k\right)}^2\right)\\ -{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{L}_i\left\{\left(\left(\mu +\frac{\lambda }{2}\right){\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}+\mu {\Vert p_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right){\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{e}_i^2\text{d}x+\frac{\lambda }{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{p}_i^2\text{d}x\right\}\end{array}$ (61)

We now group the corresponding quadratic expressions.

$\begin{array}{l}{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{e}_i\left(n_k\right)+\mu {\beta }_i{\partial }_x{p}_i\left(n_k\right)\right)\\ +{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{p}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{p}_i\left(n_k\right)-\mu {\beta }_i{\partial }_x{e}_i\left(n_k\right)\right)\\ \ge {\displaystyle \underset{i\in \mathcal{I}}{\sum }}(\lambda {\alpha }_i+{\kappa }_i\left(\mu -\frac{\lambda }{\delta }\right)-L_i\left(\left(\mu +\frac{\lambda }{2}\right){\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}+\mu {\Vert p_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right)-\lambda \delta c_i^2\left(n\right)){\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{e}_i^2\text{d}x\\ +\left(\lambda {\alpha }_i-\lambda \delta -L_i\frac{\lambda }{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right){\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{p}_i^2\text{d}x+{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left\{\lambda \left(1-\delta c^1\left(n\right)\right){\left({\partial }_x{e}_i\right)}^2+\lambda {\left({\partial }_x{p}_i\right)}^2\right\}\text{d}x\\ +\frac{1}{\nu }{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}\left(\left(\mu -\frac{\lambda }{\delta }\right){\left(p_i\left(n_k\right)\right)}^2\right),\end{array}$ (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 $\mu -\frac{\lambda }{\delta }\le 0$ , we need to absorb the corresponding boundary term again using the estimate [13] . It is obvious that $b_i$ can be positive iff $\left({\alpha }_i-\delta -L_i\frac{1}{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(\mathrm{0,}{\mathcal{l}}_i\right)}\right)$ and $\lambda$ are positive. The only

parameters that we can select in order achieve positive, respectively, non- negative coefficient are the two parameters ${\kappa }_i\mathrm{,}{\nu }_i$ coming from the cost function and the parameters $\lambda \ge \mathrm{0,}\mu \ge 0$ provided for the algorithm. Moreover, the coefficient ${\alpha }_i>0$ becomes relevant. We recall the meaning of ${\alpha }_i$ : it is $\frac{1}{\Delta t}$ ! So it becomes obvious from (63) that the norm of the reference solution to the adjoint equation ${\Vert {\rho }_i\Vert }_{\infty }$ and the Lipschitz-constant $L_i$ , reflecting the stiffness of the nonlinear term come into play. We thus need small $\Delta t>0$ to compensate the remaining terms. The question to be discussed below then is as to whether the maximum-norm of the solution ${\rho }_i$ of the adjoint equation which, in turn, involves ${\alpha }_i$ is small against ${\alpha }_i$ . Only in this case, we can choose $\lambda >0$ in order to have $b_i>0$ . If, on the other side, $\lambda =0$ , we have to compensate $L_i\Vert p_i\Vert$ , in this case the adjoint error, by choosing ${\kappa }_i$ sufficiently large and $\mu >0$ in order to have $a_i>0$ . The appearance of the adjoint error, in case $\mu >0$ , necessitates an a posteriori error estimate. We discuss the following cases $\lambda =0,\mu >0$ , $\lambda >0,\mu =0$ and $\lambda >0,\mu >0$ :

[Case 1.] $\lambda =0,\mu >0$ :

$a_i=\mu \left({\kappa }_i-L_i\left({\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}+{\Vert p_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right)\right),b_i:=0$ (64)

In this case

$\begin{array}{l}{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{e}_i\left(n_k\right)+\mu {\beta }_i{\partial }_x{p}_i\left(n_k\right)\right).\\ +{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{p}_i\left(n_k\right){\beta }_i\left(\lambda {\partial }_x{p}_i\left(n_k\right)-\mu {\partial }_x{e}_i\left(n_k\right)\right)\\ \ge {\displaystyle \underset{i\in \mathcal{I}}{\sum }}\mu \left({\kappa }_i-L_i\left({\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}+{\Vert p_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right)\right){\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{e}_i^2\text{d}x+\frac{\mu }{\nu }{\displaystyle \underset{k\in {\mathcal{J}}^M}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{\left(p_i\left(n_k\right)\right)}^2.\end{array}$ (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 ${\kappa }_i$ has to be sufficiently large, while ${\alpha }_i$ plays no role.

[Case 2.] $\lambda >0,\mu =0$ :

$\begin{array}{l}{a}_i:=\lambda \left({\alpha }_i-{\kappa }_i\frac{1}{\delta }-L_i\left(\frac{1}{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right)-\delta c_i^2\left(n\right)\right)\\ b_i:=\lambda \left({\alpha }_i-\delta -L_i\frac{1}{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(0,{\mathcal{l}}_i\right)}\right).\end{array}$ (66)

In this case we obtain

${\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{e}_i\left(n_k\right)\left(\lambda {\beta }_i{\partial }_x{e}_i\left(n_k\right)+\mu {\beta }_i{\partial }_x{p}_i\left(n_k\right)\right)$

$+{\displaystyle \underset{k\in \mathcal{J}}{\sum }}{\displaystyle \underset{i\in {\mathcal{I}}_k}{\sum }}{d}_{ik}{p}_i\left(n_k\right){\beta }_i\left(\lambda {\partial }_x{p}_i\left(n_k\right)-\mu {\partial }_x{e}_i\left(n_k\right)\right)$ (67)

$\ge \lambda {\displaystyle \underset{i\in \mathcal{I}}{\sum }}\{\left({\alpha }_i-{\kappa }_i\frac{1}{\delta }-L_i\frac{1}{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(\mathrm{0,}{\mathcal{l}}_i\right)}-\delta c_i^2\left(n\right)\right){\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{e}_i^2\text{d}x$

$+\left({\alpha }_i-\delta -L_i\frac{1}{2}{\Vert {\rho }_i\Vert }_{L^{\infty }\left(\mathrm{0,}{\mathcal{l}}_i\right)}-\delta c_i^2\left(n\right)\right){\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}{p}_i^2\text{d}x$ (68)

$+{\displaystyle \underset{i\in \mathcal{I}}{\sum }}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left(1-\delta c^1\left(n\right)\right){\left({\partial }_x{e}_i\right)}^2+\left(1-\delta c^1\left(n\right)\right){\left({\partial }_x{p}_i\right)}^2\text{d}x\}$

In this case ${\alpha }_i$ has to absorb all negative terms and, in fact, the penalty parameter ${\kappa }_i$ acts in the reverse sense than in Case 1. One needs to balance $\kappa$ , $\delta$ and $n$ in the coefficients $c^1\left(n\right)\mathrm{,}{c}_i^2\left(n\right)$ in order to make the two coefficients of all quadratic terms positive. We are then left with question as to whether the norm $\Vert {\rho }_i\Vert$ is small against ${\alpha }_i$ for suitably large ${\alpha }_i$ . Now, the adjoint ${\rho }_i$ of the full network problem has as data the right hand side $-{\kappa }_i\left(q_i-q_i^0\right)$ and as $q_i$ -dependent coefficient ${g^{\prime }}_i\left(q_i\right)$ . For a given $q_i$ , the corresponding network equation is linear and by Lax-Milgram’s lemma, the solution ${\rho }_i$ satisfies an estimate against the data, which, in turn, depend on the original solution $q_i$ .

[Case 3.] $\lambda >0,\mu >0$ : In this case $a_i\mathrm{,}{b}_i$ in (63) need to be positive. This can be achieved in general if ${\alpha }_i$ large and ${\kappa }_i$ small and $\mu$ is large compared to $\lambda$ . 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={\left(q_i^l\right)}_{i\in \mathcal{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\mathrm{,}{p}_i^l$ converge to $q_i^0\mathrm{,}{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 $\alpha =10$ , $\kappa =10$ , $\nu =1$ and $f_i=\alpha \ast x$ and Neumann controls at all simple nodes. Clearly, the exact solution of the linear problem, i.e. $\gamma \ast g\left(q\right)$ , with $\gamma =0$ and $g\left(q\right)=\left|q\right|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 $\kappa$ . We compute the solution with nonlinearity $\gamma =0.1$ using the MATLAB routine bvp4 with tolerance $e\mathrm{.}-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 $\lambda \mathrm{,}\mu$ 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 ${\alpha }_i=1000$ , ${\kappa }_i=100$ , $\mu =1$ , $\lambda =0$ , $\gamma =0.1$ with relaxation parameter $ϵ=0.1$ . Figure 13 reveals the fact that due to the optimality of the functions $q_i\left(x\right)=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).

Conflicts of Interest

The authors declare no conflicts of interest.

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.