Adaptive Refinement Strategies for the Simulation of Gas Flow in Networks using a Model Hierarchy

TL;DR

This paper introduces novel adaptive refinement strategies for gas flow simulation in pipeline networks, combining spatial, temporal, and model adaptivity to improve accuracy and reduce computational costs.

Contribution

It presents new greedy-like strategies for model adaptivity, extending existing spatial and temporal adaptivity methods, and demonstrates their superior performance in simulations.

Findings

Novel strategies outperform current methods in computational efficiency

The combined adaptivity approach effectively reduces simulation error

Strategies are validated through realistic gas flow simulations

Abstract

A model hierarchy that is based on the one-dimensional isothermal Euler equations of fluid dynamics is used for the simulation and optimisation of gas flow through a pipeline network. Adaptive refinement strategies have the aim of bringing the simulation error below a prescribed tolerance while keeping the computational costs low. While spatial and temporal stepsize adaptivity is well studied in the literature, model adaptivity is a new field of research. The problem of finding an optimal refinement strategy that combines these three types of adaptivity is a generalisation of the unbounded knapsack problem. A refinement strategy that is currently used in gas flow simulation software is compared to two novel greedy-like strategies. Both a theoretical experiment and a realistic gas flow simulation show that the novel strategies significantly outperform the current refinement strategy with respect to the computational cost incurred.

Tables

Table 1. Table 1: Cost functional constants in ( 4 ) and ( 5 ).

$m$	$C_m$	${\alpha }_m$	${\beta }_m$
$1$	$8.45\text{⋅}{10}^{-5}$	$0.952$	$0.937$
$2$	$1.06\text{⋅}{10}^{-4}$	$0.908$	$0.925$
$3$	$5.49\text{⋅}{10}^{-5}$	$0.694$	$0.857$

Table 2. Table 2: Mean total cost values in CPU seconds for strategies S ​ 1 𝑆 1 S1 - S ​ 3 𝑆 3 S3 and savings with respect to S ​ 1 𝑆 1 S1 . The subscripts 1 , 2 , 3 1 2 3 1,2,3 denote the different values of ϕ = 0.8 , 0.9 , 1 italic-ϕ 0.8 0.9 1 \phi=0.8,0.9,1 , respectively.

strategies	$S1$	$S{2}_1$	$S{2}_2$	$S{2}_3$	$S{3}_1$	$S{3}_2$	$S{3}_3$
Mean Total Cost	$36.1$	$8.30$	$8.04$	$7.72$	$7.64$	$7.39$	$7.12$
Savings (%)	-	$77.0$	$77.7$	$78.6$	$78.8$	$79.5$	$80.3$

Table 3. Table 3: Relative error ( tol = ​ 10 − 4 tol E-4 \text{tol}=${10}^{-4}$ ) and computational costs (in s s \mathrm{s} ) of a network simulation using strategies S ​ 1 𝑆 1 S1 - S ​ 3 𝑆 3 S3 and savings with respect to S ​ 1 𝑆 1 S1 . The subscripts 1 , 2 , 3 1 2 3 1,2,3 denote the different values of ϕ ∈ { 0.8 , 0.9 , 1 } italic-ϕ 0.8 0.9 1 \phi\in\{0.8,0.9,1\} , respectively.

Strategy	Relative Error	CPU time [s]	Savings
$S{1}_{}$	$1.384\text{⋅}{10}^{-5}$	7.53	-
$S{2}_1$	$3.730\text{⋅}{10}^{-5}$	2.42	$67.9\%\text{/}$
$S{2}_2$	$3.091\text{⋅}{10}^{-5}$	2.27	$69.9\%\text{/}$
$S{2}_3$	$3.091\text{⋅}{10}^{-5}$	2.29	$69.9\%\text{/}$
$S{3}_1$	$5.076\text{⋅}{10}^{-5}$	3.12	$58.6\%\text{/}$
$S{3}_2$	$5.190\text{⋅}{10}^{-5}$	2.83	$62.4\%\text{/}$
$S{3}_3$	$5.098\text{⋅}{10}^{-5}$	2.83	$62.4\%\text{/}$
reference	-	1013

