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

$$\begin{array}{l}{\alpha }_i{q}_i^{l+1}\left(x\right)-{\beta }_i{\partial }_{xx}{q}_i^{l+1}\left(x\right)+g_i\left(q_i^{l+1}\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\\ d_{ik}{q}_i^{l+1}\left(n_k\right)+\sigma {\partial }_x{q}_i^{l+1}\left(n_k\right)\\ =\sigma \left(\frac{2}{d_k}{\displaystyle \sum _{j\in {\mathcal{T}}_k}}{\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 \sum _{j\in {\mathcal{T}}_k}}{d}_{ik}{q}_j^l\left(n_k\right)-d_{ik}{q}_i^l\left(n_k\right)\right)\\ =g_{ik}^{l+1}\mathrm{,}k\in {\mathcal{J}}^M\mathrm{.}\end{array}$$ (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)

$$+{\displaystyle \sum _{i\in \mathcal{I}}}{\displaystyle \underset{0}{\overset{{\mathcal{l}}_i}{\int }}}\left({\alpha }_i{e}_i^{l+1}{\varphi }_i+{\beta }_i{\partial }_x{e}_i^{l+1}{\partial }_x{\varphi }_i+\left(g_i\left(e_i^{l+1}+q_i\right)-g_i\left(q_i\right)\right){\varphi }_i\right)\text{d}x.$$ (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

$$\begin{array}{c}{\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 \sum _{i\in \mathcal{I}}}{\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+\left(g_i\left(e_i^l+q_i\right)-g_i\left(q_i\right)\right)e_i^l\right)\text{d}x\mathrm{.}\end{array}$$ (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

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

$$\begin{array}{l}\underset{u_j\mathrm{,}{v}_{ik}}{min}I\left(q_i\mathrm{,}{u}_j\mathrm{,}{v}_{ik}\right)\mathrm{<mo>:</mo><mo>=</mo>}\frac{\kappa }{2}{\Vert q_i-q_i^0\Vert }^2+\frac{\nu }{2}{u}_j^2+\frac{1}{2{\mu }_k}{v}_{ik}^2+\frac{1}{2{\mu }_k}{\left({\mu }_k{\partial }_x{q}_i\left(n_k\right)+h_{ik}^l\right)}^2\\ \text{subject}\text{to}\\ {\alpha }_i{q}_i-{\beta }_i{\partial }_{xx}{q}_i+g_i\left(q_i\right)=f_i\mathrm{,}x\in \left(\mathrm{0,}{\mathcal{l}}_i\right)\\ d_{ik}{q}_i\left(n_j\right)=u_j\mathrm{,}{d}_{ik}{q}_i\left(n_k\right)=-{\lambda }_k{\beta }_i{\partial }_x{q}_i\left(n_k\right)+g_{ik}^l+v_{ik}\mathrm{.}\end{array}$$ (42)

Case b.):

$$\begin{array}{l}\underset{u_j\mathrm{,}{v}_{ik}}{min}I\left(q_i\mathrm{,}{u}_i\mathrm{,}{v}_i\right)\mathrm{<mo>:</mo>}=\frac{\kappa }{2}{\Vert q_i-q_i^0\Vert }^2+\frac{\nu }{2}{u}_i^2+\frac{1}{2{\mu }_k}{v}_{ik}^2+\frac{1}{2{\mu }_k}{\left({\mu }_k{\beta }_i{\partial }_x{q}_i\left(n_k\right)+h_{ik}^l\right)}^2\\ \text{subject}\text{to}\\ {\alpha }_i{q}_i-{\beta }_i{\partial }_{xx}{q}_i+g_i\left(q_i\right)=f_i\mathrm{,}x\in \left(\mathrm{0,}{\mathcal{l}}_i\right)\\ d_{ij}{\beta }_i{\partial }_x{q}_i\left(n_j\right)=u_j\mathrm{,}{d}_{ik}{q}_i\left(n_k\right)=-{\lambda }_k{\partial }_x{q}_i\left(n_k\right)+g_{ik}^l+v_{ik}\mathrm{.}\end{array}$$ (43)

Case c.):

$$\begin{array}{l}\underset{v_{ij}\mathrm{,}{v}_{ik}}{min}I\left(q_i\mathrm{,}{v}_{ij}\mathrm{,}{v}_{ik}\right)\mathrm{<mo>:</mo>}=\frac{\kappa }{2}{\Vert q_i-q_i^0\Vert }^2+\frac{1}{2{\mu }_j}{v}_{ij}^2+\frac{1}{2{\mu }_k}{v}_{ik}^2\\ +\frac{1}{2{\mu }_k}{\left({\mu }_k{\beta }_i{\partial }_x{q}_i\left(n_k\right)+h_{ik}^l\right)}^2+\frac{1}{2{\mu }_j}{\left({\mu }_j{\beta }_i{\partial }_x{q}_i\left(n_j\right)+h_{ij}^l\right)}^2\\ \text{subject}\text{to}\\ {\alpha }_i{q}_i-{\beta }_i{\partial }_{xx}{q}_i+g_i\left(q_i\right)=f_i\mathrm{,}x\in \left(\mathrm{0,}{\mathcal{l}}_i\right)\\ d_{ij}{q}_i\left(n_j\right)=-{\lambda }_j{\beta }_i{\partial }_x{q}_i\left(n_j\right)+g_{ij}^l+v_{ij}\mathrm{,}{d}_{ik}{q}_i\left(n_k\right)=-{\lambda }_k{\beta }_i{\partial }_x{q}_i\left(n_k\right)+g_{ik}^l+v_{ik}\mathrm{.}\end{array}$$ (44)