Chance-constrained optimal inflow control in hyperbolic supply systems with uncertain demand
CChance-constrained optimal inflow control inhyperbolic supply systems with uncertain demand
Simone G¨ottlich (cid:42) , Oliver Kolb (cid:42) , Kerstin Lux (cid:42)
January 13, 2020
Abstract
In this paper, we address the task of setting up an optimal production plan taking into ac-count an uncertain demand. The energy system is represented by a system of hyperbolic partialdi ff erential equations (PDEs) and the uncertain demand stream is captured by an Ornstein-Uhlenbeck process. We determine the optimal inflow depending on the producer’s risk prefer-ences. The resulting output is intended to optimally match the stochastic demand for the givenrisk criteria. We use uncertainty quantification for an adaptation to di ff erent levels of riskaversion. More precisely, we use two types of chance constraints to formulate the requirementof demand satisfaction at a prescribed probability level. In a numerical analysis, we analyzethe chance-constrained optimization problem for the Telegrapher’s equation and a real-worldcoupled gas-to-power network. AMS subject classifications.
Keywords.
Stochastic optimal control, chance constraints, Ornstein-Uhlenbeck process, hyper-bolic conservation law
In recent years, significant attention has been paid to the energy market. On the one hand, thisis due to climate protection policies. On the other hand, the decision of the German governmenton the nuclear phase-out by the end of 2022 (cid:42) will cause significant changes in the energy sector.Those changes are to a large extent triggered by the specification to have 65% of the gross elec-tricity consumption generated out of renewable energy sources by the end of 2030. One majorchallenge is the handling of uncertainty in the renewable energy production. It heavily dependson the weather conditions and is subject to large fluctuations. Those fluctuations can heavilya ff ect the grid operation or even lead to outages ([5]).Another source of uncertainty in the power grid operation are power demands (see [50]). Thismakes it di ffi cult to guarantee a su ffi cient supply thereby influencing reliability and profitabilityin power system operation (see [19]). It is crucial for all players in the electricity sector to copewith stochastic demand fluctuations. In the context of gas-to-power, those fluctuations mightcarry over to the gas network being coupled to the electricity grid. The amount of gas convertedto power and feeded into the electricity system depends on the power demand and thereforeinherits its stochasticity. Furthermore, the gas demand itself is uncertain (see e.g. [27]).There are several ways to formulate optimization problems in the presence of uncertainty. As in[21], we distinguish between robust optimization and chance constrained optimization (CCOPT) . A (cid:42) University of Mannheim, Department of Mathematics, 68131 Mannheim, Germany ([email protected],[email protected]). (cid:42) , last checked:4 th of April, 2019 a r X i v : . [ m a t h . O C ] J a n obust formulation is appropriate if the distribution of the uncertain parameter in the system isunknown because it cannot be measured or observed (no historical data), whereas a probabilisticprogramming (CCOPT) is suitable in the presence of historical observations where a distributionof the unknown quantity can be derived (see [21]).In this paper, we focus on uncertain demands. In this context, it is reasonable to assume accessto historical demand data and we therefore focus on CCOPT for the remainder of the presentmanuscript. Chance constraints have been introduced in 1958 in [7] by Charnes, Cooper andSymonds in the context of production planning, and further developed in [6]. A lot of workhas been done in CCOPT. It is worthwhile to particular mention the contributions of Pr´ekopa.His monograph on stochastic programming from 1995 [38] nowadays still serves as a standardreference. A general overview on properties, solution methods, and fields of application such ashydro reservoir management ([2, 46]), optimal power flow ([5, 50]), energy management ([45]),and portfolio optimization ([39]) can be found in [19].In this work, we are interested in a constrained optimization problem composed of three con-straints: a stochastic di ff erential equation (SDE) that models the uncertain demand, a system ofhyperbolic balance laws to describe the energy system (electricity or gas), and a chance constraintensuring demand satisfaction with a certain probability (risk level). Chance constraints in thecontext of PDE-constrained optimization are also considered in [34]. However, the requirementof a continuously Fr´echet-di ff erentiable solution of the PDE rules out most hyperbolic PDEs dueto the possibility of discontinuous solutions even for continuously Fr´echet-di ff erentiable initialdata.In [15], they apply their results on semi-continuity, convexity, and stability in an infinite-dimen-sional setting to a simple PDE-constrained control problem. Again, the PDE is not of hyperbolicnature.There already exist a few investigations of CCOPT for systems of hyperbolic nature such as gasnetworks. Accounting for the stochasticity of demand, [25] presents a method to compute theprobability of feasible loads in a stationary, passive gas network. Note that the steady state as-sumption on the network entails constraints formulated in terms of purely algebraic equationsinstead of full hyperbolic dynamics. In [21], they extend the aforementioned work by also allow-ing for uncertainty with respect to the roughness coe ffi cient. They analyze the question of themaximal uncertainty allowed such that random loads can be satisfied at a prescribed probabilitylevel within a stationary, passive gas network. The feasibility of random loads, again in a station-ary gas model, this time with active elements in terms of compressor stations is considered in[26]. In [1], they also account for the transient case. However, they do not take into account thefull hyperbolic dynamic but the linear wave equation instead.In contrast to the above-mentioned contributions, here, we consider the full isentropic gas dy-namics. Note however that the full stochasticity of the demand in the power-to-gas setting pre-sented in Section 4.3 does not enter the nonlinear dynamics of the gas flow directly but is coupledto it via the objective function.Another important di ff erence to the aforementioned contributions is the modelling of the uncer-tain demand. In many cases, (truncated) Gaussian random vectors are used as proposed in [27]and used in for example [1, 21, 26]. Using stochastic processes to model the demand enables tocapture time dynamics. Sometimes the demand is modeled by a discrete time stochastic processas e.g. in [2]. In [46], they use a continuous time stochastic process additively decoupled in adeterministic trend and a causal process generated by Gaussian innovations. In contrast to that,we consider a dynamic demand model for the random loads via a continuous time stochastic pro-cess, the so-called Ornstein-Uhlenbeck process, where the deterministic trend and the stochasticevolution of the process are coupled. This is also a popular choice to model demand in vari-ous fields (see e.g. [3] in the context of electricity). In [37], and [48], they use a multivariateOrnstein-Uhlenbeck process from the supply-side point of view. They model uncertain injections2nto the power network and assess the probability of outages. A desirable feature for both supplyand demand modeling is the mean reverting property. The process is always attracted to a cer-tain predefined mean level, which may result from historical supply/demand data. In contrast tothose contributions, we use a time-dependent mean level, which enables us to depict seasonal anddaily patterns appearing in historical data also in our mean-demand level (no longer assumed tobe constant) and set up an optimization framework based thereon.The paper is organized as follows: in Section 2, we introduce the stochastic optimal control set-ting. We define the energy system as well as the stochastic demand process, and introduce theconsidered cost functional. We distinguish between a single chance constraint (SCC) as well as ajoint chance constraint (JCC). In Section 3, we consider a deterministic reformulation of both thecost functional and the SCC, and present a stochastic reformulation of the JCC, which allows toincorporate it in our numerical framework. In Section 4, we validate the numerical frameworkfor a simple case of a hyperbolic supply system, i.e. the scalar linear advection with source term.It represents a suitable test case setting as we are able to derive the corresponding analytical so-lution in case of an SCC. We then apply our numerical routine to the Telegrapher’s equations,a linear system of hyperbolic balance laws. We conclude with a numerical investigation of anonlinear system of hyperbolic balance laws in terms of a real-world gas-to-power application. In this section, we set up the mathematical framework for the task of finding an optimal injectionplan taking into account an uncertain demand for the time period from t = 0 to final time T .This has been done for a linear transport equation in [23]. Here, we use the stochastic optimalcontrol framework set up in [23], and extend it to more complex supply dynamics on a network.We model a general energy system by a system of hyperbolic balance laws on a network, and thestochastic demand is described by an Ornstein-Uhlenbeck process (OUP). We consider di ff erent types of 2-dimensional energy systems. The network is modeled by a finite,connected, directed graph G = ( V , E ) with a non-empty vertex (node) set V and a non-empty set ofedges E . For v ∈ V , we define the set of all incoming edges by δ − ( v ) = { e ∈ E : e = ( · , v ) } , and the setof all outgoing edges by δ + ( v ) = { e ∈ E : e = ( v, · ) } . In the sequel, we identify an edge e = ( v in , v out )by the interval [ a e , b e ], where a e denotes the starting point of the edge, and b e its end point. Wefurther define the set of inflow vertices by V in = { v ∈ V : δ − ( v ) = ∅} , and the set of outflow verticesby V out = { v ∈ V : δ + ( v ) = ∅} . As a simplification, here, we restrict our network to |V in | = |V out | = 1.Note however that an extension to several inflow nodes ( |V in | > |V out | > u ( t ) acts on the vertex v in ∈ V in , and the demand Y t is realized at the vertex v d ∈ V out . We require b e = L for all e ∈ δ − ( v d ).The dynamics of the energy system on one edge e are given by the following hyperbolic balancelaw with initial condition (IC) and boundary conditions (BCs) ∂ t q e + ∂ x f e ( q e ) = s ( q e ) , x ∈ [ a e , b e ] , t ∈ [0 , T ] q e ( x,
0) = q e ( x ) , Γ ea ( q e ( a e , t )) = γ ea ( t ) , Γ eb ( q e ( b e , t )) = γ eb ( t ) . (1)Thereby, f : R → R is a given flux function and s : R → R the source term. ρ e : [ a e , b e ] → R describes the initial state of the system on edge e . The functions Γ ea/b : R → R m l/r enable to pre-scribe a certain evaluation of a certain number of components of the density at the boundary.3 in v v v v v d S u ( t ) Demand Y t Inflow u ( t ) e e e e e e e Figure 1: Example of an energy network with uncertain demand at v d Thereby, 0 ≤ m l , m r ≤ m l/r = 0 has to be interpreted in a way that no left/right boundary condition isprescribed. For a scalar conservation law, choosing m l = 1, m r = 0, and Γ ea = f corresponds toprescribing only the flow at the left boundary.The boundary conditions (BCs) themselves are given in terms of the functions γ ea : [0 , T ] → R m l at the left boundary a e of edge e , and γ eb : [0 , T ] → R m r at the right boundary b e of edge e . Thenumbers of BCs m l/r , and the functions γ ea/b need to be chosen carefully such that they are consis-tent with the characteristics of the conservation law. This will be further specified below for eachsetting 4.2-4.3.Moreover, for the non-degenerate network case, i.e. ∃ v ∈ V : | δ − ( v ) | + | δ + ( v ) | ≥
2, suitable couplingconditions c v : R × ( | δ + ( v ) | + | δ − ( v ) | ) → R m c for each vertex need to be imposed. Thereby, m c denotesthe number of coupling conditions. This will be made explicit in Subsection 4.3.For e ∈ δ + ( v in ), the function γ ea ( t ) depends on the inflow control u ( t ).To simplify notation, we do not explicitly write down the dependence of the supply on the densityat the end points of all ingoing edges and denote the supply at v d at time t simply by S u ( t ) meaning S u ( t ) = S u (cid:17) e ∈ δ − ( v d ) q e ( b e , t ) . Discretization scheme
For the numerical investigation in Section 4, the considered hyperbolic energy systems need anappropriate discretization scheme. Motivated by our real-world example, we choose an implicitbox scheme (IBOX) [32] for all considered scenarios. For a general system of balance laws (on anyedge) ∂ t q + ∂ x f ( q ) = s ( q ) , the considered scheme reads Q n +1 j − + Q n +1 j Q nj − + Q nj − ∆ t ∆ x (cid:16) f ( Q n +1 j ) − f ( Q n +1 j − ) (cid:17) + ∆ t s ( Q n +1 j ) + s ( Q n +1 j − )2 . (2)Here, ∆ t and ∆ x are the temporal and spatial mesh size, respectively, and the numerical approxi-mation is thought in the following sense: Q nj ≈ q ( x, t ) for x ∈ X j = (cid:104) a + ( j − ) ∆ x, a + ( j + ) ∆ x (cid:17) ∩ (cid:104) a, b (cid:105) , t ∈ I i = (cid:104) i ∆ t, ( i + 1) ∆ t (cid:17) . (3)4o avoid undesired boundary e ff ects, the discretization of the initial condition on bounded do-mains is done pointwise, i.e., Q nj = q ( a + j ∆ x ) . (4)As remarked in [32], for a discretization x l < x l +1 < · · · < x r − < x r , we obtain r − l equations for r − l +1 variables (in the scalar case). This entails the need to prescribe BCs at exactly one boundaryspecified by the characteristic direction. This also explains the above mentioned assumption ofno change in the signature of the characteristic directions on the considered domain. The discreteversion of the BCs is given by Γ a ( Q i ) = γ a ( t i ) , Γ b ( Q iN ∆ x ) = γ b ( t i ) ∀ i ∈ { , · · · , N ∆ t } . (5)Note that the implicit box scheme has to obey an inverse CFL condition [32], which is beneficialfor problems with large characteristic speeds whereas the solution is merely quasi-stationary.This is usually the case for daily operation tasks in gas networks and therefore motivates thechoice within this work. Uncertain Demand
The uncertain demand stream is modeled by an Ornstein-Uhlenbeck process (OUP), which is theunique strong solution of the following stochastic di ff erential equation dY t = κ ( µ ( t ) − Y t ) dt + σ dW t , Y t = y t (6)on the probability space ( Ω , F , P ). W t is a given one-dimensional Brownian motion on the sameprobability space, and y describes the demand at time t = 0. The constants σ > κ > µ ( t ), called the mean demand level. This is due to the sign of the drift term κ ( µ ( t ) − Y t ), whichensures that being above (below) the mean demand level, the process experiences a reversionback to it. By the OUP (6), we are able to capture daily or seasonal patterns in the demand. Froma mathematical point of view, this process has some nice analytical properties. For example, wecan derive the solution to equation (6) explicitly via the Itˆo formula. It reads as Y t = y t e − κ ( t − t ) + κ (cid:90) tt e − κ ( t − s ) µ ( s ) ds + σ (cid:90) tt e − κ ( t − s ) dW s . (7)Moreover, it is possible to derive its distribution explicitly as Y t ∼ N y t e − κ ( t − t ) + κ t (cid:90) t e − κ ( t − s ) µ ( s ) ds (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) m OUP ( t ) , σ κ (cid:16) − e − κ ( t − t ) (cid:17)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) v OUP ( t ) . (8)For further details on the demand process, we refer the reader to [23]. Requiring demand satisfaction for every realization of the demand process might be too restric-tive as, in some cases, it might lead to an infeasible optimization problem ([41, p. 5]).5ne possibility to overcome problems of infeasibility and to reduce the average undersupply is tointroduce an undersupply penalty term in the cost function. The e ff ect of an undersupply penaltyon the optimal output has been analyzed in [24]. A comparison of di ff erent types of undersupplycan be found in [36].Another approach is to guarantee with a certain probability that there is no undersupply withina prescribed time interval I CC ⊂ [ t ∗ , T ], where t ∗ is the first time that a supply is realized at v d . Mathematically this is formulated in terms of a chance constraint (CC). One possibilty is torequire at each point in time that the probability of a demand satisfaction is at least equal to oneminus a given risk level θ (see (9a)). This results in a so called single chance constraint (SCC).Another possibility is a joint chance constraint (JCC), that is, we require that the probabilityof a demand satisfaction is at least equal to one minus a given risk level θ on a whole interval simultaneously (see (9b)). P (cid:16) Y t ≤ S u ( t ) (cid:17) ≥ − θ ∀ t ∈ I CC , (9a) P (cid:16) Y t ≤ S u ( t ) ∀ t ∈ I CC (cid:17) ≥ − θ. (9b) Having formulated all the optimization constraints, we now address the objective function. Weaim at minimizing the arising costs. We construct our cost function out of several components.One component consists of deterministic costs C u ( t ) such as operating costs for a gas compressor( C S u realized at the demand vertex v d . We measure thetracking type costs in terms of the expected quadratic deviation between the demand and supply( C C R ).Putting the deterministic costs C
1, the tracking costs C
2, the undersupply costs C
3, and theexcess revenue C OF ( Y t , t , y t , S u ( t ) ) = w · C u ( t ) (cid:124) (cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32) (cid:125) C + w · E (cid:20)(cid:16) S u ( t ) − Y t (cid:17) | Y t = y t (cid:21)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) C − w · E (cid:104)(cid:16) S u ( t ) − Y t (cid:17) − | Y t = y t (cid:105)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) C − w · E (cid:104)(cid:16) S u ( t ) − Y t (cid:17) + | Y t = y t (cid:105)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) R , (10)where (cid:16) S u ( t ) − Y t (cid:17) − = (cid:40) S u ( t ) − Y t if S u ( t ) < Y t S u ( t ) ≥ Y t All together, we come up with the stochastic optimal control (SOC) problemmin u ∈U ad (cid:90) Tt ∗ OF ( Y t ; t ; y t ; S u ( t ) ) dt subject to (1) , (6) , and (9) . (11)To determine the optimal control, we need to specify measurability assumptions on the controlby defining the space of admissible controls U ad . U ad = { u : [ t , T ] → R | u ∈ L ([ t , T ]) , u ( t ) ≥ , and u ( t ) is F t –predictable for t ∈ [ t , T ] } . (12)6esults for this control method and two other control methods with an objective function of puretracking type ( w = w = w = 0, w = 1) for the linear advection equation without imposing CCscan be found in [23]. Having set up the SOC problem (11), the question of how to solve this minimization problemnaturally arises. One way is to trace back the SOC problem (11) to a deterministic setting, inwhich we can apply well-known methods from deterministic PDE-constrained optimization suchas adjoint calculus. In order to do so, in 3.1, we analytcially treat both types of CCs presented inSubsection 2.2. The deterministic expression of the objective function introduced in 2.3 is derivedin 3.2.
The reformulation of the CC heavily depends on the type of CC. Whereas the SCC (9a) can bereformulated via quantiles of a normal distribution, the reformulation of the JCC (9b) is notobvious. Therefore, we need to treat the di ff erent types of CCs separately. Single chance constraint
For the SCCs (9a), we use a quantile-based reformulation as mentioned in [2, p. 580]. By usingthe known distribution (8) of the OUP, this results in the deterministic state constraints S u ( t ) ≥ m OUP ( t ) + v OUP ( t ) Φ − (1 − θ ) ∀ t ∈ I CC , (13)where m OUP ( t ) = y t e − κ ( t − t ) + κ t (cid:82) t e − κ ( t − s ) µ ( s ) ds , v OUP ( t ) = σ t (cid:82) t e − κ ( t − s ) ds , and Φ is the standardnormal cumulative distribution function. Joint chance constraint
JCCs (9b) are mathematically by far more involved than SCCs (see [19, p. 1213]). No longerconsidering the constraint pointwise in time, we now have to deal with a joint probability dis-tribution. Unfortunately, we can no longer make use of the deterministic reformulation as stateconstraint as for (9a).As the integration of the JCC into the optimization framework is a core issue to tackle the SOCproblem (11), we need to come up with a di ff erent approach. As in [48], we now use a non-deterministic reformulation of the JCC as a first passage time problem. We denote by T y = inf t ≥ t (cid:16) t | Y t > S u ( t ) (cid:17) for Y t = y t < S u ( t ) the first passage time of the OUP, and obtain the equivalent formulation ofthe JCC (9b) as a first passage time problem of the form P ( T y > t ) ≥ − θ. (14)For further details on first passage times, we refer the reader to [13, Chapter 4]. We will see thatreformulation (14) enables to include constraint (9b) into our SOC framework.However, this is not obvious. The drawback is that we have to deal with the distribution of thefirst passage time of the OUP with time-dependent mean demand level for a time-dependentabsorbing boundary. The time-dependency rules out some classical approaches. Furthermore,even for a constant boundary, the task turns out to be by far more complicated than deriving thedistribution of the first passage time of a Brownian motion (see [35, 49]). A closed-form solutionis only known for a few special cases. 7ne reason making the present case particularly hard is that there is no pure di ff usion equa-tion any more that would allow to apply the formula presented in the paper [11]. So far, to thebest of our knowledge, no closed-form solution of the first passage time density in our case isknown. Several semi-analytical approaches exist: In [28], they derive an integral representationof the first-passage time of an inhomogeneous OUP with an arbitrary continuous time-dependentbarrier and extend their results to continuous Markov processes. The idea of exploiting transfor-mations among Gauss-Markov processes to relate the problem to the known first passage timedensity of a Wiener process. However, as stated in [12], the transformation of the OUP to theBrownian motion entails exponentially large times. In an alternative approach presented in [12],they consider this broader class of real continuous Gauss-Markov process with continuous meanand covariance functions. The latter approach is the one that we tailor to our time-dependentOUP.In a first step, we show that the approach in [12] is applicable to our case of the time-dependentOUP, and in a second step, we introduce the method itself. Verification of prerequisites
As we want to apply a result on the first passage time density formulated in a general Gauss-Markov setting, we first need to verify that the time-dependent OUP given by equation (6) isindeed a Gauss-Markov process.
Lemma 3.1.
The OUP given by equation (6) is a Gauss-Markov process.Proof.
We first verify that the OUP is a Gauss process, i.e. that for any integer n ≥ ≤ t < t < · · · < t n ≤ T , the random vector (cid:16) Y t , Y t , · · · , Y t n (cid:17) has a joint normal distribution. From(8), we know that, for any n ≥
1, and arbitrary k ∈ { , · · · , n } , Y t k is normally distributed. From [4,p. 259], we know that, to verify the joint normal distribution of (cid:16) Y t , Y t , · · · , Y t n (cid:17) , it is su ffi cient toshow that any linear combination of Y t , Y t , · · · , Y t n is normally distributed. This holds true dueto the linearity of the Itˆo integral.Furthermore, the drift coe ffi cient b ( t, x ) = κ ( µ ( t ) − x ), and the di ff usion coe ffi cient σ ( t, x ) ≡ σ of(6) satisfy the assumptions for existence and uniqueness of a strong solution in [18, Thm. 3.1].Hence, [18, Thm. 3.9] is applicable, which states that the corresponding SDE, in our case (6), is aMarkov process (see [33, Def. 4.6]) on the interval [0 , T ].We are now in the Gauss-Markov setting of [12]. To prepare the numerical computation of thefirst passage time density, we need to calculate some basic characteristics of our OUP. We recall,that Y t is normally distributed with mean m OUP ( t ) = y t e − κ ( t − t ) + κ t (cid:90) t e − κ ( t − s ) µ ( s ) ds, and variance given by v OUP ( t ) = σ κ (cid:16) − e − κ ( t − t ) (cid:17) , both being C ([0 , T ])-functions. The probability density function of the OUP starting at time t in y t coincides with a normal density with mean m OUP ( t ), and variance v OUP ( t ).We proceed with the covariance function. Note that the covariance is determined by the stochas-tic integral term I t = σ (cid:82) tt e − κ ( t − s ) dW s in (7) and the deterministic part can be neglected for its8alculation. Moreover, note that E [ I t ] = 0.Cov( Y s , Y t ) = σ E (cid:34)(cid:90) st e − κ ( s − u ) dW u (cid:90) tt e − κ ( t − v ) dW v (cid:35) = σ e − κ ( s + t ) E (cid:34)(cid:90) st e κu dW u (cid:90) tt e κv dW v (cid:35) = σ e − κ ( s + t ) E (cid:34)(cid:90) tt e κu [ t ,s ] ( u ) dW u (cid:90) tt e κv dW v (cid:35) = σ e − κ ( s + t ) E (cid:34)(cid:90) st e κu du (cid:35) = σ e − κ ( t − s ) − e − κ ( s + t − t ) κ (15)The covariance function (15) can be decomposed inCov ( Y s , Y t ) = e κs (cid:16) − e − κ (2 s − t ) (cid:17)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) h ( s ) · σ κ e − κt (cid:124) (cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32) (cid:125) h ( t ) . Also, the functions h
1, and h are elements of C ([0 , T ]), and their derivatives are h (cid:48) ( t ) = κe κt + κe − κ ( t − t ) , and h (cid:48) ( t ) = − σ e − κt . Numerical calculation of the first passage time density
In [12], it is shown in a first step that the first passage time density for a C -barrier satisfies anon-singular Volterra second-kind integral equation. In a second step, this equation is iterativelysolved by a repeated Simpson’s rule yielding an approximation to the desired first passage timedensity in discretized form.Below, we state Theorem 3.1. of [12] adapted to our OUP (6) and our notation. Theorem 3.2.
Let S ( t ) be a C ([0 , T ]) -function. Then, the first passage time density g (cid:16) S ( t ) , t | y t , t (cid:17) = ∂∂ t P (cid:16) T y < t (cid:17) solves the non-singualr second-kind Volterra integral equation given by g (cid:16) S ( t ) , t | y t , t (cid:17) = − Ψ (cid:16) S ( t ) , t | y t , t (cid:17) + 2 (cid:90) tt g (cid:16) S ( t ) s, s | y t , t (cid:17) Ψ ( S ( t ) , t | S ( s ) , s ) ds, y t < S ( t ) . (16) Thereby, the function Ψ is defined via Ψ ( S ( t ) , t | y, s ) = (cid:32) S (cid:48) ( t ) − m (cid:48) OUP ( t )2 − S ( t ) − m OUP ( t )2 h (cid:48) ( t ) h ( s ) − h (cid:48) ( t ) h ( s ) h ( t ) h ( s ) − h ( t ) h ( s ) − y − m OUP ( t )2 h (cid:48) ( t ) h ( t ) − h ( t ) h (cid:48) ( t ) h ( t ) h ( s ) − h ( t ) h ( s ) (cid:33) · p y,s ( S ( t ) , t ) . We adopt the notational short cuts introduced in [12, p. 466]: g ( t ) := g (cid:16) S ( t ) , t | y t , t (cid:17) , t, t ∈ [0 , T ] , t < t Ψ ( t ) := Ψ (cid:16) S ( t ) , t | y t , t (cid:17) , t, t ∈ [0 , T ] , t < t Ψ ( t | s ) := Ψ ( S ( t ) , t | S ( s ) , s ) t, s ∈ [0 , T ] , t < s ≤ t.
9e introduce a grid t < t < ... < T N , where t k = t + k · ∆ t, k ∈ { , · · · , N } . The iterative procedurebased on the repeated Simpson’s rule to obtain an approximation ˜ g ( t k ) of the first passage timedensity g ( t k ) reads as follows:˜ g ( t ) = − Ψ ( t ) , ˜ g ( t k ) = − Ψ ( t k ) + 2 ∆ t k − (cid:88) j =1 w k,j ˜ g ( t j ) Ψ ( t k | t j ) , k = 2 , , · · · , N (17)where the weights are specified by w n, j − = 43 , j = 1 , , · · · , n ; n = 1 , , · · · , N ,w n, j = 23 , j = 1 , , · · · , n − n = 2 , , · · · , N ,w n +1 , j − = 43 , j = 1 , , · · · , n − n = 2 , , · · · , N − ,w n +1 , j = 23 , j = 1 , , · · · , n − n = 3 , · · · , N − ,w n +1 , n − = 1724 , n = 2 , , · · · , N − ,w n +1 , n − = w n +1 , n = 98 , n = 1 , , · · · , N − . The iterative procedure has been proven to converge in [12].
Theorem 3.3.
We shall be given the above discretization t < t < ... < t N , where t k = t + k · ∆ t, k ∈{ , · · · , N } , where ∆ t is the discretization step. The first passage time density obtained by the iterativeprocedure (17) converges to the true first passage time density as the step size tends to zero, i.e. lim ∆ t → | g ( t k ) − ˜ g ( t k ) | = 0 for all k ∈ { , · · · , N } To use the iteratively approximated first passage time density to obtain the risk level correspond-ing to S ( t ), we apply Algorithm 3.1. It enables to integrate the JCC (9b) in our optimization Algorithm 3.1
Algorithm to calculate the risk level corresponding to S ( t ) Input:
OUP characteristics: mean m OUP ( t ), variance v OUP ( t ), covariance decomposition in termsof h ( t ) , h ( t ), and its derivatives h (cid:48) ( t ) , h (cid:48) ( t ); boundary S ( t ); discretization step ∆ t , final activetime T CC of CC; risk level θ Output:
First passage time density g , and cumulative distribution function G Define parameters of OUP. Choose time discretization ∆ t . Calculate m OUP ( t ) , h ( t ) , h ( t ) , S ( t ) , m (cid:48) OUP ( t ) , h (cid:48) ( t ) , h (cid:48) ( t ) , S (cid:48) ( t ) for discretized time interval[0 , T CC ] with step size ∆ t . Calculate approximation of first passage time density via (17). Calculate the corresponding discrete values of the cumulative distribution function G . if G ( end ) < = θ then JCC fulfilled. else JCC violated. end if return g and G . 10able 1: Comparison of risk-MC and risk values based on Algorithm 3.1 T ∆ t M risk-MC risk-fptd di ff · (cid:42) -0.008110000 0.1435 0.1479 (cid:42) (cid:42) (cid:42) (cid:42) C -boundary S ( t ). As this control is a result of the optimization procedure, the cal-culation of the first passage time density needs to be repeated in every optimization iteration. Itis therefore worthwhile to mention that the above introduced iterative procedure is well-suitedfor computational e ffi ciency. This is because the algorithm only requires the characteristics ofthe OUP in terms of its initial data ( t , y ), its mean m OUP ( t ), its variance v OUP ( t ), its covariancedecomposition in terms of the functions h and h as well as a prespecified boundary S ( t ) and achosen discretization step ∆ t . No Monte Carlo (MC) methods, and no high-dimension integralcomputations are involved and no particular software packages are necessary (see [12]). Validation of first passage time simulation
For a validation of Algorithm 3.1, we consider the following test case setting: t = 0 , κ = / , σ =0 . , µ ( t ) = 0 . . · sin( / πt ) , y = 0 . , S ( t i ) = m OUP ( t k ) + 0 . . · t k / T , where t k , k ∈{ , · · · , N ∆ t } are the discretization points, and N ∆ t + 1 denotes the number of discretization pointsfor the interval [0 , T ] that corresponds to the chosen time step ∆ t . We denote by M the numberof MC repetitions used.In Table 1, we show the risk of hitting the boundary S ( t ) based on a MC simulation (risk-MC) andbased on the evaluation of the cumulative distribution function G from Algorithm 3.1 (risk-fptd),and calculate the di ff erences (di ff ). We observe that our algorithm 3.1 gives rather precise resultsalready for large step sizes. This is beneficial when using it within the optimization of large gasnetworks in Subsection 4.3. The MC-risk values are obtained using the plain MC method. Notethe rare event character of undersupply. This is even more pronounced in real world settings wheremost likely the chosen risk tolerance is 5% or lower instead of values around 15% in our test case.Therefore, a rare event simulation technique as e.g. in [48] in the context of power flow reliability,where the probability of an outage is very small, might lead to more accurate risk estimates evenfor larger step sizes.However, even with the plain MC method, we observe that the MC-risk and the risk-fptd valuesapproach up to a precision in the range of 10 indicating the correct functioning of the algorithm.11 .2 Reformulation of the cost function We benefit from an analytical treatment of the cost function that has been set up in [23] for thetracking type costs C t instead of only t = 0: E (cid:20)(cid:16) Y t − S u ( t ) (cid:17) | Y t = y t (cid:21) = y t e − κ ( t − t ) + 2 y e − κ ( t − t ) (cid:90) tt e − κ ( t − s ) κµ ( s ) ds + (cid:32) κ (cid:90) tt e − κ ( t − s ) µ ( s ) ds (cid:33) + σ κ (cid:16) − e − κ ( t − t ) (cid:17) − y ( t ) · (cid:32) y e − κ ( t − t ) + κ (cid:90) tt e − κ ( t − s ) µ ( s ) ds (cid:33) + y ( t ) . (18)It remains to extend the approach for the undersupply cost component C R . To this end, we make use of the truncated normal distribution. Thefollowing definition is taken from [29, p. 81] and adapted to our notation. Definition 3.1.
Let X be a random variable on a probability space ( Ω , A , P ). We say that X followsa doubly truncated normal distribution with lower and upper truncation points a and b respectively( X ∼ N ba ( ξ, σ )) if its probability density function is given by ρ X,a,b ( x ) = σ ϕ (cid:18) x − ξσ (cid:19) (cid:32) Φ (cid:32) b − ξσ (cid:33) − Φ (cid:18) a − ξσ (cid:19)(cid:33) − [ a,b ] ( x ) , where ϕ is the density, and Φ the cumulative distribution function of a standard normally dis-tributed random variable. We denote by ξ and σ the mean and the variance of the non-truncatednormal distribution.We call the distribution singly truncated from above respectively from below if a is replaced by −∞ respectively b is replaced by ∞ . Proposition 3.4 (cf. [29, p. 81]) . Let X ∼ N ba ( ξ, σ ) . Then, the expected value of X reads as E [ X ] = ξ + ϕ (cid:16) a − ξσ (cid:17) − ϕ (cid:16) b − ξσ (cid:17) Φ (cid:16) b − ξσ (cid:17) − Φ (cid:16) a − ξσ (cid:17) σ . (19)We use equation (19) denoting the expectation of a truncated normally distributed random vari-able to reformulate the expectations in C R in (10). For C
3, we obtain E (cid:104)(cid:16) S u ( t ) − Y t (cid:17) − | Y t = y t (cid:105) = (cid:90) ∞ S u ( t ) (cid:16) S u ( t ) − y (cid:17) ρ Y t ( y ) dy = S u ( t ) · (cid:16) − Φ Y t ( S u ( t ) ) (cid:17) − (cid:90) ∞ S u ( t ) yρ Y t ( y ) dy = S u ( t ) · (cid:16) − Φ Y t ( S u ( t ) ) (cid:17) − (cid:90) R yρ Y t ,S u ( t ) , ∞ ( y ) dy · P (cid:16) Y t > S u ( t ) (cid:17) = S u ( t ) · (cid:16) − Φ Y t ( S u ( t ) ) (cid:17) − (cid:32) m OUP ( t ) + ρ Y t ( S u ( t ) )1 − Φ Y t ( S u ( t ) ) (cid:112) V OUP ( t ) (cid:33) · (cid:16) − Φ Y t ( S u ( t ) ) (cid:17) , (20)where ρ Y t denotes the density, and Φ Y t the cumulative distribution function of Y t , and ρ Y t ,S u ( t ) , ∞ denotes the density of the singly truncated random variable from below by S u ( t ) . The last equationresults from formula (19) with a = S u ( t ) , and b = ∞ .12n a similar manner, by setting a = −∞ , and b = S u ( t ) in (19), for the expectation in R , we have E (cid:104)(cid:16) S u ( t ) − Y t (cid:17) + | Y t = y t (cid:105) = S u ( t ) Φ Y t ( S u ( t ) ) − (cid:32) m OUP ( t ) − ρ Y t ( S u ( t ) ) Φ Y t ( S u ( t ) ) (cid:112) V OUP ( t ) (cid:33) Φ Y t ( S u ( t ) ) . (21)With equations (18), (20), and (21), we have a completely deterministic reformulation of our costfunction (10) at hand and denote it by OF detReform ( Y t , t , y t , S u ( t ) ). This reformulation of the costfunction together with the reformulation of the CC as state constraint (13) allows us to drop theOUP (6) as constraint in the optimization problem (11).We are left with a completely deterministic PDE-constrained optimization problem, i.e. (11)without the OUP (6) as constraint and with (10) replaced by OF detReform ( Y t , t , y t , S u ( t ) ), and (9) re-placed by (13). Hence, we are free to apply a suitable method from deterministic PDE-constrainedoptimization of our choice to solve the problem. Here, we use a first-discretize-then-optimize ap-proach, and numerically calculate the discrete optimal solution based on discrete adjoints. Alltogether, we come up with the deterministic reformulation of the stochastic optimal control (SOC)problem given bymin u ∈U ad (cid:90) Tt ∗ OF detReform ( Y t ; t ; y t ; S u ( t ) ) dt subject to (1) , (13) , and/or (14) . (22)Note that the JCC (14) is still a probabilistic constraint, which can, however, be handled in adeterministic way by applying Algorithm 3.1. We apply deterministic discrete adjoint calculus (see e.g. [10, 20, 31]) to solve the determinis-tically reformulated SOC problem (22). Numerically, this is implemented in ANACONDA, amodularized simulation and optimization environment for hyperbolic balance laws on networksthat allows the inclusion of state constraints, which has been developed in [30]. Note that theinclusion of state constraint is important for the inclusion of the SCC (13). Further note thatANACONDA allows for regularization terms in the objective function to tackle finite horizone ff ects or oversteering for coarse discretizations. Both e ff ects can be reduced by penalizing con-trol variations, which has been activated within the computations in Sections 4.1 and 4.2 withpenalty parameter 10 − . To be more precise about the applied schemes, we use the IBOX scheme(2) to discretize the energy system (1), and the IPOPT solver [47] as well as the DONLP2 method[42, 43] within ANACONDA for the optimization procedure.We start this section with a validation of our numerical routine in Subsection 4.1 for a specialcase of (1), for which we are able to derive an analytical solution of the SOC problem in caseof an SCC. In Subsections 4.2 and 4.3, we apply our numerical routine in case of two di ff erentparameter settings. The choices in the energy system (1), the uncertain demand (6), and the CCs(9) stated in Table 2 correspond to the settings in Sections 4.2, and 4.3. For a numerical validation of our optimization routine, we consider the special case of the linearadvection on one edge with velocity λ ∈ R , and nonzero, linear source term sρ . As we onlyconsider the dynamics on one edge, we omit the dependence on the edge e . The latter equationresults from (1) by setting the dimension n = 1, the number of prescribed left, and right BCs to m l = 1, and m r = 0, a e = 0 and by choosing the functional relation of the left BC as Γ a = id, where13 arameter setting (PS) Tele
GtP s
GtP l T y . . µ ( t ) 1 + sin( πt ) y + 0 . πt ) 0 . . πt )Speed of mean reversion κ / Intensity of demand fluctuations σ n m l m r Γ ea ρ e p ( ρ e ) p ( ρ e )Functional relation of right BC Γ eb ρ e ρ e ρ e Flux function f ( ρ e ) (cid:32) − C − ρ e L − ρ e (cid:33) (cid:32) ρ e p ( ρ e ) + ( ρ e ) / ρ e (cid:33) (cid:32) ρ e p ( ρ e ) + ( ρ e ) / ρ e (cid:33) CC active I CC [1 . ,
4] [6 ,
12] [0 , θ
5% 5% 5%
Table 2: Specification of energy systemsid denotes the identity function, and by setting the flux function to f ( ρ ) = λρ . It reads as ∂ t ρ + λ∂ x ρ = sρρ ( x,
0) = ρ ( x ) ,ρ (0 , t ) = u ( t ) , x ∈ [0 , b ] , t ∈ [0 , T ] . (23)This IBVP has the advantage that we are able to derive an analytical solution to the correspondingSOC problem (22) with SCC (13) and can compare our numerical implementation against it. Analytical solution of the linear advection with source term
To derive the analytical solution of the SOC problem, we first need the analytical solution ofthe IBVP (23). Therefore, we divide the time interval [0 , T ] in two subintervals I = [0 , x / λ ], and I b = ( x / λ , T ] depending on the position x ∈ [0 , b ]. The solution on I is determined by the initialcondition ρ ( x ), and the solution on I b is determined by the inflow control u ( t ) acting as a left BC.For the solution of (23) on I b , we can change the role of x , and t as we deal with an empty system atthe beginning. Then, we can use the explicit solution of the IVP for the linear advection equationwith source term (see [44]). All together, the solution of (23) reads as ρ ( x, t ) = (cid:40) e st ρ ( x − λt ) t ∈ I e s x / λ u ( t − x / λ ) t ∈ I b (24)The solution of (23) is illustrated graphically in Figure 2 based on the choices λ = 4, ∆ x = 0 . ∆ t = ∆ x / λ (exact CFL-condition), T = 1, s = − ρ ( x ) = sin( x ), and u ( t ) = sin(10 · t ). Note thatthe solution at t = 0 is the prescribed initial condition. For x ∈ [0 , . t = 0 . I b , i.e. the inflow control determines the solution, and for x ∈ [0 . , t = 0 . I ,which means that, thereon, the solution results from the shifted initial condition. The solution at t = 0 .
275 depends only on the inflow control.For the particular case of (23), we are able to derive an analytical expression for the optimalcontrol for the particular cost function OF detReform with w = w = w = 0, and w = 1. Theorem 4.1.
Let the supply dynamics (1) be given by (23) , and the uncertain demand by the OUP (6) with existing second moment. Furthermore, set w = w = w = 0 , and w = 1 in (10) . Then solving (11) with (9) being an active SCC (9a) on the whole interval [0 , T ] yields the optimal control u ∗ ( t ) = m OUP ( t ) + v OUP ( t ) Φ − (1 − θ ) ∀ t ∈ [0 , T − / λ ] (25)14 Space
Figure 2: Regions of influence of IC ( I ) and inflow control ( I b ) in solution (24) Proof.
From [23], we know that the optimal supply without imposing a CC in the linear advectionsetting (23) with s = 0 for CM1 is given by the conditional expectation S u ( t ) = E [ Y t + / λ | Y ]. As thesource term, does not influence the characteristics of the linear advection (see equation (24)), weagain have a constant time delay between inflow and outflow. That enables us to still considerthe solution pointwise in time. Due to the analytical solution (24), the optimal supply is alwaysreachable as long as S u ( t ) ∈ [ u min , u max ], where u min , and u max are lower and upper bounds on thecontrol ( u min = 0, and u max = + ∞ in (11)). Since we have for the chance constraint level m OUP ( t + / λ ) + v OUP ( t + / λ ) Φ − (1 − θ ) ≥ m OUP ( t + / λ ) = E [ Y t + / λ | Y ] , and the objective function is quadratic with vertex in ( t, E [ Y t + / λ | Y ]), the optimal control is ob-tained at the left boundary of the feasible region at each point in time. Numerical validation
We now validate our numerical procedure by considering the particular SOC problem given bymin u ( t ) ,t ∈ [0 ,T ] ,u ∈ L (cid:90) Tt ∗ OF detReform ( Y t ; t ; y t ; S u ( t ) ) dt subject to (9a) , and (23) , (26)where we set w = w = w = 0, and w = 1 in (10). We consider the interval [0 , I CC = [0 . , λ = 4 and we consider a rate forthe source term of s = − .
1. In the numerical implementation, we have included a nonnegativityconstraint u ( t ) ≥ ff ect ourtest case (see Figure 3). Therefore, a direct comparison of our numerical solution against theanalytical one is possible.In our validation, we used the parameters y = 1, µ ( t ) = 1 + 2 sin(8 πt ), κ = 3 and σ = 0 . . . . . . . . . . . . . . . . . Analytical control without SCCOptimal control in x=0 without SCCAnalytical control with SCCOptimal control with SCC
Figure 3: Numerical versus analytical solution of (26)the activation time of the SCC in the numerical optimal control is inherited from the jump in theanalytical solution.
Having seen that our numerical routine provides very good approximations to the analyticalsolution in the case of the linear advection with source term, in this section, we now considerthe Telegrapher’s equation on a network, and apply our routine thereon. The correspondingenergy system is obtained by specifying the quantities in the IBVP (1) according to the valuesin Table 2 according to the parameter setting
Tele q e , q e ) T = ( U e , I e ) T , where U e is the voltage, and I e the current on edge e . We deal with a system of linear hyperbolic ba-lance laws with linear source terms, where we assume the same constant parameters R, L, C, G forresistance, inductance, capacitance, and conductance on each edge. We endow the system withinitial conditions (ICs) U e , and I e . On the bounded domain X = [ a, b ], where a = min { a e | e ∈ E} ,and b = max { b e | e ∈ E} , we prescribe boundary conditions (BCs) v ext ( t ) as an externally givenfunction to prescribe the voltage at the end of the line, and u ( t ) as the voltage control at thebeginning of the line, taking the role of the inflow control into the energy system. ∂ t U e + C − ∂ x I e = − GC − U e ∂ t I e + L − ∂ x U e = − RL − I e , (27a) U e ( x,
0) = U e ( x ) , I e ( x,
0) = I e ( x ) , (27b) U e ( a, t ) = u ( t ) , U e ( b, t ) = v ext ( t ) , (27c)Before we present our numerical results for the Telegrapher’s equation, we first address the ques-tion of well-posedness of the system on one edge, and in a second step the well-posedness of thesystem in the network case as well as the existence of an optimal control.Note that the system (27a) is diagonalizable. Hence, we can rewrite it in characteristic variablesdenoted by ξ = ( ξ + , ξ − ) T (see [22]). In case of lossless transmission ( R = G = 0), we trace the16ystem back to a decoupled system of two classical linear advection equations by splitting thedynamics into left- and right-traveling waves with characteristic speeds λ − and λ + . In the losslesscase for one edge normed to a length of 1, the IBVP (27) reads as ∂ t ξ + + λ + ∂ x ξ + = 0 ∂ t ξ − + λ − ∂ x ξ − = 0 ξ + ( x,
0) = 0 , ξ − = 0 ξ + (0 , t ) = g + ( t ) , ξ − (1 , t ) = g − ( t ) , where λ + − = + − (cid:16) √ LC (cid:17) − . Note that voltage U , and current I can be expressed in terms of thecharacteristic variables as U ( x, t ) = (cid:113) LC ( ξ + − ξ − ), and I ( x, t ) = ξ + + ξ − .As we deal with an empty system at the beginning, we can change the role of x , and t , and obtainfor λ + >
0, and λ − < / λ + ∂ t ξ + + ∂ x ξ + = 0 / λ − ∂ t ξ − + ∂ x ξ − = 0 ξ + (0 , t ) = g + ( ξ + (1 , t )) , ξ − (1 , t ) = g − ( ξ − (0 , t )) . whose solution is ξ + ( x, t ) = ξ + (0 , t − / λ + x ) = g + ( t − / λ + x ) ξ − ( x, t ) = ξ − (1 , t − / λ − x ) = g − ( t − / λ − x )By interchanging the role of x and t again (empty system at the beginning), the existence of aunique weak solution of (27) with nonzero source term is ensured as a special case of the IVPfor the system of hyperbolic balance laws with dissipative source term studied in [8]. For furtherdetails, we refer the reader to [44].As we consider the dynamics given by the Telgrapher’s equation (27a) on each edge of our net-work, we need to impose suitable coupling conditions at the inner nodes v ∈ V int = V \ ( V in ∪ V out ).Therefore, we define the set of all edges connected to node v ∈ V as E v = { e v , · · · , e vk in } ∪ { e vk in +1 , · · · , e vk in + k out } . We define the coupling function as c v : R × ( | δ + ( v ) | + | δ − ( v ) | ) → R m c , c ( q e , e ∈ δ − ( v ) ∪ δ + ( v )) = U e v ( b e v , t ) − U e v ( b e v , t ) ...U e v ( b e v , t ) − U e vk in ( b e vk in , t ) U e v ( b e v , t ) − U e vk in+1 ( a e vk in+1 , t ) U e v ( b e v , t ) − U e vk in+ k out ( a e vk in + k out , t ) (cid:80) i ∈ δ − ( v ) I ( b i , t ) − (cid:80) j ∈ δ + ( v ) I ( a j , t ) for each vertex v ∈ V int . Thereby, m c denotes the number of coupling conditions. Now, we simplystate the coupling conditions as c v = 0 for all v ∈ V int . Remark.
The coupling conditions with respect to U e v , e v ∈ E v , state the equality of the voltages atthe corresponding boundaries of all edges with common node v . The coupling condition for I e atnode v models the conservation of the flow of current.17e are interested in the existence of an optimal inflow control u ( t ) for the system (27). Similarproblems have been tackled in [9]. There, they derive the well-posedness of a system of nonlinearhyperbolic balance laws on a network with edges represented by the interval [0 , ∞ ), and the ex-istence of an optimal control minimizing a given cost functional under certain conditions in thecontext of gas networks and open canals.As our flow function f e ( ξ ) = Λ ξ is linear with λ − < < λ + , it satisfies the assumption (F) in[9]. Moreover, our source term s ( ξ ) = B is a constant function implying its Lipschitz propertyand its bounded total variation. Hence, condition (G) from [9] is also fulfilled. This gives us thewell-posedness of our IBVP (27) on a network on the positive half-line (see [9, Theorem 2.3.]).We now turn our attention to the full SOC problem (22) in consideration. An existence result ofan optimal control in a gas network setting in the presence of state constraints has been derivedin [14], where the state constraints enter the optimization problem via a barrier method. Notethat the SCC (9a) for our SOC problem can be incorporated as a state constraint. The existence ofa minimizer in the presence of a JCC (9b) is more involved as a reformulation as state constraintis no longer possible.Therefore, we focus on the numerical analysis of the SOC problem (11) with particular focus onthe influence of di ff erent types of chance constraints. Telegrapher’s equation with single chance constraint
We start with a single chance constraint (9a) in the Telegrapher’s setting on the network depictedin Figure 1 with parameters specified in Table 2 (
Tele w = w = w = 0, and w = 1. We consider the network topology from Figure 1. The initialconditions are set to U e ( x,
0) = 1 on all edges and I e ( x,
0) = 1 for e and e , I e ( x,
0) = 0 . e , e , e and e as well as I e ( x,
0) = 0 for e . We consider [ a e , b e ] = [0 ,
1] for edges e ∈ { e , . . . , e } and [ a e , b e ] = [0 ,
2] for e . We set the model parameters in (27) to R = 0 .
01 (resistance), L = 0 . C = 0 .
125 (capacitance), and G = 0 .
01 (conductance).Our numerical results in Figure 4 show that the optimally available current I (purple solid line)at node v d matches the expected value of the OUP (blue dashed line) till t = 1 .
5, and followsthe course of the optimally available current without SCC (black dotted line). For t ∈ I CC , theoptimal available current matches the given CC-Level (turquoise solid line). This goes along withour intuitive understanding: the optimal available current matches the 95%- confidence level.This is because a higher supply results in an even higher tracking-type error in the cost function,and a lower supply violates the SCC at the risk level of 5%. This numerical finding is in line withthe theoretical result in Theorem 4.1. The energy transition phase comes with a lot of challenges in its realization. Aiming for a highpercentage of energy provided from renewable energy sources, one has to somehow cope withthe large volatility in the energy generation. One possibility to react to fluctuations to still ensurea stable demand satisfaction are gas turbines. The possibility of gas-to-power, i.e. the withdrawalof gas from the gas system and its transformation to power, seems to be a promising approach. Agas turbine can be booted to full performance within several minutes. Thus, a gas turbine mightbe appropriate to overcome short-term bottlenecks in energy generation. Therefore, we finallyfocus on the modeling of gas-to-power in the context of uncertain demands. We particularly payattention how the withdrawal of gas a ff ects the behavior of the gas network.One major mathematical challenge that the gas network entails is that its governing equationsare no longer a linear system of hyperbolic balance laws. The nonlinear gas transport can bedescribed by the so-called isentropic Euler equations . They can be obtained from equation (1) by18 . . . . . . . Conditional mean OUP for CM1Available supply in x = 1 without SCCSCC levelAvailable supply in x = 1 with SCC Figure 4: Optimal supply for Telegrapher’s equation on network of Figure 1putting in the parameters of
GtP (cid:32) ρ e q e (cid:33) t + q e p ( ρ e ) + ( q e ) ρ e x = (cid:32) g ( ρ e , q e ) (cid:33) (28a) ρ e ( x,
0) = ρ e ( x ) , q e ( x,
0) = q e ( x ) . (28b)Thereby, q e is denoted by ρ e and describes the density, q e is identified with q and gives the flow, g a given source term chosen as in [17], and p e is the pressure on edge e , which is specified by thepressure law p e ( ρ e ) = d · ρ β for all edges e (with β = 1 and d = 340 ms in the examples below). Notethat the exponent β = 1 means that we deal with the isothermal Euler equations, a special case ofthe isentropic Euler equations. However, the analysis is not limited to this choice and results forvarious pressure laws can for example be found in [16].To ensure the well-posedness of the system, we also need coupling conditions at the nodes. Asin [17], we use pressure equality and mass conservation at all nodes at all times (Kirchho ff -type-coupling), i.e. for all v ∈ V and for all t ∈ [0 , T ], we have p vin ( t ) = p vout ( t ) , (cid:88) e ∈ δ − ( v ) q e ( t ) = (cid:88) ˜ e ∈ δ + ( v ) q ˜ e ( t ) . To couple the gas network to the power system, we adapt the deterministic coupled gas-to-power-system described in [17], and extend it by an uncertain power demand given by the OUP (6). Theadapted setting is sketched in Figure 5. In our case, the power system is shrinked to only one node v d , where the aggregated uncertain demand is realized. This is because our focus is on matchingthis demand Y t best possible by the power S u ( t ) provided by the gas-to-power turbine, whereasthe distribution within the power system is considered as a separate task.The power demand is attained preferably by the gas-to-power conversion amount S u ( t ) . The miss-ing power ˜ Y t needs to be covered by an external power source.The gas consumption to generate the power is described by a quadratic function (cid:15) ( S u ( t ) ) = a + a S u ( t ) + a ( S u ( t ) ) . d Y t External supply: ˜ Y t = Y t − S u ( t ) Supply through gas q = (cid:15) ( S u ( t ) )= a + a S u ( t ) + a S u ( t ) v in v v v v v v e e e e e e e compressor Figure 5: Coupled gas-to-power systemAs the amount of gas withdrawn from the network is our control, we have u ( t ) = (cid:15) ( S u ( t ) ) . Pressure bounds play an important role in gas networks. In our case, we add a lower pressurebound at node v , which is p v ( t ) ≥ bar ] (29)in the first scenario below, appearing as an additional algebraic constraint in the SOC (22).There is a pressure drop in the gas network caused by the gas withdrawal for the conversion topower. It can be compensated by a compressor station. The modeling of the compressor station istaken from [17]: the compressor is modeled as edge e with time-independent in- and outgoingpressure p v and p v , and flux values q v , and q v through the nodes v and v . We assume thatthe compressor is run via an external power source only linked to the scenario in Figure 5 viathe cost component C q v = q v . Note that the operating costs increase if the ratio p v / p v increases. Thecompressor can be controlled by the gas network operator. Therefore, a second control variable,i.e. the control of the compressor station u compr ( t ), is added to the SOC problem (22) in terms of u compr ( t ) = p v ( t ) − p v ( t ) . It appears in the operating costs C compr in the cost component C OF detReform . Moreover,costs occur for the gas consumption to satisfy the power demand. These costs complete the de-terministic costs C C C compr + 0 . · (cid:90) Tt u ( t ) dt. For our purpose, it is important to note the deterministic nature of the compressor costs C compr and the costs for the gas consumption within the objective function (10).Due to the active element in terms of the compressor station, we need to adapt two of the abovecoupling conditions at v and v as p v out ( t ) = p v in ( t ) + u compr ( t ) , and (30a) q v out ( t ) = q v in ( t ) − u ( t ) . (30b)20 . . . . . . Conditional mean OUP for CM1SCC levelAvailable supply in x = 1 with SCC Figure 6: Optimal amount of gas-to-power conversionEquation (30a) describes the gas coupling accounting for the possible pressure increase at v . Themass conservation at v except for the gas withdrawn for the gas-to-power conversion is modeledby equation (30b).It might well be the case that the power demand cannot completely be served by the gas networkdue to pressure bounds in combination with the maximal performance of the compressor station,or other physical or economical reasons. In this case, the external power supplier comes into play.Hence, the di ff erence Y t − S u ( t ) is covered by external power supply. Gas-to-power system with single chance constraint
For the numerical investigation of the impact of an SCC on the optimal amount of gas withdrawnfrom the network, we consider the parameter setting GtP s in Table 2. In the cost function (10)of the SOC (22), we consider a pure tracking type functional by setting w = w = w = 0, and w = 1.Our numerics refer to the setting schematically represented in 5. The gas network is a subgrid ofthe GasLib-40 network, which approximates a segment of the low-calorific gas network locatedin the Rhine-Main-Ruhr area in Germany [40]. An extension of the gas network by a compressorstation in a purely deterministic setting has already been analyzed in [17]. Here, we presentresults including the compressor station in the presence of an uncertain demand stream as wellas an SCC. To the best of our knowledge, this has not been considered before. We choose adiscretization of dt = 900[s] and dx ≈ a e , b e ] for the edges are due to the network specificationsin [40] and can be found in [17, Table 1] for the subgrid considered here.In Figure 6, we see that the optimal amount of gas-to-power conversion (purple solid line) followswell the course of the expected value of the OUP (6) (blue dashed line), and from t = 6 on matchesthe given SCC level (turquoise solid line) until T = 12. This again coincides with the theoreticalresult for the linear advection with source term in Theorem 4.1.In Figure 7, the necessity of the compressor station as an active element to keep the lower pres-sure bound is illustrated. There is a pressure increase at the compressor station with time-shiftto ensure the lower pressure bound at v , which is set to p v ( t ) ≡ bar ] here. This behavior hasalready been observed and investigated in a deterministic setting in [17]. As we would expect,21 . . . . P r e ss u r e i n c r e a s e [ b a r ] (a) Pressure increase at compressor station 0 2 4 6 8 10 124343 . . P r e ss u r e [ b a r ] (b) Pressure at the sink Figure 7: Pressure evolutionthis increase is particularly pronounced while the SCC is active. The lower bound on the supplyinduces the need for a higher power level. This entails a higher gas withdrawal inducing a pres-sure drop in the gas network. The latter drop needs to be compensated by the compressor stationto keep the lower pressure bound at the sink.
Gas-to-power system: comparison of single and joint chance constraint
In a next step, we consider the general objective function (10) being able to depict real costs. Todo so, we include operating costs of the compressor station (C1), costs for external power supply(C3), as well as profit of selling excess power from the gas conversion (R), and set w = w = 10 − , w = 0, and w = 10 − . Since the cost component C w = 0. To handle the terms in the cost function, we use the deterministic reformulationof the additional cost components from Subsection 3.2. Furthermore, we impose a JCC on theinterval I CC = [0 , dt = 60[s] and dx ≈ ff ers from the one in thesimulation. The amount of gas-to-power conversion can be adapted every 15 minutes. We workwith a lower pressure bound of p v ( t ) ≡ .
5. In Figure 8, we compare the optimal amount of gaswithdrawal for the JCC with the one that would be obtained for an SCC at the same risk level of θ = 5% within the same interval I CC . The grey scale shows the pointwise quantile levels of theOUP and the blue dashed line indicates the mean of the OUP. Our numerical results reveal thatthe optimal supply for the JCC (red solid line) evolves in the same structural manner as the SCC(purple dotted line) but at a higher level. This is in line with our intuitive understanding: theJCC acts pathwise and thus represents a more restrictive constraint. Note that a path that onceviolated the constraint cannot be counted as save path , i.e. a path below the imposed CC, at anylater point in time anymore. However, this is possible for an SCC.We conclude our numerical investigation of the impact of chance constraints on the optimal sup-ply with a Monte Carlo (MC) investigation of the obtained optimal supply levels for the JCC. Weintroduce two time-dependent sets of paths distinguishing those paths that at least once hit theprescribed supply boundary S ( t ) until time t from those that stay below this level until time t , i.e. Ω hit ( t ) = { ω ∈ Ω | ∃ s ∈ [0 , t ] : Y s ( ω ) > S ( t ) } , and Ω save ( t ) = { ω ∈ Ω | ∀ s ∈ [0 , t ] : Y s ( ω ) ≤ S ( t ) } .
22 0 . . . . . . . . . . Mean OUP demandOptimal supply for JCCOptimal supply for SCC . . . . . Figure 8: Optimal amount of gas to power conversion with JCC compared to SCCNote that Ω hit ∪ Ω save = Ω . We denote the elements in Ω hit hit paths , and the elements in Ω save save paths . In Figure 9a, we depict the optimal supply level with JCC as the solid blue line. Theconfidence intervals of the OUP are shown in grey scale. Note that they are calculated pointwisein time corresponding to an SCC. We investigate the running maximum of the save paths (greendotted line) and over all paths (light blue solid line), which we define via the evaluations of thenumerical approximation ˆ Y j of the OUP obtained for M = 10 MC samples. This leads to r ( j ) = max { k ∈ ,...,M } ˆ Y j ( ω k ) , j ∈ { , . . . , N ∆ t } and r save ( j ) = max { k ∈ ,...,M | ˆ Y j ( ω k ) ≤ S ( t j ) } ˆ Y j ( ω k ) , for j ∈ { , . . . , N ∆ t } , where N ∆ t + 1 is the number of time grid points. Note that the running maximum r save ( j ) isonly well-defined if the index set is nonempty. Furthermore, we depict the instances of the firstpassage time of each hit paths as black asterisks.We observe that the running maximum r save ( t ) stays close below the optimal supply S ( t ) whereasthe running maximum r hit ( t ) evolves slightly above S ( t ). This can be interpreted in the followingway: we minimize the expected quadratic deviation between our supply and the demand at themarket as a major component in our cost function. At the same time, a guarantee that no un-dersupply occurs with a 95% probability is given. Hence, it appears logical that we control ourenergy system in a way being close to this lower probabilistic undersupply bound to avoid costlyoversupply occurrences.With respect to the range where we expect the demand to evolve within (confidence intervals ingrey scale), our optimal supply with JCC appears by far too large. However, looking at the savedistance d save ( j ) = S ( t j ) − r save ( j ) in Figure 9b, we observe primarily values below 0 . ff er explaining the comparably large optimal supply.23 1 2 3 40 .
51 Time[h]
Optimal supply for JCCFirst hitting timesEvolution of r ( j ) Evolution of r save ( j ) . . . . . (a) Evolution of the running maximum . . d save (b) Evolution of the save distance d save Figure 9: MC analysis of optimal supply with JCC
In this work, we analyzed the optimal inflow control in hyperbolic energy systems subject touncertain demand and particularly focused on quantifying the related uncertainty. By imposingchance constraints, we were able to limit the risk of an undersupply to a chosen probability level θ . Thereby, we distinguished between SCCs and JCCs. Whereas it turned out that there is aquantile-based reformulation of the SCC enabling us to include the SCC as a state constraintinto the optimization, the JCC case was more involved. We came up with a first passage timereformulation of the JCC, which we solved by using an algorithm set up in the general context ofGauss-Markov processes [12].To show that our numerical procedure can be applied to real-world phenomena, we took a realgas network from the GasLib-40 library [40] and calculated the optimal amount of gas with-drawn from the network to cover an uncertain power demand stream described by an OUP in thepresence of a JCC.An interesting aspect for further research is the theoretical existence of an optimal control for theabove analyzed hyperbolic energy systems in the presence of a JCC. Acknowledgement
The authors are grateful for the support of the German Research Foundation (DFG) within theproject “
Novel models and control for networked problems: from discrete event to continuous dynam-ics ” (GO1920/4-1) and the BMBF within the project “
ENets ”. Moreover, the author Kerstin Luxwould like to thank Hanno Gottschalk for suggesting the consideration of joint chance constraintsin this context.
References [1]
D. Adelh ¨utte, D. Aßmann, T. Gonz `alez Grand `on, M. Gugat, H. Heitsch, F. Liers, R. Hen-rion, S. Nitsche, R. Schultz, M. Stingl, and D. Wintergerst , Joint model of probabilis-tic/robust (probust) constraints with application to gas network optimization . available at https://opus4.kobv.de/opus4-trr154/frontdoor/index/index/docId/215 , 2017.242]
L. Andrieu, R. Henrion, and W. R ¨omisch , A model for dynamic chance constraints in hydropower reservoir management , European J. Oper. Res., 207 (2010), pp. 579–589.[3]
M. T. Barlow , A di ff usion model for electricity prices , Math. Finance, 12 (2002), pp. 287–298.[4] H. Bauer , Probability theory , vol. 23 of De Gruyter Studies in Mathematics, Walter de Gruyter& Co., Berlin, 1996. Translated from the fourth (1991) German edition by Robert B. Burckeland revised by the author.[5]
D. Bienstock, M. Chertkov, and S. Harnett , Chance-constrained optimal power flow: Risk-aware network control under uncertainty , SIAM Rev., 56 (2014), pp. 461–495.[6]
A. Charnes and W. W. Cooper , Chance-constrained programming , Management Sci., 6 (1959),pp. 73–79.[7]
A. Charnes, W. W. Cooper, and G. H. Symonds , Cost horizons and certainty equivalents: Anapproach to stochastic programming of heating oil , Management Science, 4 (1958), pp. 235–263.[8]
C. C. Christoforou , Uniqueness and sharp estimates on solutions to hyperbolic systems withdissipative source , Comm. Partial Di ff erential Equations, 31 (2006), pp. 1825–1839.[9] R. Colombo, G. Guerra, M. Herty, and V. Schleper , Optimal control in networks of pipes andcanals , SIAM J. Control Optim., 48 (2009), pp. 2032–2050.[10]
C. D’Apice, S. G ¨ottlich, M. Herty, and B. Piccoli , Modeling, simulation, and optimizationof supply chains , Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,2010. A continuous approach.[11]
D. A. Darling and A. J. F. Siegert , The first passage problem for a continuous Markov process ,Ann. Math. Statistics, 24 (1953), pp. 624–639.[12]
E. Di Nardo, A. G. Nobile, E. Pirozzi, and L. M. Ricciardi , A computational approach to first-passage-time problems for gauss-markov processes , Advances in Applied Probability, 33 (2001),pp. 453–482.[13]
E. B. Dynkin , Markov processes. Vol. I , Die Grundlehren der Mathematischen Wissenschaften,Band 121, Academic Press Inc., Publishers, New York and Springer-Verlag, Berlin-G¨ottingen-Heidelberg, 1965. Translated with the authorization and assistance of the authorby J. Fabius, V. Greenberg, A. Maitra, G. Majone.[14]
H. Egger, T. Kugler, and W. Wollner , Numerical optimal control of instationary gas transportwith control and state constraints , preprint, (2017). last checked 01/02/2019.[15]
M. H. Farshbaf-Shaker, R. Henrion, and D. H ¨omberg , Properties of chance constraints ininfinite dimensions with an application to PDE constrained optimization , Set-Valued Var. Anal.,26 (2018), pp. 821–841.[16]
E. Fokken, S. G ¨ottlich, and O. Kolb , Modeling and simulation of gas networks coupled to powergrids , J. Engrg. Math., 119 (2019), pp. 217–239.[17]
E. Fokken, S. G ¨ottlich, and O. Kolb , Optimal control of compressor stations in a coupledgas-to-power network , in Advances in Energy System Optimization, V. Bertsch, A. Ardone,M. Suriyah, W. Fichtner, T. Leibfried, and V. Heuveline, eds., Cham, 2020, Springer Interna-tional Publishing, pp. 67–80.[18]
T. C. Gard , Introduction to stochastic di ff erential equations , vol. 114 of Monographs and text-books in pure and applied mathematics, Marcel Dekker, Inc., New York, 1988.2519] A. Geletu, M. Kl ¨oppel, H. Zhang, and P. Li , Advances and applications of chance-constrainedapproaches to systems optimisation under uncertainty , Internat. J. Systems Sci., 44 (2013),pp. 1209–1232.[20]
P. Goatin, S. G ¨ottlich, and O. Kolb , Speed limit and ramp meter control for tra ffi c flow net-works , Eng. Optim., 48 (2016), pp. 1121–1144.[21] T. Gonz ´alez Grand ´on, H. Heitsch, and R. Henrion , A joint model of probabilistic/robust con-straints for gas transport management in stationary networks , Comput. Manag. Sci., 14 (2017),pp. 443–460.[22]
S. G¨ottlich, M. Herty, and P. Schillen , Electric transmission lines: control and numericaldiscretization , Optim. Control Appl. Methods, 37 (2016), pp. 980–995.[23]
S. G ¨ottlich, R. Korn, and K. Lux , Optimal control of electricity input given an uncertain de-mand , Mathematical Methods of Operations Research, (2019).[24] ,
Optimal inflow control penalizing undersupply in transport systems with uncertain de-mands , in Progress in Industrial Mathematics at ECMI 2018, I. Farag ´o, F. Izs´ak, and P. L.Simon, eds., Cham, 2019, Springer International Publishing, pp. 485–490.[25]
C. Gotzes, H. Heitsch, R. Henrion, and R. Schultz , On the quantification of nominationfeasibility in stationary gas networks with random load , Math. Methods Oper. Res., 84 (2016),pp. 427–457.[26]
M. Gugat and M. Schuster , Stationary gas networks with compressor control and random loads:optimization with probabilistic constraints , Math. Probl. Eng., (2018), pp. Art. ID 7984079, 17.[27]
H. Heitsch, R. Henrion, H. Le¨ovey, R. Mirkov, A. M¨oller, W. R ¨omisch, and I. Wegner-Specht , Empirical observations and statistical analysis of gas demand data , in Evaluating gasnetwork capacities, T. Koch, ed., MOS-SIAM Series on Optimization, Society for Industrialand Applied Mathematics, Philadelphia, 2015, pp. 273–290.[28]
Y. Jiang, A. Macrina, and G. Peters , Multiple barrier-crossings of an ornstein-uhlenbeck di ff u-sion in consecutive periods , SSRN Journal, (2019).[29] N. L. Johnson and S. Kotz , Distributions in statistics. Continuous univariate distributions. 1 ,Houghton Mi ffl in Co., Boston, Mass, 1970.[30] O. Kolb , Simulation and Optimization of Gas and Water Supply Networks , Verlag Dr. Hut,M ¨unchen, 2011.[31]
O. Kolb and J. Lang , Simulation and continuous optimization , in Mathematical optimizationof water networks, vol. 162 of Internat. Ser. Numer. Math., Birkh¨auser/Springer Basel AG,Basel, 2012, pp. 17–33.[32]
O. Kolb, J. Lang, and P. Bales , An implicit box scheme for subsonic compressible flow withdissipative source term , Numerical Algorithms, 53 (2010), pp. 293–307.[33]
R. Korn, E. Korn, and G. Kroisandt , Monte Carlo Methods and Models in Finance and In-surance , Chapman & Hall/CRC Financial Mathematics Series, CRC Press, Boca Raton, FL,2010.[34]
D. P. Kouri and T. M. Surowiec , Existence and optimality conditions for risk-averse PDE-constrained optimization , SIAM/ASA J. Uncertain. Quantif., 6 (2018), pp. 787–815.[35]
A. Lipton and V. Kaushansky , On the first hitting time density of an ornstein-uhlenbeck process .arXiv preprint: 1810.02390, 2018. 2636]
K. Lux and S. G ¨ottlich , Optimal inflow control in transport systems with uncertain demands -a comparison of undersupply penalties , PAMM, 19 (2019), p. e201900252.[37]
T. Nesti, J. Nair, and B. Zwart , Reliability of dc power grids under uncertainty: a large devia-tions approach . arXiv preprint: 1606.02986v2, 2018.[38]
A. Pr´ekopa , Stochastic programming , vol. 324 of Mathematics and its Applications, SpringerNetherlands, 1995.[39]
R. T. Rockafellar and S. Uryasev , Optimization of conditional value-at-risk , Journal of Risk,2 (2000), pp. 21–41.[40]
M. Schmidt, D. Assmann, R. Burlacu, J. Humpola, I. Joormann, N. Kanelakis, T. Koch,D. Oucherif, M. E. Pfetsch, L. Schewe, R. Schwarz, and M. Sirvent , Gaslib - a library of gasnetwork instances , Data, 2 (2017).[41]
A. Shapiro, D. Dentcheva, and A. P. Ruszczy ´nski , Lectures on stochastic programming: Mod-eling and theory , vol. 9 of MPS-SIAM series on optimization, SIAM Soc. for Industrial andApplied Mathematics, Philadelphia, Pa., 2009.[42]
P. Spellucci , A new technique for inconsistent QP problems in the SQP method , Math. MethodsOper. Res., 47 (1998), pp. 355–400.[43] ,
An SQP method for general nonlinear programs using only equality constrained subprob-lems , Math. Programming, 82 (1998), pp. 413–448.[44]
C. Teuber , Optimal Control Methods for Transmission Lines , Dissertation, University ofMannheim, Mannheim, 2017.[45]
W. van Ackooij , Chance Constrained Programming : with applications in Energy Management ,PhD thesis, Ecole Centrale Paris, 2013.[46]
W. van Ackooij, R. Henrion, A. M ¨oller, and R. Zorgati , Joint chance constrained program-ming for hydro reservoir management , Optim. Eng., 15 (2014), pp. 509–531.[47]
A. W¨achter and L. T. Biegler , On the Implementation of a Primal-Dual Interior Point FilterLine Search Algorithm for Large-Scale Nonlinear Programming , Mathematical Programming,106 (2006), pp. 25–57.[48]
W. S. Wadman, D. T. Crommelin, and B. P. Zwart , A large-deviation-based splitting estimationof power flow reliability , ACM Trans. Model. Comput. Simul., 26 (2016), pp. 1–26.[49]
M. C. Wang and G. E. Uhlenbeck , On the theory of the Brownian motion. II , Rev. ModernPhys., 17 (1945), pp. 323–342.[50]