Equations

$$\underset{r_i}{\min }$$

$$\eta -\sum _{i=1}^{3N_p}\sum _{k=1}^{r_i}{w}_{ik}\le {\text{tol}}_{\eta }.$$

$$\begin{array}{rl}\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|-g\rho h^{\prime },\end{array}$$

$$\begin{array}{rl}\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|-g\rho h^{\prime }.\end{array}$$

$$\begin{array}{rl}\rho v& =\text{constant},\\ p(x)& =\sqrt{p_{\text{in}}^2-\frac{\lambda c^{2}x}{D}\rho v\left|\rho v\right|}.\end{array}$$

$$\frac{{\sum }_{j\in {\mathcal{J}}_p}\left(e_{m,j}+e_{x,j}+e_{t,j}\right)}{\left|M({\mathbf{u}}^h)\right|}<\text{tol}.$$

$$e_{x,j}\left(r_{x,j}\right)\approx \frac{e_{x,j}}{2^{s_x\cdot r_{x,j}}},e_{t,j}\left(r_{t,j}\right)\approx \frac{e_{t,j}}{2^{s_t\cdot r_{t,j}}},\text{mbox}forallj\in {\mathcal{J}}_p,$$

$${\text{tol}}_m$$

$${\text{tol}}_x={\text{tol}}_t$$

$$b_j=\underset{j\in {\mathcal{J}}_p}{\max }\left\{\Delta e_{m,j},\Delta e_{x,j},\Delta e_{t,j}\right\}$$

$$b_j=\max \left\{\frac{\Delta e_{m,j}}{\Delta c_{m,j}},\frac{\Delta e_{x,j}}{\Delta c_{x,j}},\frac{\Delta e_{t,j}}{\Delta c_{t,j}}\right\},$$

$$F\left(m,n_x,n_t\right)=C_m\cdot n_x^{{\alpha }_m}\cdot n_t^{{\beta }_m},$$

$$F_c\left(m,r_x,r_t\right)=C_m{\left(2^{r_x}\cdot n_{x,0}\right)}^{{\alpha }_m}\cdot {\left(2^{r_t}\cdot n_{t,0}\right)}^{{\beta }_m}.$$

$$\sum _{j\in {\mathcal{J}}_p}\left(e_{m,j}+e_{x,j}+e_{t,j}\right)<\text{tol}\cdot M({\mathbf{u}}^h)=3.$$

Full text

††journal: Applied Mathematics and Computation

1 Introduction

The simulation of gas flows in pipeline networks is a topic of research that has been studied at various scales: from the individual pipeline to the entire network. Studies in control and optimisation of gas supply in a dynamic supply-demand environment strongly depend on large scale simulations of pipeline networks. In the last decades, considerable research on the modelling, simulation and optimisation of gas flow through pipeline networks has been conducted, see e.g. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Depending upon requirement, there exist multiple models to predict the system behaviour with varying levels of accuracy. Generally, more accurate models are computationally more expensive. Hence, in order to make real-time decisions, an appropriate trade-off between accuracy and computational complexity should be made. This can be achieved by using a hierarchy of models, where the models can be adaptively switched during the simulation process. Beside the models, the discretisation mesh may be varied in space and time, which places the demand for a strategy to automatically steer the simulation by changing the models and the discretisation meshes. This steering is based on simulation error estimates, which have been studied in detail in [14].

Since the simulations are the basis for decisions in the optimisation and control of the gas flow, the reliability of the simulation is of prime importance. The simulation is to be carried out such that the relative error in the state or in a functional of interest is below a specified tolerance. Starting with a coarse simulation, an adaptive strategy is used to bring the error below the tolerance by refining the discretisation in time and space or refining models, i.e., shifting to a model of higher accuracy. Hence, we have three different refinement possibilities for each pipe $j\in\mathcal{J}_{p}$ of the pipeline network, where $\mathcal{J}_{p}$ denotes the set of pipes in the network. These refinement possibilities are indexed by $i=1,\ldots,3N_{p}$, where $N_{p}\coloneqq\left|\mathcal{J}_{p}\right|$ is the number of pipes. Refinements are to be chosen such that the computational costs are kept low. We define an optimal refinement strategy as a strategy which returns the solution of the constrained optimisation problem

