Analysis of District Heating Networks

In the context of district heating networks we consider a model for the distribution of energy through water (for heating or cooling) from a central power station to the consumers. We prove the well posedness of the system, by using the Banach Fixed Point Theorem together with stability estimates for reduced systems. Eventually we consider optimal control problems motivated by applications and we provide the existence of optimal controls in special situations.

station or combined heat and power plant (briefly CHP), to residential houses, commercial and industrial sites. The heat is transported along a network of pipes, filled with an incompressible fluid, typically water. Despite the initial costs due to the construction of the whole network, this system typically results in a high efficiency, in a smaller production of pollution compared to localized heat generators, and in reduced costs for final users in a long time scale.
Here we present, from a mathematical point of view, a complete teleheating system, modeled by ordinary and partial differential equations on a graph; see also [3,32]. Such equations describe both the fluid dynamics inside the pipes' network and the energy distribution towards the final consumers. Our main result consists in proving the existence of a unique solution to such a system and the continuous dependence of the solution from the initial and boundary conditions. This result opens the road for the study of optimal ways for heat distribution. In simple cases we also present optimal control problems and we show that optimal controls do exist; see [32] for finding optimal controls in a reduced case and Remark 3.11.
The model considered here is constructed starting from by the classic Euler system [29], which describes the evolution of a fluid through macroscopic quantities: density, linear momentum and energy. For results on Euler systems on networks we refer to [7-9, 12, 23] and references therein. In the present framework, since the fluid's incompressibility, in each pipe the velocity of the fluid is described by an ordinary differential equation. Furthermore at junctions coupling conditions can be given in a fixed form, since no change of characteristics can occur. For compressible flows this is a truly challenging problem.
In the present case, the heat energy satisfies a hyperbolic advection equation, since it is transported along each tube according to the fluid velocity. At junctions suitable coupling conditions provide the mass and energy conservation and distribution rules over the pipes. We use coupling conditions slightly different from those considered in [32] with the aim of avoiding low regularity of the solution. Note that although the energy equations are linear, nonlinear effects can occur due to the coupling. In [31,Sect. 3.2 ] a minimal example is given, how discontinuous solutions can occur even with smooth input data. In this example two neighboring consumers are directly supplied by the power plant, but also mutually connected by an additional pipe. The flow in this additional pipe is always pointing towards the consumer with the larger demand. Thus, if the demands change, the flow changes its direction and a different temperature can enter the pipe.
Partial differential equations on networks attracted a lot of interest in the last years, mainly due to applications, which range from traffic flow [25] to gas distribution [11,26], irrigation channels [17,21], sewer networks [2] and flexible strings [15,22,28]. In this setting, coupling conditions at nodes play an important role, since they strongly influence the existence, the well posedness, and the qualitative properties of the solution.
The proof of the main result is obtained through the classical Banach Fixed Point Theorem and through ad-hoc stability estimates. The most delicate part consists in deducing stability and total variation estimates for the energy subsystem. This is due to the fact that the energy evolves according to a transport equation whose characteristic curves move with speed given by the fluid velocity. Even for simple networks the direction of the flow can change arbitrarily often. This fact produces solutions for the energy functions whose traces at the nodes have low regularity. Such stability estimates e.g. are missing in [33].
The paper is organized as follows. In Sect. 2 we construct in detail the model starting from reduced systems, which describe the hydrodynamics part and the energy part. In Sect. 3 we state and prove the main result of the paper. Moreover we introduce the control problems and we show the existence of optimal controls in various situations. Finally, in the appendix, there are the detailed proofs of the stability of the advection equation w.r.t. the velocity of the flow and several auxiliary results.
Moreover we underline that the incoming and the outgoing edges are not determined by the sign of the velocity of the fluid inside the pipe, but only by the choice of the parametrization.
As stated in the introduction, in each pipe the starting point is the Euler system see [29,Chap. 1] where t is time, x ∈ [a, b] the spatial parametrization, ρ denotes the density of the fluid, v is the velocity, p is the pressure, and e is the energy density. In the case of an incompressible fluid, the macroscopic variable ρ is constant, so that the first equation in (2.1) simply becomes ∂ x v = 0 and (2.1) reduces to where the right hand side f = f (v) models the friction as well as gravitational effects, while the source term g = g(e) models the energy loss due to not perfect insulation.

