On the Numerical Treatment of Interlaced Target Values – Modeling, Optimization and Simulation of Regulating Valves in Gas Networks

,


Introduction
The physical and technical properties of gas transport in networks have been studied for many decades already.They have gotten increasing attention in the last years since countries worldwide attempt a turnaround in their national energy policies towards carbon dioxide neutral means of energy production.
Reasons for this are the potential usage of gas networks as flexible energy storage for balancing the highly volatile generation of renewable energy sources [Federal Ministry for Economic Affairs and Energy, 2021] or the prospect of shifting to hydrogen transport in the future [European Commission, 2020].In both cases, controlling the network is expected to become more dynamic and complex, either due to the increasing variability of supply and demand or the gas properties of hydrogen compared to natural gas [Hoppmann-Baum et al., 2020].To make the task of operating gas networks manageable, the decisionmaking is split into multiples levels.It ranges from the complete view on the network featuring only abstract representations of the technical elements to detailed but local control decisions, for example, the degree of openness of a single pressure regulator or the rotational speed of a gas compressor [Stelter, 1988].The decisions on the different levels are usually connected via the concept of set-point values or target values.These represent a control requirement determined in the full network view, which then serves as an objective for the local operation of the actual machinery.
Controlling an element with set-points is a well-known concept in the community of automated control systems.The difference between the measured (actual) and the desired (target) value of an observed process variable is the basis for error-controlled regulation with negative feedback (closed-loop) in automatic control [Porter and Khaki-Sedigh, 1988].Proportional-Integral-Derivative (PID) control is the most common feedback control algorithm in engineering systems.In process industries, most of the control loops (≈ 95%) are of PID type [ Åström and Murray, 2008].The advanced control methods including model based control techniques (model predictive control [Camacho and Alba, 2007;Darby and Nikolaou, 2012;Haber et al., 2014] or multivariable control [Skogestad and Postlethwaite, 2005]), intelligent control techniques (fuzzy control [Nguyen et al., 2019] or neural network based control [Meireles et al., 2003]), and adaptive control [ Åström and Kumar, 2014], are often used as a high level control procedures for designing high performance controllers that can be applied to high-order and multivariable processes, typically nonlinear and subject to constraints.
When looking at work regarding the simulation or optimization of gas networks, the usage of setpoints to control active elements has a long history as well.In [Mallinson et al., 1993], the authors present the optimization problem of finding a stationary control for the British gas network.The control of the compressors is defined as the corresponding set-point for the gas pressure after compression, also called the outlet or outgoing pressure of the element.The same holds for [Rachford and Carter, 2000], which examine a transient optimization problem on a gunbarrel network, given a desired final state of the network.In a second paper, the same authors extended the model to prepare for not one but multiple possible future scenarios in a robust way [Carter and Rachford, 2003].In recent years, more set-point or target value models have been used.In their description of a stationary gas network control optimization problem, [Pfetsch et al., 2014;Schmidt et al., 2015;Koch et al., 2015] introduced a model for a regulator without remote access having a fixed set-point for the outgoing pressure.However, depending on whether the actual state value exceeds the set-point, matches it, or falls below it, the regulator is closed, active, or fully open.In the last case, it acts as a bypass.Finally, there are also models including a multitude of different set-point or target values.In [Pambour et al., 2016], a model used for gas network simulation was introduced, in which regulators and compressors are given a set-point either for the inlet pressure, the outlet pressure, the flow, or the ratio of both the in-and outlet pressure.Alternatively, they can be chosen to be closed or in bypass mode.The simulation model of [Benner et al., 2019] features two different modes for each compressor and regulator.These either try to keep the inlet or the outlet pressure at a given target value.
In this article, we introduce a much more involved description for the target value control of active elements.It features up to 6 different set-point values, which are in simultaneous use and hence competing.As a result, the element can react to changing conditions in the surrounding network without the need to switch the target value objective manually.For example, it is able to autonomously switch between holding the inlet pressure, holding the flow, or even closing the element, in the case of a regulator.This model captures the behavior of complex elements used in real-world networks and industry-standard simulators to precisely express the desired control given dynamically changing conditions in the corresponding part of the network.
The description can, in general, be applied to regulators and compressors.However, we will focus on regulators here, i.e., the elements reducing the pressure in flow direction by reducing the element's throughput.In literature, there exist several models for regulators.For example, [Huck, 2018] models an idealized regulator by demanding that the pressure difference between its outlet and inlet pressure must be positive and inside a given interval.In [Benner et al., 2019], the regulator is described as a resistor with variable diameter.Hence, it is given a degree of openness ranging between 0 and 1.This allows them to include the edge case of a fully closed regulator, in which the two pressures are decoupled and the flow is zero.The model presented in [Pfetsch et al., 2014;Schmidt et al., 2015;Koch et al., 2015] distinguishes two types of regulators: Those with and without remote access.The ones without remote access are controlled as already described above using the outlet pressure target value.When having remote access, the regulators can either be closed or again keep their pressure difference inside a given interval.Another and precursor to one of our models presented here is presented in [Hennings et al., 2021].Within there, the authors consider three different modes for a regulator: The regulator is either closed, active, or open.A closed regulator is defined as for the other models, an open one does assume the inlet and outlet pressures to be equal, while an active regulator allows for arbitrary pressure reduction.In addition, the regulator features a check valve or flap trap, which prevents flow against the regulator's orientation.
We will present two models for regulators using the target value control: One model suitable for dynamic simulation using a given set of target values over time and one for discrete optimization problems, in which target values should be determined that induce the desired element behavior while changing as little as possible.
The rest of the paper is structured as follows: In Section 2, we will first introduce the fundamental gas flow dynamic in pipes and our representation of a gas network as a directed graph.Based on this, we describe our new regulator model based on the 6 different target values in Section 3. In the following two sections, we first introduce a corresponding model suitable for the transient simulation of regulators given the target value settings for a future time horizon.Afterwards, we describe a discrete optimization model that can determine the control of a regulator in terms of the target values in a transient context.
In Section 6, we present two experiments on a minimal gas network and compare the performance of the two models against each other as well as a commercial simulation software.We conclude with the outlook in Section 7.