[TABLE]

Here, for each refinement possibility $i$, $r_{i}$ is the number of refinements, $w_{ik}$ is the relative error reduction due to the $k^{\text{th}}$ subsequent refinement and $v_{ik}$ is the corresponding cost addition. The constants $c$ and $\eta$ denote the cost and relative error of the starting simulation, respectively. We note that if $v_{ik}$ and $w_{ik}$ are constant for all $k$, then this problem is equivalent to the unbounded knapsack problem which is NP-hard, see e.g. [15]. In this paper we aim to find a good approximation to the solution of this generalisation of the knapsack problem (1). For this, we propose three adaptive refinement strategies which return approximate solutions of (1). The ideas presented here are for the example of pipeline networks. By generalising pipes to functional sub-domains, the principles of adaptive refinement can be extended to simulations for other applications which use a model hierarchy, e.g., power grids and water supply networks.

This paper is organised as follows. Section 2 introduces a model hierarchy for the simulation of gas flow through a single pipe as well as the aim of the adaptive refinement strategies. Section 3 describes the three proposed refinement strategies and Section 4 introduces the design of a synthetic experiment. The results of this synthetic experiment are contained in Section 5 and an application of the refinement strategies to a realistic gas network simulation is given in Section 6. Finally, some conclusions are contained in Section 7.

2 Simulation of Gas Flow

In this section, we give an example of a model hierarchy for gas flow simulations. We then outline a framework for the adaptive simulation of pipeline networks using the given model hierarchy. Finally, we highlight the aims of the refinement strategies in the simulation framework.

2.1 The Model Hierarchy

As an example, we take a three-model hierarchy for the gas flow simulations, discussed in detail in [16]. At the highest level, we have the isothermal Euler equations relating the gas density $\rho$ and gas flow $\rho v$ by

[TABLE]

together with the equation of state for real gases $p=\rho z(p)RT$ with compressibility factor $z(p)=1-\alpha p$ for a constant $\alpha\in\mathbb{R}^{+}$. Here, $p$ denotes the pressure, $v$ the velocity of the gas, $R$ the specific gas constant, and $T$ the temperature. Further, $\lambda>0$ is the Darcy friction coefficient, $D$ the pipe diameter, $g=$9.81\text{,}\mathrm{m}\text{/}{\mathrm{s}}^{2}$$ the acceleration due to gravity, and $h^{\prime}$ the slope of the pipeline. If the term $\frac{\partial}{\partial x}\left(\rho v^{2}\right)$ is small, it can be dropped resulting in a semilinear model

[TABLE]

A further simplification of assuming a stationary state and zero slope $h^{\prime}=0$ yields a system of two ODEs, which can be solved analytically and are referred to as the algebraic model

[TABLE]

Here, $c=\sqrt{p/\rho}$ denotes the speed of sound within the gas and $p_{\text{in}}=p(0)$ the input pressure. The three models are shown in hierarchical form in Fig. 1. The model hierarchy is set in the decreasing order of accuracy for our purpose of gas flows. Each pipe in the network is simulated using one of these three models and varying discretisation stepsizes in space and time.

2.2 Adaptive Gas Network Simulation