Remark 2.3
Common choices for the source terms f and g in (2. 2) are The term − λ 2d v|v| corresponds to the Darcy-Weisbach friction law of the fluid with the pipe's walls (d is the pipe's diameter and λ is the friction factor). The source term −γ ∂ x h represents static pressure due to different elevations at the ends of the tube, where the constant γ is the gravitational acceleration and ∂ x h is the slope of the pipe. Finally, the term − 4κ d (T (e) − T ∞ ) models the energy loss due to heat flux over the pipe wall, where κ is the heat transmission coefficient of the pipe, T (e) is the fluid temperature depending on its energy density e and T ∞ is the external temperature.
In a district heating network, the heating power is distributed through water in a system of pipes. This can be described by considering a copy of (2.2) for each pipe I i , i ∈ I, of the network via suitable coupling conditions at nodes. The full system then can be divided into three different parts, namely the hydrodynamics part (briefly HYD), the consumer part, and the energy part. We describe them separately.

The Hydrodynamics Part
The flow in each pipe of the network, except those connected to a house, can be modeled by a copy of the first two equations of the complete incompressible Euler system (2.2), i.e. for every i ∈ I \ I H , by the system Note that, since the incompressibility assumption holds, namely ∂ x v i = 0, each v i is a function depending only on time; hence ∂ x p i depends only on time, so that p i is affine in x. For every i ∈ I \ I H , the initial condition for v i is given by Moreover, for every junction J ∈ J and t > 0, we consider the coupling conditions (2.6) The first equation states the conservation of mass. For constant density the mass flux in pipe j is the product of cross sectional area A j and fluid velocity v j . The pipes are assumed to be filled completely for all times. The remaining equations provide the continuity of the pressure at the junction. Finally, at the CHP we impose, for t > 0, the boundary condition Note that no initial condition for the pressure on the edges is needed. Indeed, thanks to the particular form of the second equation in (2.4) and to the compatibility condition (2.6), it is sufficient to provide the boundary condition at the CHP (2.7). For a later use, we consider the following assumptions:

The Consumers' Part
For every k ∈ I H , denote by J k the junction at position a k and let Q k (t) be the time dependent power demand of the consumer k. For every k ∈ I H , the velocity v k is determined by an ordinary differential equation depending on the power demand Q k .
More precisely, v k is governed by the equation where α > 0 is a relaxation parameter, e k,out = e(T k,out ) is the energy related to the fixed temperature level T k,out , to which the water inside the pipe is cooled down, and e J k is the energy at the starting node for consumer k. Note that equation (2.8) pushes We consider the following assumptions for the consumer subsystem: (C.1) There exist two positive constants Q max , Q min > 0 such that, for every k ∈ I H , Q k ∈ L 1 loc (0, +∞); Q min , Q max . Finally we give the definition of solution for the system (2.4)-(2.8).
3. v is a solution to (2.8) in the sense of Carathéodory [19], i.e. for every k ∈ I H and for every t

The Energy Part
The energy transport in the network is described by the third equation in (2.2). Thus for every i ∈ I, we consider the advection equation (2.9) supplemented with the initial condition and with suitable boundary and coupling conditions. More precisely, the boundary condition applies only at the level of CHP and, when admissible, it reads for a.e. t > 0. It prescribes the energy produced by the CHP. The coupling conditions at each internal node J ∈ J depend on the solution of a differential equation. Fix an internal node J ∈ J and t > 0. Define the (possibly time dependent) sets (see Fig. 1) Fig. 2 Schematic model of the node J with two incoming and one outgoing pipes, used in the derivation of the coupling condition (2.14) and consider the Cauchy problem where V J > 0 represents the volume inside the junction J . At the outgoing edges according to the fluid velocity, we impose that the energy is equal to e J : for t > 0 and for every i ∈ inc −,t (J ) and j ∈ out +,t (J ), we set (2.14)

Remark 2.5
The coupling conditions (2.13) can be justified in the following way. For J ∈ J , define V J > 0 the volume inside the junction J . Denote by e J (t) the temperature inside the junction J at time t, the variation e J of e J during the time step t is given by see Fig. 2. Thus, passing to the limit as t → 0, we obtain (2.13).

Remark 2.6
A reasonable alternative to the coupling condition (2.14) could be the following: stating the conservation of energy inside the node and the perfect mixture of energies; see also [3,32]. The advantage of conditions (2.13)-(2.14) is mainly technical. Indeed (2.13)-(2.14) permit to deduce total variation estimates needed in the proof of the main result.
We consider the following assumptions. We give now the definition of solution for the energy system. Definition 2.7 Fix T > 0 and assume (E.1) and (E.2). Furthermore, assume that for