Models for Gas Flow in Pipe Networks
Mathematical models of various so-called active elements, such as regulators, compressors, or valves, manipulate the gas flow through themselves by blocking, resisting, or enhancing it.However, the response to their actions is determined in the surrounding pipes that provide the volume for the gas to be.Hence, we need to model the behavior of the gas in pipes in terms of pressures p, mass-flows q, and their interaction.We will use a version of the isothermal Euler equations, which we will introduce and describe in the following subsection.

Isothermal Euler Equations
A simplified or friction dominated version of the one-dimensional set of isothermal Euler equations for the description of the behavior of fluids within perfectly cylindrical pipes of length with χ ∈ {0, 1} being a time constant model switch.These particular versions are referred to as ISO2 if χ = 1, or as ISO3 if χ = 0 according to [Domschke et al., 2021] and can be derived by a time scaling approach presented by [Brouwer et al., 2011].According to empirical observations in [Hennings, 2021], we will use χ = 0 exclusively here.The intrinsic quantities of the equations are p ≡ p(x, t) for the pressure, q ≡ q(x, t) for the mass flow, ρ(x, t) for the density, and v(x, t) for the velocity and are parameterized along the longitudinal axis by x ∈ Ω = [0, L] as well as being averaged across the cross-sectional area A = 1 4 πD 2 as enclosed by the pipe.Furthermore and throughout, let ṗ(x, t) ≡ ∂ t p(x, t) and likewise for q(x, t) be the time derivatives.Also, let R s be the specific gas constant, T be the time-constant temperature, κ = R s T /A as shorthand, g ≈ 9.81 m/s 2 the gravity constant, and δh = (h r − h )/L as secant slope defined by the heights h on the left and h r right end of the pipe.Furthermore, let z(p) = z 0 + z 1 p + • • • + z n p n be a polynomial model or truncated Virial expansion for the gas factor for real gases [Onnes, 1902], and λ be the friction coefficient of the Darcy-Weisbach equation [Brown, 2003].The latter is an empirical model of pressure loss due to friction with the pipe wall.
For the determination of coefficients of the polynomial real gas factor model z(p), we will make use of the linear AGA formula [SIMONE Research Group and LIWACOM Informationstechnik GmbH, 2021;Domschke et al., 2021] of the American Gas Association, being accurate up to pressures p ≤ 70 bar, the quadratic formula of Papay [Papay, 1968;Saleh, 2002] for pressures p ≤ 150 bar, or the constant model z 0 = c 2 /(R s T ) as suggested by [Domschke et al., 2021], where c is the speed of sound within the given gas mixture.
To approximate the friction factor λ, we will use the Nikuradse formula [Nikuradse, 1950], which is given by It is an explicit simplification of the otherwise implicit Colebrook-White [Colebrook and White, 1937;Brown, 2003] formula.Here r is the roughness of the pipe wall surface and is usually provided in millimeters mm.Two important identities that have already been utilized for the ultimate derivation of (1) are where the first one is referred to as the state equation for real gases [Menon, 2005].We call (1) the pipe equations.More in-depth details on Euler equations in the context of fluid transport can be found, e.g., in [Osiadacz, 1996;Domschke et al., 2021;Benner et al., 2019].
Alternatively, in [Hennings, 2018], a linearization of the pipe equation (1) has been proposed, which has the benefit of avoiding divisions by small pressure values.There it was observed that the equations (3) yield the identity |v| = z(p) • κ • |q|/p, which can be fixed on predetermined or forecast velocities |v|, such that Provided we use a constant model for the real gas factor, then equations (1a) and (1b") together yield a linearized model of the pipe equations.

Gas Networks as Graph Structures From a Macroscopic Perspective
We have just established the system of pipe equations ( 1), but we are interested in gas networks as a whole and in solving online tasks.In this context, online means to utilize simulation or optimization to formulate control and operation recommendations as a guiding as well as continuous process live at operation time.Thus, the overall problem results in potentially large systems that must be solved repeatedly in a very limited time frame.For this reason, we have to take a macroscopic viewpoint regarding modeling the network and its building blocks.Hence, we consider gas networks as directed graphs G = (N , A) over some time horizon [t 0 , T ] ⊆ R of interest.The set of nodes N is finite |N | < ∞, and the set of arcs A ⊆ N × N is a collection of ordered tuples of nodes.
Nodes represent the intersection points of all the arcs.Each node n ∈ N serves as a boundary to all of its adjacent arcs and hence intertwines and connects them by the so-called flow balance equation, which is a Kirchhoff typed law q a (0, t). (4) Here A [n] is the set of all arcs (i.e., ordered tuples) a ∈ A whose first node is n and A r [n] likewise is the set of all arcs a ∈ A whose second node is n, i.e.
Furthermore, q (n) is the gas inflow or outflow to the network, and therefore q (n) = 0 for most nodes.However, if n is an entry or exit, a time-dependent piecewise constant or piecewise linear flow-profile q (n) : [t 0 , T ] ⊆ R → R is assumed to be provided.A second kind of Kirchhoff law hosted at nodes mediates boundary pressures from adjacent arcs in the following sense where p (n) (t) is the so-called node pressure.Hence, this Kirchhoff law forces adjacent boundary pressures to coalesce at all times and associates the mutual value to be the pressure of the node.Arcs represent mostly pipes but may also stand for other network elements such as active elements, like valves, regulators, and compressors.For each arc a = (n , n r ) ∈ A, the flow q(x, t) ≡ q a (x, t) and pressure p(x, t) ≡ p a (x, t) on both its ends are of particular interest and, therefore, will be abbreviated by q (t) ≡ q(0, t), q r (t) ≡ q(L, t), p (t) ≡ p(0, t), and p r (t) ≡ p(L, t).Correspondingly, the same holds true for the time derivatives, such that ṗ (t) = ṗ(0, t), ṗr (t) = ṗ(L, t), q (t) = q(0, t), and qr (t) = q(L, t).

On Spatial Discretizations of the Pipe Equations
The literature offers quite a variety of potential spatial discretizations for the Equations (1).We will pick two schemes for our presentation and numerical experiments.The first discretization in our listing originates as discretization for pipes in water networks [Huck et al., 2014] and also has been studied in-depth in [Huck, 2018] for gas networks in advance.We will refer to it as left-right-discretization The second scheme, or implicit-box -scheme, has been analyzed within [Kolb et al., 2010].This particular scheme actually combines a spatial and a time discretization altogether.Dissecting the spatial discretization relates to the trapezoidal as observed in [Streubel, 2021] and reads We will refer to (7) as the trapezoidal -discretized pipe equation.

Towards the Regulator -A Target Value Model
Regulators can down-regulate or halt gas flow through themselves, gain granular control of their surrounding neighborhood, and ultimately down the network.They can also mediate and exchange gas between areas operated at different pressure levels.Their operation is controlled by target values.In this section, we will provide a detailed explanation of the underlying target value logic, the interactions between pipes and regulators but also establish notions as well as terminology.

On the Relations of Flows and Pressures
Both regulators and compressors do manipulate the gas flow through themselves.Due to the gas dynamic modeled within neighboring pipes, the pressures will react and change as a consequence.Thus, from their direct influence on gas flow, active elements also have an indirect but immediate influence on pressures.
Elaborating on this thought, either discretization as introduced so far does approximate pressures and flows linearly along a pipe, i.e., q(x, t) Thus, we find a common precursor as space continuous system of differential equations with algebraic constraints derived from (1) as where we refer to ))] as bracket term.Now provided that B z (p) > 0, we can conclude for (8a) that net-positive inflows imply positive time derivatives of pressures, i.e., and likewise the other way around for net-positive outflows.Furthermore, B z (p) > 0 is indeed true if and only if z(p) > p∂ p z(p), which is generally true for any constant and linear model of the real gas factor and also true for the quadratic formula of Papay on its recommended operational range.Hence, we could take control over local pressures in a predictable manner by manipulating the exchanged flow between two serial pipes.This is the fundamental mechanism on which we intend to base our regulator modelings upon.Thus in the sense of this paper, we consider a regulator (sometimes also referred to as control valve) RG to be a sub-network in itself, involving 1 entry pipe P RG followed by a so-called atomic regulator RG atomic followed by 1 more exit pipe P RG r .Figure 1 depicts the described sub-network configuration.Now the atomic regulator centering between both hidden pipes is allowed to change the flow q by increasing or decreasing its resistance.In other words, it exploits our observations in (9) which roughly translate into the following expected behavioral tendencies n r p r q q q q q q r q r Figure 1: Depiction of the hidden sub-network structure (above in the picture) behind a regulator or sometimes control valve RG (below in the picture), including an atomic regulator RG atomic as well as two internal pipes P RG and P RG r .
Figure 2: Visualization of a regulator controlling the right pressure; picture created by [Kuphaldt, 2019], published under CC BY 4.02