We consider a gas flow simulation over a pipeline network $\mathcal{J}_{p}$. The simulation time $\left[0,T\right]$ is divided into time intervals of equal size $\left[t_{k},t_{k+1}\right]$, $k=0,1,\dots,N\!-\!1$ and $t_{N}=T$. Given a starting model distribution over the network $\mathbf{m}_{0}=\left[m_{0,1},m_{0,2},\ldots,m_{0,N_{p}}\right]^{T}$, with $m_{0,j}\!\in\!\{1,2,3\}$, and a corresponding discretisation in space $\mathbf{n}_{x,0}$ and in time $\mathbf{n}_{t,0}$, a simulation is run for $\left[t_{0},t_{1}\right]$. We obtain error distributions along the network using a posteriori error estimations (see [14]): $\mathbf{e}_{m}=\left[e_{m,1},e_{m,2},\ldots,e_{m,N_{p}}\right]^{T}$ for the model errors and $\mathbf{e}_{x}=\left[e_{x,1},e_{x,2},\ldots,e_{x,N_{p}}\right]^{T}$ and $\mathbf{e}_{t}=\left[e_{t,1},e_{t,2},\ldots,e_{t,N_{p}}\right]^{T}$ for the spatial and temporal discretisation errors, respectively. These error estimators are derived using a dual weighted residual method with a user-defined functional of interest $M(\mathbf{u})$, where $\mathbf{u}=\left[\rho,\rho v\right]^{T}$, see [14]. The simulation error for a single pipe is the sum of all three errors. For the simulation to be valid, the relative error must be below a given tolerance tol. Hence, we require that

[TABLE]

If the tolerance is not achieved, models and discretisation meshes are refined. The task of deciding the required refinements is made by an adaptive strategy. A switch to a higher model in the hierarchy is called a model refinement and a refinement of the mesh is called a discretisation refinement. With the new models and discretisations we re-simulate for the time interval and continue the cycle. Once the solution meets the tolerance requirements, the models and discretisations are coarsened where possible, the simulation progresses to the next time interval and the cycle repeats. This simulation flow is shown in Figure 2 for an interval $\left[t_{k},t_{k+1}\right]$.

2.3 Structure and Aim of Adaptive Strategies

Our focus lies on finding adaptive strategies that control the errors and drive the simulation. For the time interval $\left[t_{k},t_{k+1}\right]$, an adaptive strategy takes as input the error distributions $\mathbf{e}_{m},\mathbf{e}_{x},\mathbf{e}_{t}$, the model distribution $\mathbf{m}_{k}$ and the number of nodes in the spatial and temporal discretisations $\mathbf{n}_{x,k},\mathbf{n}_{t,k}$. The strategy returns a refinement scheme $\mathbf{r}_{m}=\left[r_{m,1},\ldots,r_{m,N_{p}}\right]^{T}$, $\mathbf{r}_{x}=\left[r_{x,1},\ldots,r_{x,N_{p}}\right]^{T}$ and $\mathbf{r}_{t}=\left[r_{t,1},\ldots,r_{t,N_{p}}\right]^{T}$, where $r_{m,j}\!\in\!\{0,1,2\}$ and $r_{x,j},r_{t,j}\!\in\!\mathbb{N}$ denote the number of refinements to be made in the models and in the discretisations for all pipes $j\in\mathcal{J}_{p}$ such that constraint (2) is satisfied. The aim of the adaptive strategies is to achieve this constraint while keeping the computational costs that are incurred in the simulation low.

3 Refinement Strategies

In this section, we discuss three strategies for the adaptive refinements in the network simulation. Errors $\mathbf{e}_{m},\mathbf{e}_{x},\mathbf{e}_{t}\in\mathbb{R}^{N_{p}}$ and refinements $\mathbf{r}_{m},\mathbf{r}_{x},\mathbf{r}_{t}\in\mathbb{N}^{N_{p}}$ are assigned to every pipe $j\in\mathcal{J}_{p}$. One refinement in the spatial mesh or temporal discretisation is defined to be taking the new stepsize as half of its previous value. Hence, approximate relations between the initial errors $e$ and the errors after $r$ refinements $e(r)$ are given by

[TABLE]

where $s_{x}$ and $s_{t}$ are the convergence orders of the spatial and temporal discretisation schemes. The reduction for the error is only an approximation. Therefore, to have a safe upper bound for the estimated error after refinement, we multiply these approximated errors by a safety factor of refinement $f_{r}>1$ in each of the pipes that require refinements to be made. This ensures that it is very unlikely that the error overshoots its estimated value after the refinement.