The Complete System
Finally we consider the complete DHN system, composed by hydrodynamic and energy part. It is given by the system of differential equations with boundary conditions and with the coupling conditions (2.6) and (2.14).
We introduce here the concept of solution to the complete system.

Well Posedness Result
In this part we deal with the well posedness result for the district heating network, which is stated in Sect. 3.3. In Sect. 3.1 and in Sect. 3.2 we consider well posedness for the hydrodynamics and energy part respectively, which are the basic steps in the proof of the main result.
The proof is very similar to the proof of Theorem 5.3 in [27], but here we use less regularity assumptions on the inputs and the solution. We show that there is a decomposition of the system into pure algebraic and pure ODE parts, where the existence of unique solutions follows directly from standard results about Cauchy problems for solutions in Carathéodory sense [19].
Define the transformation where A T Q ∈ R N I ×N H rearranges the velocities of consumer edges v 0 ∈ R N H onto the full velocity vector. The columns of A T Q are the unit vectors of the consumer edges in I such that A Q A T Q = I d N H . The matrix A t ∈ R (M−1)×N I maps the edges of a spanning tree of G onto the full edge set and A PC ∈ R N I −(M−1)×N I is the incidence matrix of a set of fundamental cycles of G corresponding to the spanning tree.
Furthermore we define A = (A r , A p r ) as the incidence matrix of the whole graph, where A p r ∈ R N I ×1 is the column corresponding to the node of the CHP. A r is the remaining matrix for the inner nodes and consumers. A formal definition of these matrices can be found in [27]. For those matrices, we use the following properties: can be found in [27]. Property (iv) holds due to the fact that the rows of A Q only affect consumer edges while those are excluded in the definition of A t and A PC . System (2.4)-(2.8) can be compactly written as together with the corresponding initial conditions. Here, p J ∈ L 1 ((0, T ); R M ) is the vector of node pressures. Due to (2.6) this is well-defined. The last equation is the vectorized version of (2.8). We insert the transformation (3.2) and multiply by A t A PC from the left. This leads to the equivalent formulation (3.4) For system (3.4) the requirements for the existence theorem of solutions in Carathéodory sense are met, as the right hand side of the first equation is measurable in t as well as Lipschitz continuous in v 0 . We get a unique solution v 0 and due to the pure algebraic connection also a unique v 1 . Similarly, the ODE for v 2 has continuously differentiable right hand side and we get existence of a unique v 2 .
Applying Gronwall inequality provides us with bound after transforming back to the usual coordinates. The stability estimate (3.1) directly follows from basic ODE theory and Gronwall Lemma, e.g. [6,24].

Remark 3.2
Similar stability estimates can be deduced for the pressure p. We have, for t > 0, . (3.6) Using the fact that ∂ x p i does not depend on the space variable, we deduce that, for t > 0, and thus where the individual node pressures p J is taken from p J in (3.4). The form of the last equation in (3.4) immediately gives stability estimates for p J due to (H.2) and boundedness of v. Those estimates are similar to the ones for ∂ t v 2 . For the ease of notations and clarity of the proofs, we omit the estimates for p. Note that only the velocity couples to the energy transport. The triangular structure of (3.4) allows evaluation and stability estimates for v independently of the pressure, thus the results for the coupled system do not depend on p.

Remark 3.3
The constant L in Theorem 3.1 depends exponentially on T , since it is obtained through the Gronwall Lemma. More precisely where the Landau symbol O(1) denotes a suitable constant, which depends on initial data by max k∈I v k,0 , but not on T .

Energy Network
Here we consider system (2.9), supplemented with the initial conditions (2.10), with the coupling conditions (2.14), and with the boundary conditions (2.11).
Before we start with the proof, we state a lemma for the stability of the solution on a single edge.
, all with bounded total variation. Let G > 0 be fixed and let g : R → R be a Lipschitz continuous function such that |g(y 1 ) − g(y 2 )| ≤ G|y 1 − y 2 | for every y 1 , y 2 ∈ R. Let e 1 and e 2 be the solutions to the following IBVP problems: ∈ (a, b), t > 0 and P 2 : Then, for a.e. t > 0, the following stability estimate holds:

10)
where the constant K depends on the total variation and on the L ∞ -norm ofē 1 ,ē 2 , e L , e R , and on the L ∞ -norm of v 1 and v 2 .

for the fluxes at the boundaries we have
, (3.12) where . A sequential consideration with (3.10) can remove the bound on t.
The proof of Lemma 3.5 is given in Appendix A.1.

Proof of Theorem 3.4
Existence of solutions to such transport problems on networks have already been proven, e.g. in [33]. We focus on the stability estimate (3.9) that, to our knowledge, has not been shown before.
The proof consists of two steps. For both steps, we decompose the network problem into shorter time intervals. Those intervals are chosen in such a way that information can not reach one node from another, so that the network problem then decouples into a sequence of localized problems on segments. By defining the intermediate times 0 , where the constant v max denotes the maximum possible of velocity, in each interval [t i−1 , t i ] the information starting from one node can not reach a neighboring one. Therefore, a local consideration of the nodes with the adjacent edges is sufficient. To ease the notations, we assume here the topological orientation of the edges to point towards the node, such that positive velocities always represent an inflow, whereas negative velocities correspond to outgoing flows (different from the more complicated global definitions from (2.12)).
Step 1: Stability of node valuesē J andẽ J . From classical results about ODEs, we deduce that, for the junction J , the solutions e J andẽ J of (2.13) fulfill ·, a i ), andK and K J ,0 are suitable constants, not depending on time.
Inserting it into (3.13) yields Applying Gronwall inequality, we obtain and, due to the positivity ofK K J , Step 2: Stability estimates on the whole network Here, for suitable L 1 > 0 and L 2 > 0, we aim to prove: . (3.15) For each edge i ∈ I with I i = [a i , b i ] we can apply Lemma 3.5 and get , (3.16) where e L , e R are the left and right boundary values. In case the edge is connected to a node, these are the node energies e J .
Altogether we have i∈I ē i (t, ·) −ẽ i (t, ·) is the largest node degree in the network.
Consequently, (3.14) follows directly from (3.15). For a node J ∈ J and a single time interval [t i−1 , t i ], using the result of step 1, we obtain This is a recursive formula of the form with , So the explicit estimate of (3.20) is . (3.23) The result follows directly from (3.18).

Remark 3.6
The constant L in Theorem 3.4 depends exponentially on T . More precisely where the Landau symbol O(1) denotes a suitable constant, which depends on initial data, but not on T . Regarding the initial data, the constants contain only max i∈I e k,0 ∞ and the bound on the total variation TV e .

The Complete System
In this part, we deal with the well posedness result for the complete system. .

(3.24)
Proof Define the set Consider the operator where T is defined by two subsequent steps. For given initial values e 0 and external inputs e b we first we define by T e the operator producing a solution to the energy subsystem according to Theorem 3.4. The solution g = T e (v) is split into g = g I , g J , where g J denotes the node energy functions. Then, with v 0 and Q K given, the operator T v provides the solution to the hydrodynamic part (see Theorem 3.1). Its output we name By Theorem 3.4, there exists a unique g = g I , g J = T e (v), such that, in particular, g J ∈ L 1 ([0, T ]; R) for every J ∈ J (see Item 2 in Definition 2.7).
Since the set of internal nodes on pipes connected to the houses, i.e. pipes in the set I H , are a subset of J , by Theorem 3.1 there exists a unique w = T e (g J ) such that, in particular w ∈ AC([0, T ]; R N I ) and h∈J A h w h = 0 for every junction J ∈ J , that is w ∈ X v .
T is a contraction. Let T ≤ T . Fix two elementsv,ṽ ∈ X v and denotē

By Theorem 3.4 we deduce that
Therefore, using Theorem 3.1, we deduce that Thus the requirements of Banach fixed point theorem [1] are fulfilled, leading to the unique existence of a solution on [0, T ].
Since L H (T ) and L E (T ) do not depend on T , we aim to repeat this procedure to cover the full interval  with (ē,p,v) and (ẽ,p,ṽ) the corresponding solutions. By Theorem 3.1, we deduce that

By Theorem 3.4 we deduce that
. (3.28) In particular, the estimates of Theorems 3.1 and 3.4 hold point wise, such that Here we can now apply Gronwall inequality to get The monotonicity of the right hand side and (3.28) imply which concludes the proof.