Regulators Controlling the Gas Flow
The task of regulators in gas networks is to reduce the pressure along their orientation by decreasing the flow through the element.This flow decrease is accomplished by reducing the regulator's opening degree o ∈ [0, 1], which can be interpreted as changing the diameter of the regulator, assuming its cross-sectional area is cylindrical, see [Fügenschuh et al., 2015].If a regulator is fully opened at o = 1, it does not reduce the flow and p = p r holds.At o = 0, a regulator is fully closed, which disconnects the network.Hence, q = 0 and the pressures are decoupled.Values in between induce artificial resistance to reduce the flow and therefore create a pressure decreases in flow direction as derived in Section 3.1, i.e., p ≥ p r .
A schematic visualization of the principle operation of a regulator is depicted in Figure 2. It originates from [Kuphaldt, 2019], where more insights on various technical details and further information on other gas network elements are provided, too.

Target Value Control
For gas network operators, i.e., dispatchers, adjusting the opening degree of single elements by hand is far too involved and would require continuous observation.Instead, regulators have automated systems fed by desired ranges for the 3 local intrinsic quantities p , p r , q m .The 6 defining bounds for these ranges p in , p in , p out , p out , q , q are called target values (and sometimes also referred to as set-point values).
An overview of the 6 types of target values and their consequential influence is listed in Table 1.
If a target value imposed bound is violated by the state values (we then also refer to the target value as being violated), the element adjusts according to Table 1 to prevent or reduce the violation as much as possible.To resolve conflicts of multiple violated target values demanding opposing changes, the target value types are prioritized by ϕ(τ ) for the target values τ from top to bottom.In general, pressure target values and those reducing the opening degree and thereby the flow have higher priorities.Some target values may share their priority value, but only if they influence the element's control in the same way.q ≥ q closing minimal flow q 1 q ≤ q opening Table 1: Overview of all types of target values sorted by priority from high at the top to low at the bottom.We also listed the implied bound on the corresponding quantity of the regulator as well as the direction of change of the regulator's opening degree, which is triggered if the target value is violated.
Note that it is possible to force the regulator to be fully opened or fully closed by using certain target value combinations, for example, by setting The Regulator's Check Valve A regulator contains a built-in element that prevents flow against its orientation, i.e., q ≥ 0. This element is called a check valve or flap trap and automatically closes the regulator in case the right pressure rises above the left one, i.e. p < p r .This action has a higher priority than all the target values and therefore happens independently of these.If the regulator is not closed by the check valve, p ≥ p r holds.
Infinite Minimal Flow Target Value For most elements, the lowest priority target value q is fixed to ∞, which effectively is the same as choosing q = q ≤ ∞, in which case we may also denote q set ≡ q = q .Regulators were actually ∞ = q < q are called "flow band regulators" or simply "band regulators".They keep their current degree of openness and do not react to changes in the flow as long q stays within q ≤ q ≤ q and provided that all other target values are satisfied.This behavior would increase the overall complexity of the models used for the simulation and optimization of target values tremendously and is out of this paper's scope.
Target Value Existence Not all target values are always in use.This has to be taken into account by the models described below, for example by prohibiting their violation by setting maximal target values to ∞ and minimal target values to 0.
However, we assume that the fixed target value q always exists and in addition at least one of the closing target values.Otherwise, the element's control could only be changed in one direction.