In the model hierarchy presented in Fig. 1, discretisation errors feature only in the models M1 and M2. The algebraic model M3 has no discretisation errors. Thus, when models are switched from M3 to M2, discretisation errors are introduced. For pipe $j$ simulated with the most detailed model M1, we set the model error to $e_{m,j}=0$.

Strategy 1 - Individual Bounds $\left(S1\right)$

In order to meet the tolerance of the network error, we derive fixed individual error bounds for individual pipes for each of the three error types. The simulation is then carried out such that for each pipe the error is below the individual bound for all three errors.

We set the tolerance for the model errors as

[TABLE]

The remaining tolerance is equally divided between the spatial and temporal discretisations, i.e.,

[TABLE]

are the bounds for both error types for the entire network. To get the bounds for individual pipes, we uniformly distribute these bounds over the entire network, i.e., we divide them by the number of pipes ${N_{p}}$.

For the refinements, first a discretisation refinement is computed to bring the errors below the respective tolerance for each pipe. Subsequently, if the network error still exceeds the tolerance, models are refined to the next model higher up in the hierarchy. The simulation errors are then re-evaluated and this cycle is repeated until (2) is satisfied. This strategy is discussed in detail in [14]. The pseudocode for Individual Bounds is given in Algorithm 1.

Strategy 2 - Maximal Error Refinement $\left(S2\right)$

Since Strategy $1$ assigns individual tolerances, it loses the view of the network as a whole. However, the contribution of different errors to the overall network can balance each other without overshooting the total network tolerance. This is accounted for in the following strategy where we seek to make only those refinements which result in the maximal error reduction. This results in an iterative procedure. For every pipe $j$, we compute in every iteration the error reduction due to a single refinement in the model $\Delta e_{m,j}$ and in the space and time discretisations $\Delta e_{x,j}$ and $\Delta e_{t,j}$. The best option

[TABLE]

is passed to the network. On the network level, we mark those refinements for which the error reductions are larger than $\phi\cdot\max\left\{b\right\}$, with $\phi\leq 1$. This iteration is repeated until the network error is brought below the tolerance. The function maximalErrorRefinement, see Algorithm 2, represents the network controller and the function errorReduction, see Algorithm 3, represents the pipe level computations.

The spatial and temporal discretisation errors also depend on the simulation model. In lines 4 and 5 of Algorithm 3, we use $F_{x}\!\left(m_{c}\right)$ and $F_{t}\!\left(m_{c}\right)$, with $m_{c}$ the current model, as error amplification factors for the discretisation errors which account for this model dependency. Here, $m_{c}\!=\!1,2,3$ refers to models M1, M2 and M3, respectively. Since the discretisation error is absent for $m_{c}\!=\!3$, we set $F_{x}\!\left(3\right)\!=\!F_{t}\!\left(3\right)\!=\!0$. Furthermore, we set $F_{x}\!\left(2\right)\!=\!F_{t}\!\left(2\right)\!=\!1$. For model M1, the amplification factor for the discretisation errors is set with respect to the benchmark model M2. The factor $F_{m}\!\left(a,b\right)$ denotes an error reduction for the model error when models are shifted from $a$ to $b$. The notation $\left[b,z\right]=\texttt{max}\{\cdot\}$ is similar to MATLAB notation where $b$ denotes the maximal element and $z$ is the corresponding index. Note that in determining $\Delta e_{m}$ we also consider changes in the spatial and temporal discretisation errors. The central idea is to account for net error reduction.

Strategy 3 - Maximal Error-to-Cost Refinement $\left(S3\right)$