Optimal Control Problems
In this part, we consider the complete teleheating system (2.15) as a control system, where the boundary term e 1,b at the CHP acts as a control function, and we assume that, given T > 0, the solution exists on the time interval [0, T ]. In this perspective, a natural problem consists in finding a control function e 1,b , satisfying assumption (E.1), which minimizes the functional for suitable coefficients α k ≥ 0, β k ≥ 0, and γ 1 , γ 2 , γ 3 , γ 4 ≥ 0. The minimization of the first two terms in (3.29) aims at producing the required temperature in the houses respectively on the whole time interval and at the final time. The third term in (3.29) measures the total energy produced by the CHP. The minimization of the fourth term penalizes too many oscillation in the energy production. The fifth and sixth terms have similar meanings with respect to the total power provided at the CHP. The next result deals with the lower semicontinuity of J . Then the functional J : E 2 → R, defined in (3.29), is lower semicontinuous with respect to the L 1 topology. If γ 2 = 0 and γ 4 = 0, then J : E 1 → R is also continuous.
Proof We prove that J : E 2 → R is lower semicontinuous with respect the L 1 topology. The final statement can be proved similarly. Fixē 1,b ∈ E 2 and a sequence e 1,b n ∈ E 2 such that e 1,b n →ē 1,b in L 1 (0, T ) as n → +∞. Denote with (ē,p,v) and respectively with (e n , p n , v n ) the solution to the system (2.15) corresponding to the boundary dataē 1,b and respectively to e 1,b n , according to Definition 2.8. Note that these solutions exist by Theorem 3.7.
Clearly, by Theorem 3.7, we have that Therefore, for every k ∈ I H , we deduce that The proof is so concluded.
The next result deals with optimal solutions.
Since K is compact, then there existsē 1,b ∈ K and a subsequence e 1,b n k such that e 1,b n k →ē 1,b as k → +∞. By using the lower semicontinuity of J (see Proposition 3.8), we deduce that concluding the proof.

Remark 3.10
Note that the previous result holds if K is a compact subset of E 1 (or of E 2 ) with respect the L 1 -topology. Therefore Corollary 3.9 can not be applied for general closed and bounded subsets K of E 1 or E 2 . If instead K is a closed and bounded subset of a finite dimensional subspace of E 1 or E 2 , then the functional J admits a minimum.

Remark 3.11
The results in the present section deal only with the existence of optimal controls which minimize the functional J in (3.29). Conversely, no necessary conditions, that represent the main tool in the search of optimizers, are deduced here. From the analytic point of view, this is a challenging problem since it needs some differentiability properties of the functional J with respect the topology of the control set. Dealing with solutions with low regularity makes the problem hard.
However one can try to find approximate optimal controls using numerical schemes inspired for example by the gradient descent method [4] as in [31], even if the present regularity does not fully justifies it. A possible not-optimized algorithm for approximating an optimal control could be the following.
1. Replace the control set K by a finite-dimensional approximation K d . 2. Fix an initial control e 0 ∈ K d . 3. Find the numerical gradient ∇ J (e 0 ) of J at e 0 and construct the control e 1 = e 0 − γ 0 ∇ J (e 0 ), possibly with a projection on the set K d , where γ 0 > 0 is a sufficiently small learning rate. 4. Recursively, given e n ∈ K d , find the numerical gradient ∇ J (e n ) of J at e n and construct the control e n+1 = e n − γ n ∇ J (e n ), possibly with a projection on the set K d , where γ n > 0 is a sufficiently small learning rate. and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A.1 Lipschitz continuity of solutions to advection equations
In this section we give, for an advection equation on a bounded interval, a proof of the stability of the solution and the traces of the flux at the boundary w.r.t. the speed. Note that we do require a uniform upper bound on the velocity, but do not restrict the orientation of the flow, i.e. the direction of the flow can change arbitrarily often.