Target Values for Compressors
In addition to regulators, also active compressors are controlled by target values.While compressors increase the pressure in flow direction, they influence their flow throughput to achieve these local pressure changes similar to regulators.For the two most common types of compressors, the flow changes are realized by controlling the compressor's rotational speed, see [Fügenschuh et al., 2015;Schmidt et al., 2017] for a more detailed explanation.Due to the similarities, the target value description for compressors is very close to the one of regulators.However, it is out of this paper's scope.

Modeling Regulators for Dynamic Simulations
For the derivation of regulator models for simulations, we will follow [Streubel, 2021].Target values are no bounds in a classical sense.Instead, they influence and guide the behavior of a regulator over time.Among ordinary differential equations, we find as a prototype model.We can utilize (10) to model the different responses of regulators to answer various kinds of imminent or present target value violations.The solution to Equation ( 10) is and always tends towards c(t).In some sense, we may describe the behavior of x as chasing c(t), or we may see x(t) as a leaky or lazy smoothing of c(t).In case c(t) = c 0 is constant, then x(t) would also converge.With the parameter α > 0, we can adjust how loose or tight x(t) will follow c(t).In the limit case ẋ/α α→∞ − −−− → 0, we find x(t) = c(t).We adapt this model to realize a regulator following its target values.We may substitute x for the flow q and replace c(t) with the target value for flows and find given q set = q = q .On top of (11), we may add the check valve q = α max(0, −q) and non-compressing behavior q = α min(p − p r , 0) and find q = α max(−q, min(p − p r , q set − q)), which is a nesting of max − min-comparisons.This nesting is necessary and represents the prioritization of the already integrated features.Also, we find sub-modelings for all remaining target values q = α min(0, p out −p r ), q = α min(0, p r − p out ), q = α min(0, p in −p ), q = α min(0, p − p in ), where we can exploit our knowledge (9) of relations between flow and adjacent pressures of neighboring pipes P RG , P RG r of Figure 1.Finally, we combine all single-aspect models within the nesting of ( 12) and end up with the full regulator model q = α max(−q, min(p − max( p in , p r ), min( p out , p ) − p r , max(q set − q, p − p in , p out −p r ))), (13) for q set = q = q .Two possible modifications of (13) are ṗr + q − ṗ = α max(−q, min(p − max( p in , p r ), min( p out , p ) − p r , max(q set − q, p − p in , p out −p r ))), (13.i) or the limit system q/α α→∞ − −−− → 0 0 = max(−q, min(p − max( p in , p r ), min( p out , p ) − p r , max(q set − q, p − p in , p out −p r ))). (13.ii)

Modeling Regulators for Discrete Optimization
To be able to optimize over the regulator's target values in the context of time discretized transient gas network operation, we aim for a model suitable for discrete optimization solvers.Our objective is to minimize the number of changes in the target values such that the control decisions induced by them are feasible to fulfill the demands of the network, which are usually given at the entries and exits.
For improved readability, we use indicator constraints of the form stating that the constraint a T x ≤ a 0 for a set of variables x = {x 1 , . . ., x i } is active if the binary variable y attains the value b.If x is bounded, these can be reformulated using linear constraints, see for example [Bonami et al., 2015].

A Basic Regulator Model
We use the following model describing the general behavior of a regulator a = ( , r) from the literature, see for example [Koch et al., 2015;Hennings et al., 2021]: The model features the mode variables m op a for an open regulator, m cl a for a closed regulator, and m ac a for an active regulator with opening degree in (0, 1), which has also been used in [Hennings et al., 2021].We add the mode m cvcl a for a regulator closed by the internal check valve in case of p ≤ p r .In this case, the target values to not influence the regulator according to Section 3.3.Note that the open mode is often replaced by a "bypass" mode in the literature, representing a bypassing network path that also allows flow in the backwards direction, which is not included in our regulator model.