The adaptive refinements are made with an objective of reducing the computational cost without compromising on the simulation error. However, the previous two strategies do not address the computational costs explicitly. They address the error tolerance which is merely a constraint to the adaptive strategies, viewed in the optimization setting (1). In Strategy 3, however, we also take into account the computational costs that are incurred by the refinement. The idea of this strategy is similar to the greedy approximation algorithm for solving the unbounded knapsack problem, see [17]. Strategy 3, given in Algorithm 4, is similar to Algorithm 2 on the network level. On the pipe level, however, we compute the cost additions $\Delta c_{m,j},\Delta c_{x,j},\Delta c_{t,j}$, $\mbox{ for all\ }j\in\mathcal{J}_{p}$, using the cost functional $F_{c}\left(m_{c},r_{x},r_{t}\right)$, for each of the error reductions $\Delta e_{m,j},\Delta e_{x,j},\Delta e_{t,j}$. The error controller on pipe level then passes the best option

[TABLE]

i.e., the maximal error-to-cost ratio, to the network. This extension in the error controller on pipe level is given in function errorToCostReduction in Algorithm 5.

4 Design of Experiment

We tested the performance of the three refinement strategies on $10^{4}$ random samples of error distributions. The random samples represent errors in a simulation of gas flow in a given network of ${N_{p}}=12$ pipes. A better strategy will lead to lower computational costs while reducing the simulation error to a level below the tolerance. We compute the computational cost per pipe using a cost functional of the form

[TABLE]

where $m\in\{1,2,3\}$ denotes the model and $C_{m}$, $\alpha_{m}$ and $\beta_{m}$ are model-dependent constants. These constants, given in Table 1, are determined by the method of least squares. For this, gas flow simulations through a single pipe are performed using the software ANACONDA (cf. [18, 19]) with many different values of $n_{x}$ and $n_{t}$, which return the corresponding computational cost values $F$. A plot showing the computational costs of simulating a single pipe using different models depending on the number of nodes in space and time is depicted in Fig. 3. We can rewrite the functional (4) in terms of refinements assuming that the initial number of nodes $n_{x,0},n_{t,0}$ are known. Then we get

[TABLE]

The initial number of space and time discretisation nodes are chosen from the interval $\left[100,200\right]$ using a random number generator. All models are set to the most simple model M3 in the beginning. The error reduction upon model refinement also takes into account the introduction or increase of the spatial and temporal errors. This requires that the spatial and temporal errors are small when compared to the model error. Hence, for the experiment, the initial model errors are taken from $\mathcal{U}\left[0,1\right]$, where $\mathcal{U}\left[a,b\right]$ denotes a uniform probability distribution on the interval $\left[a,b\right]$, and the spatial and temporal discretisation errors are taken from $\mathcal{U}\left[0,0.2\right]$. For the model and discretisation errors after refinement, the approximate error reductions in lines 4, 5 and 8 of Algorithm 3 are used, where we choose $F_{m}(3,2)\!=\!3/4$ and $F_{m}(2,1)\!=\!1/4$. Models M1 and M2 are discretised with the implicit box scheme, as in [20], for which it holds that $s_{x}\!=\!2$ and $s_{t}\!=\!1$. The parameter $\kappa=1/3$ for Strategy 1 is chosen such that all three errors have an equal fraction of the tolerance. Strategies 2 and 3 are tested for a fraction $\phi\in\left\{0.8,0.9,1\right\}$ of the maximal best option. The strategies work for a relative error tolerance of $\text{tol}=${10}^{-1}$$ with a target functional value $M(\mathbf{u}^{h})=2.5\cdot{N_{p}}=30$. Hence, we require for the total network simulation error that

[TABLE]

The results of this experiment are given in the next section.

5 Results and Discussions

Each strategy returns a refinement scheme which brings the simulation error below the tolerance. The goal is to have low computational costs. The mean of the total computational cost values in CPU seconds over $10^{4}$ samples is shown in Table 2. We show the percentage savings in mean total computational cost of the strategies with respect to Strategy $1$. We denote the strategies as $S1-S3$. The subscripts $1,2,3$ refer to $\phi=0.8,0.9$ and $1$, respectively. We observe that strategies $S2$ and $S3$ have a percentage saving of over $77\%$ with respect to $S1$ for all values of $\phi\in\left\{0.8,0.9,1\right\}$. Among the different values, $\phi=1$ performs best for both $S2$ and $S3$.