Proof of Lemma 3.5
Observe first that the broad solution (i.e. the solution constructed using characteristics curves) to an initial-boundary value problem such as P 1 or P 2 is indeed a MV-solution, as introduced in item 3 of Definition 2.7. The proof of this claim follows the line of [10,Lemma 4.3], see also [13,Formula 4.5] for the definition of solution constructed through characteristics when the boundary datum is not zero. In order to prove the stability estimate (3.10), first we regularize the velocities v 1 and v 2 : there exist sequences Denote by e ν 1 , e ν 2 the corresponding solutions to P ν 1 and P ν 2 , constructed through characteristics. Using [34,Theorem 2.4], we deduce that, for a.e. t > 0, e ν i (t) converges to e i (t) in L 1 ((a, b); R), for i = 1, 2. The definition of MV-solution in item 3 of Definition 2.7 has the remarkable feature of being stable under convergence in L 1 , therefore we obtain that e i is the solution to P i , for i = 1, 2.
Using [34,Proposition 2] and adapting the proof of [34,Theorem 2.4] to the present linear setting, we obtain that for all t ∈ [0, T ] the following estimate holds where we used the bounds provided by Lemmas A.3 and A.4. Passing to the limit as ν → +∞, we obtain that (3.10) holds.
In the second part of this proof, we show (3.12), the proof for (3.11) being symmetric.
The case of positive velocity. Here, assuming that v 1 and v 2 are positive velocities, we prove that The proof is divided into three steps.
Step 1. Here only the initial datum is changed, namely, given v ≥ 0, e 1 and e 2 are solutions to Denote with the flux associated to the vector field g, i.e. solves d see [5]. Since v ≥ 0, for t > 0 sufficiently small, where X is the characteristic associated to v. Thus Step 2. Here we consider the same initial datum, but different velocities v 1 ≥ 0, v 2 ≥ 0, namely, givenē, e 1 and e 2 are solutions to Consider a sequenceē h ∈ C ∞ (]a, b[; R) such that Let θ ∈ [0, 1] and define the velocity Denote with X θ and e h θ respectively the characteristic associated to the velocity v θ and the solution to Observe that, for t ≥ t o and x o ∈ (a, b) We have Now, using these calculations and step 1, we deduce that Passing to the limit as h → +∞, we get Step 3. Here we consider different initial dataē 1 ,ē 2 and different velocities v 1 ≥ 0, v 2 ≥ 0, namely, e 1 and e 2 are solutions to Let e 3 be the solution to the additional problem Using the results of Step 1 and Step 2, we have that proving (A.1).

The general case.
For ε > 0 we define the sets and split the estimate for into the contributions coming from the different subsets of the whole domain For C ε , all involved velocities have an absolute value smaller than ε, leading to the estimate For the set B ε , one of the velocities is negative, meaning the corresponding contribution is zero and thus The crucial part is the estimate for A ε , where we can use the results for positive velocities obtained above. The set A ε is closed: due to the bounded total variation of v, it can be written as a finite union of closed intervals The goal is to get an estimate similar to (A.1) for each I ε i . Here it is important that the dependencies on the initial datum do not overlap, as that would lead to a boundlessly growing right hand side as ε tends to zero and the number of intervals grows. This is why we use a combination of above estimate and (3.10), where we take special care of the domains of dependency in the spatial integrals.
We start be considering the interval I ε n ε : by (A.1) we get is actually more precise than (A.1): indeed, we take care, in the L 1 -norm in space, of the true area of dependency induced by the flow velocity. In order to bound the first term on the right hand side of (A.3), we exploit (3.10): with c n ε = a n ε − v max (t ε n ε − t ε n ε −1 ). Due to v j > 0 in I ε n ε −1 , j = 1, 2, the last term in the inequality above can be estimated as follows: Using (A.6) for the time interval I ε n ε and (A.3) for I ε n ε −1 , leads to We remark that the coefficient of the norm of the difference between the velocities is computed taking care of the time intervals on which the L 1 -norm is evaluated. We can proceed by adding the time interval I n ε −2 : in particular, we exploit (A.6) for I ε n ε ∪ I ε n ε −1 and add (A.3) for I n ε −2 , noting that we can bound the last line on the right hand side of (A.6) exploiting (A.4) and (A.5). What we obtain is entirely analogous to (A.7). Proceeding in this way, we get Proof The proof follows almost directly using the fact that the solution to (A.8) along the characteristics lines is a solution of an ordinary differential equations with vector field given by g. Thus the amplification, obtained by the Gronwall Lemma, produced by the solution at time t is given by e Gt . Therefore the total variation of e(t) can be control by the total variation of the initial and boundary conditions multiplied by e Gt and by the jumps between the boundary data and the initial condition, again multiplied by the factor e Gt ; see also [34,Theorem 2.3] for a similar situation.
Here, deg(G) is the maximal node degree in the network and V the minimal volume of a junction. Using A.