Stable-Pushing Target Value Combinations
In preparation for the model construction, we determine for each target value the cases in which its implied bound has to be obeyed.While in reality, the adjustments of active elements caused by violated target values happen with a certain delay, we make the following assumption for our target value model: Assumption 5.1.The control of each active element, i.e., the opening degree adjustment for regulators, reacts to changing target values or operation point conditions immediately and with perfect precision according to the given list of target value priorities ϕ.Hence, we can assume the control is always fully adjusted to the given network situation and target values.
Note that this is equivalent to the simulation model (13.ii) assuming α → ∞.We observe the following: Observation 5.2.Since q always exists and is violated, there always is at least one violated target value.From these target values, we call the one with the highest priority ϕ the pushing target value.
Note that there might be multiple violated target values of highest priority.However, these always have the same direction, i.e., either opening or closing altogether, see Table 1.Now, we can deduce the following theorem.
Theorem 5.3.If a regulator a is not closed due to its check valve, i.e., p < p r , one of the following statements hold: 1) There is at least one target value τ , called the stable target value, with higher priority than the pushing target value as well as opposing direction in terms of change in case of violation, such that τ is equal to its corresponding operation point quantity.
2) The regulator a is fully open, i.e., in open mode, and the pushing target value is an opening target value.
3) The regulator a fully closed, i.e., in closed mode, and the pushing target value is a closing target value.
Proof.Observation 5.2 states the existence of the pushing target value as a violated target value of highest priority at any given point in time.According to the description given in Section 3.3, the regulator tries to adjust the opening degree to reduce the corresponding violation.However, according to Assumption 5.1, the opening degree is already perfectly adjusted and therefore does not change.The table shows all possible stable-pushing target value combinations, one for each line.In addition to the stable target value, the pushing target value, and the possible regulator modes, it shows the status of the bound for each target value encoded by the color of the cell: Blue stands for satisfied bounds, red for violated bounds, and orange for tight bounds in which the target value is equal to the operation point value.The text in the cells describes the necessary target value choice for each combination: The stable target value is marked by the quantity symbol combined with a "ˆ".For the violated and satisfied target values, the two possible values are "min" and "max".The value "min" encodes a target value that needs to be smaller than the state value, while "max" represents a target value that needs to be larger than the state value.The cells without coloring and filled with an asterisk do not have any relevant condition for the given row since the pushing target value has the same or a higher priority ϕ.
a) The regulator cannot adjust the opening degree in the desired direction since it is already fully opened or closed.b) There is some higher priority target value that would be violated if the opening degree would be changed by any amount in the desired direction due to the induced operation point changes.
For case a), either statement 2) or 3) holds, depending on the direction of the pushing target value's opening degree change.For case b), the higher priority target value has to have the opposite opening degree change direction compared to the pushing target value.Furthermore, the current target value must be equal to the corresponding operation point quantity since any change in the opening degree, no matter how small, would cause a violation.Hence, statement 1) is true.
Note that while the situations of 2) and 3) are coupled to a specific regulator mode, this is not the case for 1), where the two opposing target values may fix the opening degree at a position that just happens to be fully open or fully closed.
Using Theorem 5.3, we are able to compile a complete list of possible configurations in terms of stable and pushing target value combinations, which cover all possible states of a fully adjusted regulator according to Assumption 5.1.It is displayed in Table 2.Each row represents either a possible stablepushing target value combination with opposite directions or a single pushing target value without a stable target value counterpart.Since the pushing target value is by definition the highest priority target value that is violated, all bounds implied by target values having a higher priority have to be satisfied.On the other side, all constraints of target values with priority lower than or equal to the pushing target value are irrelevant for that stable-pushing combination.
Note that in the case of non-existing target value types, the corresponding lines with non-existing stable or pushing target values, as well as the columns of the missing target value types, would be removed.

A Mixed Integer Linear Formulation
We now present our mixed integer linear programming (MILP) target value model using the characterization derived in the previous section.However, in contrast to the overview presented in Table 2, we will only allow the usage of the active mode for all target value combinations featuring a stable target value.This is based on the fact that we only use ≤ and ≥ relations to represent violated bounds since strict inequalities cannot be used in the MILP context.As a result we formulate the following theorem: Theorem 5.4.Given a stable-pushing target value combination c 1 as well as a set of target values and state values satisfying the constraints implied by c 1 .Then these values also satisfy the constraints for an open mode combination c 2 and a closed mode combination c 3 .The two different pushing target values of c 2 and c 3 are equal to the stable and the pushing target value of c 1 .
Proof.Since we only use ≤ and ≥ to represent violated bound constraints, the stable target value of c 1 can be alternatively interpreted as a violation as well as a fulfillment of the corresponding target value bound.Hence, the values satisfy the constraints of that non-active mode combination c 4 , which has the stable target value of c 1 as pushing target value, i.e., that violated target value of highest priority.The same is true for that non-active mode combination c 5 , which has the pushing target value of c 1 as pushing target value without having an explicit stable target value.Since the stable and pushing target values of c 1 have opposite directions, one is a closing target value, and one is an opening one.Hence, one of the combinations c 4 and c 5 has the open mode and is the aforementioned combination c 2 , and the other has the closed mode and is c 3 .
As a result of Theorem 5.4, we assume that the regulator is in active mode iff it has a stable target value since all the open and closed mode cases with stable target value are covered by corresponding combinations without a stable target value.This enables us to formulate the following MILP model for a target value control of a single regulator a = ( , r) at time t.As the model extends the basic regulator model ( 14), we will only state the additional variables and constraints.First, we introduce some sets: All target values existing for regulator a Θ ac := {( q , q ), ( p in , p in ), ( p in , p out ), All stable-pushing target value combinations for active ( p in , q ), ( p in , q ), ( p out , p in ), ( p out , p out ), ( p out , q ), ( p out , q )} Θ op := {(−, p in ), (−, p out ), (−, q )} All stable- For the following set of constraints, we introduce the function o mapping from a given type of target value to the corresponding point of operation: In addition, we use ϕ as a function from target value types to the priorities given in Table 1 above.Also, we introduce the parameter ε, which represents a relative tolerance used in the equality constraint of the stable target value.As explained at the beginning of Section 5.2, the adjustment of the regulator's control to changing target values or operating point conditions happens with a certain delay in reality.Hence, there usually is an offset between the used target values and the current operation point.To better reflect this behavior and prevent too many target value changes caused by minor adjustments of the stable target value, we introduced the relative ε value.
The constraints for our target value model are then given as Equation ( 15a) ensures the choice of exactly one stable-pushing target value combination fitting to the current mode.Note that for the check-valve-closed mode, no combination is chosen.Constraints (15b), (15c), and (15d) establish the consequences of the choice of a stable-pushing target value combination and are based on the illustration of Table 2. Structure-wise these constraints collect for a given target value type τ all stable-pushing target value combinations in which this type of target value is the stable one for (15b), is violated for (15c), or is satisfied for (15d).If one of these combinations is chosen, the corresponding relations between the target value and the corresponding operation point are enforced.Finally, constraint (15e) defines the target value change variables, forcing the actual target value variable τ x,a,t to keep their value if the change variable does not indicate a change, i.e., is zero.
As objective function, we minimize the number of target value changes: min a,t x∈Γa δ x,a,t .
Note that [Pfetsch et al., 2014;Schmidt et al., 2015;Koch et al., 2015] introduced a model for regulators without remote access that have a set-point control for the maximum right pressure with value p set .Depending on whether the right pressure is larger than, equal to, or less than p set , the mode is closed, active, or open, respectively.In addition, these regulators also feature a check valve behavior, which in their case allows flow against the regulator's orientation (bypass instead of open mode).Except for this back flow, the described behavior can as well be achieved with our formulation, for example, by setting τ p out = p set , τ p in = 0, and τ p out = ∞.The target value p in is always satisfied and p out is always violated.Hence, we are either in active mode and have a stable p out , with p r equal to the set-point value, or p out is violated, leading to a closed regulator, or p r is below p out and satisfied, leading to the open mode.This shows that our optimization model covers regulators with and without remote access.