Thus, by working with a greedy-like strategy for error control, an adaptive process can reduce the computational cost significantly. Furthermore, accounting for the computational cost explicitly in our estimates, we find even better refinement schemes that result in lower computational costs.

6 Application to a Realistic Network Simulation

We now apply the three different strategies to a simulation of a gas supply network, which is shown in Fig. 4. The considered network consists of twelve pipes (P01 – P12, with lengths between 30km and 100km), two sources (S01 – S02), four consumers (C01 – C04), three compressor stations (Comp01 - Comp03) and one control valve (CV01). Starting with stationary initial data, the boundary conditions and the control for the compressor stations and the control valve are time-dependent. The simulation time is $14\,400\text{\,}\mathrm{s}\text{/}$. The target functional $M(\mathbf{u})$ is given by the total fuel gas consumption of the three compressors and the error estimators are evaluated using a dual weighted residual method as developed, for example, in [21, 22, 23, 24]. For details of the derivation of the error estimators we refer to [14, 25]. The simulation is performed using the software ANACONDA (cf. [18, 19]).

Remark 1**.**

For the strategies proposed in Section 3, the temporal error was considered individually for each pipe. Since in the implementation of ANACONDA, the time stepping is uniform for the entire network, the temporal error $\eta_{t}$ was computed globally and divided by the number of pipes in order to get a local temporal error. However, if a best option $B_{j}$ was supposed to be the temporal error and had to be refined then all pipes were refined in time uniformly and the best options $B_{j}$ were updated.**

A reference solution was computed using model M1 and a very fine discretisation. The strategies were run with a relative tolerance of $\text{tol}=${10}^{-4}$$. Table 3 shows the relative error of the simulation compared to the reference solution, the CPU time taken and the percentage savings of strategy $S2$ and $S3$ in relation to strategy $S1$.

Compared to the synthetic experiment in Section 5, we see that the savings of strategies $S2$ and $S3$ applied to the simulation are in a similar range. The choice of the parameter $\phi\in\{0.8,0.9,1.0\}$, however, does not seem to have a significant influence on the saving. Moreover, the Maximal Error-to-Cost strategy $S3$ does not result in a larger saving of CPU time. What is noticeable is that the relative errors of the strategies $S2$ and $S3$ are closer to the proposed relative tolerance, which shows that they are not as restrictive as the Individual Bounds strategy $S1$.

7 Conclusions

In this paper we address the problem of automatic error control for large scale gas flow simulations that use a model hierarchy. The simulation needs to be reliable, i.e., keeping the total relative error below a specified tolerance, while retaining low computational costs. The problem of finding an optimal refinement strategy is a generalisation of the knapsack problem. We present three strategies for adaptive simulation error control via spatial and temporal discretisation mesh and model refinements. The strategy Individual Bounds, which is currently implemented in ANACONDA, sets a uniform bound for each error type and each pipe, Maximal Error Refinement iteratively chooses those refinements that result in the largest error reduction and has a network overview, and Maximal Error-to-Cost Refinement also accounts for the increase in computational cost inflicted by the refinement.

We constructed a synthetic experiment to test the three strategies. From this experiment we see that the two greedy-like strategies significantly reduce the computational cost as compared to the Individual Bounds strategy. This result is largely reflected in an actual gas flow simulation using ANACONDA for a $12$ pipe network including compressor stations and a control valve. Especially when the simulation process is a key component in a gas flow optimisation problem, the novel refinement strategies lead to considerable computational savings without compromising on the simulation accuracy.

Acknowledgements

The authors would like to thank the Deutsche Forschungsgemeinschaft for their support within the projects B01 and B03 of the Collaborative Research Centre/Transregio 154 Mathematical Modelling, Simulation and Optimization using the Example of Gas Networks. Jens Lang was also supported by the Excellence Initiative of the German Federal and State Governments via the Darmstadt Graduate School of Excellence Energy Science and Engineering.