Numerical Evaluation
To show and verify the accuracy of our models, we present in this section two experiments using the target value control on a single regulator.It is embedded in a simple path network between two pipes, each with L = 10 km, D = 0.9 m, r = 0.012 mm, see Figure 3.
For both scenarios, we determine the solution using the dynamic simulation model presented in Section 4 as well as the optimization model presented in Section 5.The reference solution is determined using the SIMONE simulator3 , which has established itself as a de facto standard due to its wide use in industry.We will refer to solutions computed with SIMONE as industry simulations as a convenient shorthand.The initial steady state of the scenarios features a flow of 10 kg/s on all arcs as well as inflow at the entry and outflow at the exit.The pressure ranges from 50 bar at n in to 49.992 bar at n out , the regulator is fully open, the constant gas temperature is set to 283.15 K, and the gas mixture is, for the sake of simplicity, assumed to be pure methane.
As pipe model for the simulation, we use the left-right-discretization (6), equipped with the linear AGA formula for the real gas factor.Note that we chose α = 1.0e3 for the simulation model.
For the optimization, we take the linearized pipe model (1b") discretized with the implicit box scheme.Furthermore, we use a constant real gas factor z, which is determined for each node using the formula of Papay and computed from the initial state.The linearized model is not reliable in the case of volatile flow conditions but works fine for the following scenarios featuring similar flows for most of the time.As a benefit, the resulting overall model is linear, which enables us to use a MILP solver for the optimization, in our case Gurobi [Gurobi Optimization, 2020].Furthermore, we picked ε = 0 for the optimization model.

First Scenario
For the first scenario, we use a set of target values obtained from our project partner, the gas network operator Open Grid Europe GmbH (OGE).These target values should lead to different opening degree and throughput changes of the regulator while assuming a steady inflow and outflow at the boundaries of 10 kg/s over the whole time horizon of 12 h.This scenario is a stress test.It has been designed to provoke relevant interferences from the regulator involving almost each available target value.An overview of the different switching actions of the target values is provided in Table 3.For the optimization, the TV Time as hh:mm Value p in 00:00 48.0 bar p out 00:00 55.0 bar p in 00:00 100.0 bar p out 00:00 40.0 bar q 00:00 9.0 kg/s q 01:00 15.0 kg/s q 02:00 6.0 kg/s q 02:30 10.At the initial start, all the pressure target values are fulfilled, but q = 9 kg/s is violated.Hence, the regulator closed until q = q .At 01:00, q is set to 15 kg/s.Hence, the regulator fully opens again, the pressures equalize, and we go back to the initial steady state.Next, q is set to 6 kg/s for 30 minutes and then goes to 10 kg/s.The regulator enforces corresponding values of q, which first create a pressure imbalance by reducing the flow and then stabilizes this imbalance.At 3:30, we start using the pressure target values by setting p out to 47.0 bar, which is violated.Since p out has the highest priority, the regulator closes until p r = 47.0 bar and then keeps this value, i.e., q = 10 kg/s is established again.After setting back p out , the situation does not change since q still enforces q = 10 kg/s.Now, p in = 55.0 bar is set, which again is violated and of high priority.Therefore, we further increase the pressure imbalance until reaching p = 55.0 bar and keep this situation stable, even though p in is reduced to 53 bar again shortly after.The next new setting is p out to 46.0 bar at 06:30, which is a violation of the target value.Since p in and p out are currently fulfilled, this forces the regulator to open until p r = 46.0bar is reached.Here, we stabilize the flow at q = 10 kg/s again.The reduction of q to 6 kg/s does not change this since p out is of higher priority than q .As final changes, we first set p out to 46.5 bar and then to 47.5 bar.The first change shifts the stable situation to p r = 46.5 bar.However, we do not reach p r = 47.5 bar afterwards since this would violate p in = 53 bar.Hence, we stabilize at p = 53.0bar and keep this situation until the end of the time horizon.
All three solutions are very similar to each other in terms of the development of the regulator's operation point over time.Hence both derived models, which are the simulation model and the optimization model, do catch the regulator's behavior compared to the industrial simulation.As an overview, we listed below the maximal relative error and the relative error at the end of the time horizon for the pressure values, where we always choose the bigger of the two errors given for the left and right pressure, as well as the flow values.The values show that the differences in the pressures are, in general, smaller than those in the flows.This is due to the fact that each of the three methods has a different internal implementation of the reaction of the regulator's opening degree to changing target values conditions.Hence, even though the time discretization is the same for all approaches, the flow values are different in the unstable periods shortly after a target value change leading to a new violation.However, these differences only occur for brief periods of time.We also note that the differences between the two simulation procedures are, in general, smaller than those between the optimization and one of the simulation runs.This is no surprise, as we assume an idealized target value behavior for the optimization model, see Assumption 5.1.However, we are glad to see how close the solution of the idealized model is to both simulations.

Second Scenario
In the second scenario, we tested both models in a different setup.Here, we determine the regulator's target values using the optimization model and verify them using the simulation model.Since this task is quite demanding in a volatile environment, even for only a single regulator, we discretized the time horizon for the optimization model into steps of 15 minutes, leading to 48 time steps in total, whereas the simulation still solves in 3 min time step resolution.
To ensure and provoke the necessity of target value adjustments, we prescribed not only the inflow at the entry and outflow at the exit but also the pressure values over time at both nodes.We took these pressure values of the OPT solution from scenario 1 with time step resolution of 15 min.This enables us to compare the number of target value changes between the OPT solutions of scenarios 1 and 2 and provide an estimate of the optimization potential.
There is a fundamental difference in the interpretation of value updates of target values between the optimization and the simulation model.A simulation starts with the control adjustment at some point in time, let it be t, and proceeds over the following time steps with the adjustment to the new target value(s).In contrast to this and due to Assumption 5.1, the optimization model assumes a control which has already adjusted to the new target value directly at time t.Hence, it basically assumes that the target values have been changed at some point within the 15 minute time step preceding time t.Here, we have a degree of freedom and could eliminate this by a subsequent event detection approach to narrow down the event-interval of the switch to fit the time-resolution of the following simulation.Instead, we use a heuristic and will communicate changes in the target values that occurred at time t in the optimization model to both simulation processes with time t−∆ t /2, where ∆ t is the time resolution of the optimization model, i.e., ∆ t = 15 min.This minimizes the maximum possible error between the used and the actual change time.Note that we round 7.5 min to the integer value of 8 min here.
The target values determined by the optimization are listed in Table 4.The comparisons of the different solutions over time are shown in Figure 7 for the simulation (SIM) and the industry simulation (IND), in Figure 8 for optimized solution (OPT) and the industry simulation, and in Figure 9 for the optimized and simulated solution.
When comparing the listed target value changes in Tables 3 and 4, we observe that we only needed 8 target value changes instead of the previous 11 to archive the regulator control.In the optimized solution, we do not use the changes of p out at 04:30 and q at 06:30 as they do not affect the regulator's control.In addition, the control change at 07:30, which was a consequence of changing p in at 05:30 and p out at 07:30, could be realized by a single change of p out at 07:30.Similar to this last case, some of the other control changes are now triggered by different target value changes.Another example would be the change at 05:00, which was previously caused by setting p in to 55 bar and is now triggered by setting p out to 45 bar.
When looking at the solutions over time, the development is again quite similar, which can also be observed from the table below.It shows the maximal relative error and the relative error at the end of the time horizon for the pressure values, where we always choose the bigger of the two errors given for the left and right pressure, as well as the flow values.max p error max q error end p error end q error SIM vs IND 0. 00:00 00:00 48.0 bar p out 00:00 00:00 55.0 bar p in 00:00 00:00 100.0 bar p out 00:00 00:00 40.0 bar q 00:00 00:00 9.0 kg/s q 01:00 00:52 13.0 kg/s q 02:00 01:52 6.0 kg/s p out 02:  The error values of the comparison between our simulation and the one of the industry software are consistent with those of the first scenario.This is not surprising, as only the simulated instance changes while the general setup for both stays the same.When comparing one of the two and the optimized solution, we clearly see the difference to the values of the first scenario.Especially the maximum flow error, but also the maximal and overall errors in the pressure have increased significantly.This is an expected change, as the optimization now uses a more coarse discretization.Nevertheless, the pressure differences are still very small in general and not too far away from the values of either simulation.
Regarding the errors in the flow, the main difference occurs in the very first minutes of the scenario.Since the initial target values are violated, the control adjusts immediately, which happens during the  first time step for all of the three approaches.Due to the different time discretizations, the adjustment is 5 times faster for the simulations than for the optimization.However, the overall error regarding the flow values is surprisingly small and consistent to the values of the first scenario.We summarize that both simulations yield very similar results, while the solution of the target value optimization may produce different values during the times of adjusting the control to new target values, but gives very similar results regarding the stabilized control situations.

Outlook
In this article, we introduced models for regulators in gas networks and embedded target values for their control.The target values form a system of 6 simultaneously active and hence competing set-points, which enables the regulating element to react to changing conditions in the surrounding network without the need of continuous human observation or manual interactions.Our description lives up to expectations and needs compared to features found in industry-standard simulators.
To capture the behavior of target value controlled regulators, we introduced two different models.The first one is a model suitable for dynamic or transient simulation.Here, the target values are assumed to be piecewise linear functions, and their competing logic is encoded in terms of nested max − min-comparisons.A second model was designed to be used in discrete optimization problems and determines an optimized set of target values changing as little as possible while enabling the demanded network control.It is based on a characterization of the target value behavior as a set of specific element states, the stablepushing target value combinations.These assume a regulator reacting to changes in its environment or of its target values immediately and perfectly precise.Depending on the chosen combination, different constraints are active that relate the target values to the values of the element's point of operation.
Both models have been evaluated in comparison to a commercial simulator, actively used in industry.We showed that the behavior of our simulation model matches the expectations in that the simulated solutions stick close in terms of the resulting pressure and flow values.The relatively higher differences in the flow values are short-lived and their integral equaling the total amount of transported gas is negligible.These differences in the flow values represent slightly different timings and adjustment speeds of regulator actions.Possible strategies to address this could be the introduction of delay or slope limiter.However, this is more a question of technical specification rather than mathematical consideration.For the optimization model, we could show that it works very well when using the same time discretization as for the simulation together with already fixed target values.Here, the differences to both simulated results are just slightly larger than both simulations among each other, even though we used a model based on more simplifying assumptions.In a more realistic scenario, in which the target values are subject to optimization, we need to use fewer and larger time steps.As a result, the relative differences to the simulated results have increased.However, they still stay very close and bounded, especially the components representing the pressures.This indicates that the optimization model is able to reproduce the target value behavior reasonably well.
There are multiple paths to pursue for future research.Firstly, we could replace or add other kinds of target values, e.g., bounding values for the ratio or difference of in-and outlet pressures.Secondly, the regulator model for simulation could be converted into a model for band regulator if q ≤ q by replacing (11) with q = α[min(0, q −q) + max(0, q −q)], just before adding the remaining target value mechanisms.Thirdly, we may apply our techniques to embed target-value-based controls into compressor unit, group, and station modeling.In that regard, we have already discussed it briefly in Section 3.4, but it would need to be elaborated and established further in full detail.Furthermore, the optimization model can probably be improved to enhance the performance when combined with state-of-the-art optimization frameworks.While the model captures the target value behavior very well already, it is quite challenging to solve for one single element already.The overall goal is to make it usable in a transient gas network operation model using multiple regulators and hundreds of other network elements.Finally, the concept of a target-value-based element control can be applied to active elements in other fluid-based flow networks.As an example, we mention pumps in water networks, which from a macroscopic viewpoint similarly control their network as compressors do for gas networks.

Figure 3 :
Figure 3: Benchmark network, consisting of one regulator arc and two pipes surrounding the regulator.
(OGE)  in our benchmark scenario.Each row represents a target value change of the given type of target value to the given value at the given time.The time is specified as hours and minutes past the initial state.simulation, as well as the reference industry simulation, we use a time discretization of 3 min.Note that the target values are fixed for the optimization here, which effectively transforms the problem into a plain feasibility problem.The results for the first scenario are shown in Figure4comparing the simulated solution (SIM) with the one computed by the industry simulator (IND), Figure5compares the optimized solution (OPT) with the one from the industry simulator, and Figure 6 compares the optimized and simulated solution.Each figure shows the development of the left and right pressure values and target values of the regulator in the

Figure 4 :
Figure 4: Comparison of the simulated solution (SIM) with the reference simulation from the industry simulator (IND) for the first scenario with predefined target values based on Table 3.The pictures show the left and right pressure values and target values of the regulator, the flow values and target values, and the difference of both solutions in the three quantities integrated over time and normalized by the quantity integral of the IND values.

Figure 5 :
Figure 5: Comparison of the optimized solution (OPT) with the reference simulation from the industry simulator (IND) for the first scenario with predefined target values based on Table 3.The pictures show the left and right pressure values and target values of the regulator, the flow values and target values, and the difference of both solutions in the three quantities integrated over time and normalized by the quantity integral of the IND values.
The target values determined by the optimization model.Each row represents a target value change of the given type of target value to the given value at the given time.The time is specified as hours and minutes past the initial state.The OPT time is the one determined by the optimization, and the SIM time is the one communicated to both simulation processes.

Figure 7 :
Figure 7: Comparison of the simulation (SIM) with the reference simulation from the industry simulator (IND) for the second scenario with optimized target values, see Table 4.The pictures show the left and right pressure values and target values of the regulator, the flow values and target values, and the difference of both solutions in the three quantities integrated over time and normalized by the quantity integral of the IND values.

Figure 8 :
Figure 8: Comparison of the optimized solution (OPT) with the reference simulation from the industry simulator (IND) for the second scenario with optimized target values, see Table 4.The pictures show the left and right pressure values and target values of the regulator, the flow values and target values, and the difference of both solutions in the three quantities integrated over time and normalized by the quantity integral of the IND values.

Figure 9 :
Figure 9: Comparison of the optimized solution (OPT) with the reference simulated solution (SIM) for the second scenario with optimized target values, see Table 4.The pictures show the left and right pressure values and target values of the regulator, the flow values and target values, and the difference of both solutions in the three quantities integrated over time and normalized by the quantity integral of the SIM values.

Table 2 :
The τ variables represent the target values chosen by the model.They can be set inside the variable bounds τx,a and τx,a for each target value x.These are given for each individual regulator a and can considerably restrict the feasible operating range of the regulator.If, for example, the minimal left pressure target value p in has a lower bound value of 40 bar, then the regulator has to operate at a left pressure value of at least 40 bar in active mode and open mode, since according to Table2, the target value p in can only be violated for the closed mode or if the regulator is closed due to the check valve.The variables θ represents the stable-pushing target value combination chosen for each regulator a. Finally, the δ variables indicate a change in the target value variables τ between the values of the previous and the current time point.
pushing target value combinations for open Θ cl := {(−, p in ), (−, p out ), (−, q )} a All existing combinations for regulator a Now, we can specify the new variables as τ x,a,t ∈ [ τx,a , τx,a ]

Table 3 :
The target values specified by Open Grid Europe GmbH max p error max q error end p error end q error