Chance-constrained optimal inflow control in hyperbolic supply systems with uncertain demand

In this paper, we address the task of setting up an optimal production plan taking into account an uncertain demand. The energy system is represented by a system of hyperbolic partial differential 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 preferences. The resulting output is intended to optimally match the stochastic demand for the given risk criteria. We use uncertainty quantification for an adaptation to different levels of risk aversion. More precisely, we use two types of chance constraints to formulate the requirement of demand satisfaction at a prescribed probability level. In a numerical analysis, we analyze the chance-constrained optimization problem for the Telegrapher's equation and a real-world coupled gas-to-power network.


INTRODUCTION
In recent years, significant attention has been paid to the energy market. On the one hand, this is due to climate protection policies. On the other hand, the decision of the German government on the nuclear phase-out by the end of 2022 1 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 electricity consumption generated out of renewable energy sources by the end of 2030. One major challenge is the handling of uncertainty in the renewable energy production. It heavily depends on the weather conditions and is subject to large fluctuations. Those fluctuations can heavily affect the grid operation or even lead to outages. 1 Other sources of uncertainty also play an important role in the planning of energy systems. In Reference 2, various types of uncertainties have been taken into account to solve an expansion planning problem for electricity and gas networks. One source of uncertainty in the power grid operation of particular interest in the present manuscript are power demands (see Reference 3). This makes it difficult to guarantee a sufficient supply thereby influencing reliability and profitability in power system operation (see Reference 4). It is crucial for all players in the electricity sector to cope with stochastic demand fluctuations.
The article is organized as follows: in Section 2, we introduce the stochastic optimal control setting. We define the energy system as well as the stochastic demand process, and introduce the considered cost functional. We distinguish between a single chance constraint (SCC) as well as a joint chance constraint (JCC). In Section 3, we consider a deterministic reformulation of both the cost functional and the SCC, and present a stochastic reformulation of the JCC, which allows to incorporate it in our numerical framework. In Section 4, we validate the numerical framework for a simple case of a hyperbolic supply system, that is, the scalar linear advection with source term. It represents a suitable test case setting as we are able to derive the corresponding analytical solution 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 a nonlinear system of hyperbolic balance laws in terms of a real-world gas-to-power application.
To wrap up, the novelty of our chance-constrained stochastic optimal control framework consists in a joint consideration of the very generic choice of hyperbolic balance laws to describe the supply system on the one hand and the exact deterministic reformulation of the SCC and in particular the JCC for the OUP modeling a time-dependent uncertain demand stream on the other hand. This combination of the hyperbolic nature of the time-dependent supply system (no steady-state assumption), and advantageous stochastic process modeling choice for the uncertain demand stream has not been investigated yet to the best of our knowledge.

STOCHASTIC OPTIMAL CONTROL SETTING
In this section, we set up the mathematical framework for the task of finding an optimal injection plan taking into account an uncertain demand for the time period from t 0 = 0 to final time T. This has been done for a linear transport equation in Reference 28. Here, we use the stochastic optimal control framework set up in Reference 28, 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 the stochastic demand is described by an Ornstein-Uhlenbeck process (OUP).

Energy system with uncertain demand
We consider different types of 2-dimensional energy systems. The network is modeled by a finite, connected, directed graph  = (, ) with a nonempty vertex (node) set  and a nonempty set of edges . For v ∈ , we define the set of all incoming edges by − (v) = {e ∈  ∶ e = (⋅, v)}, and the set of all outgoing edges by + (v) = {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. We further define the set of inflow vertices by  in = {v ∈  ∶ − (v) = ∅}, and the set of outflow vertices by  out = {v ∈  ∶ + (v) = ∅}. As a simplification, here, we restrict our network to | in | = | out | = 1. Note however that an extension to several inflow nodes (| in | > 1), and outflow nodes (| out | > 1) is straightforward by considering a vector-valued inflow control, and a multivariate Ornstein-Uhlenbeck process as used in Reference 27. An exemplary network structure with | in | = | out | = 1 is depicted in Figure 1. This particular network topology will later be used in Section 4.2 for the network of transmission lines and Section 4.3 for the gas network (see lower part of Figure 5). Each edge then represents a power transmission line or a gas pipeline, respectively. As can be seen in Figure 1, the inflow control u(t) acts on the vertex v in ∈  in , and the demand Y t is realized at the vertex v d ∈  out . We require b e = L for all e ∈ − (v d ). The dynamics of the energy system on one edge e are formulated in terms of a hyperbolic balance law. 29,30 This type of equation is often used when it comes to model the flow of a conserved quantity (see, eg, Reference 31). The presence of a source term in those equations accounts for a gain or a loss in this otherwise conserved quantity (eg, due to friction). The hyperbolic balance law with initial condition (IC) and boundary conditions (BCs) reads as Thereby, f ∶ R 2 → R 2 is a given flux function and s ∶ R 2 → R 2 the source term. e 0 ∶ [a e , b e ] → R 2 describes the initial state of the system on edge e. The functions Γ e a∕b ∶ R 2 → R m l∕r enable to prescribe a certain evaluation of a certain number of components of the density at the boundary. Thereby, 0 ≤ m l , m r ≤ 2 denote the number of prescribed BCs at the left respectively right boundary. A value of m l/r = 0 has to be interpreted in a way that no left/right boundary condition is prescribed. For a scalar conservation law, choosing m l = 1, m r = 0, and Γ e a = f corresponds to prescribing only the flow at the left boundary.
The boundary conditions (BCs) themselves are given in terms of the functions e a ∶ [0, T] → R m l at the left boundary a e of edge e, and e b ∶ [0, T] → R m r at the right boundary b e of edge e. The numbers of BCs m l/r , and the functions e a∕b need to be chosen carefully such that they are consistent with the characteristics of the conservation law. This will be further specified below for each setting (Sections 4.2-4.3).
Moreover, for the nondegenerate network case, that is, for each vertex need to be imposed. Thereby, m c denotes the number of coupling conditions. This will be made explicit in Subsection 4.3.
For e ∈ + (v in ), the function e a (t) depends on the inflow control u(t). To simplify notation, we do not explicitly write down the dependence of the supply on the density at the end points of all ingoing edges and denote the supply at v d at time t simply by S u(t) meaning ) . ( Clearly, the supply at v d depends on the inflow but due to the finite speed of propagation the proper formulation is via the function S u ∶ R 2×| − (v d )| → R. Thereby, the symbol ⨉ in Equation (2) generates the | − (v d )|-tuple of elements in R 2 given by the densities evaluated at the end points of corresponding ingoing edges at time t.
Formulation (1) of the energy system in terms of a hyperbolic balance law has the advantage that time-dependency can be appropriately taken into account, see, eg, Reference 32 for their usage in controlling a network of power lines. In the gas context (see, eg, Reference 33), it allows to account for nonlinear phenomena inherent to gas flow such as shocks and rarefaction waves. With model (1), we have a very generic building block for the energy system that can also be seen in the general light of the study of flows on networks. For example, the stochastic optimal control approach might also be useful for control tasks related to water flow considerations (see References 34,35).

Discretization scheme
For the numerical investigation in Section 4, the considered hyperbolic energy systems need an appropriate discretization scheme. Motivated by our real-world example, we choose an implicit box scheme (IBOX) 36 for all considered scenarios. For a general system of balance laws (on any edge) the considered scheme reads Here, Δt and Δx are the temporal and spatial mesh size, respectively, and the numerical approximation is thought in the following sense: To avoid undesired boundary effects, the discretization of the initial condition on bounded domains is done pointwise, that is, As remarked in Reference 36, for a discretization x l < x l + 1 < … < x r − 1 < 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 boundary specified by the characteristic direction. This also explains the abovementioned assumption of no change in the signature of the characteristic directions on the considered domain. The discrete version of the BCs is given by Note that the implicit box scheme has to obey an inverse CFL condition, 36 which is beneficial for 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 the choice within this work.

Uncertain demand
The uncertain demand stream is modeled by an Ornstein-Uhlenbeck process (OUP), which is a well-established choice for modeling demand uncertainty (see, eg, Reference 25). As emphasized in the Introduction, this choice of a dynamic demand model allows for a coupling of a deterministic trend and the stochastic evolution of the process. Here, we consider the OUP in terms of the unique strong solution of the following stochastic differential equation on the probability space (Ω,  , P) as it has been previously done in Reference 28. We restrict ourselves to state some of its crucial properties for our upcoming treatment here and refer the reader to Reference 28 for further assertions on the particular choice of (7). W t is a given one-dimensional Brownian motion on the same probability space, and y 0 describes the demand at time t 0 = 0. The constants > 0 and > 0 describe the speed of mean reversion and the intensity of demand fluctuations. By mean reversion, we refer to the property of the process that it is always attracted by a certain time-dependent level (t), called the mean demand level. This is due to the sign of the drift term ( (t) − Y t ), which ensures that being above (below) the mean demand level, the process experiences a reversion back to it. The latter property qualifies it to depict fluctuations around a given mean demand level while at the same time staying close (in terms of ) to it. This is a crucial property in the modeling of the demand as one usually has some deterministic prediction (eg, based on historical data) of its evolution but the stochastic nature leads to deviations therefrom. The role of the deterministic forecast in the model is taken by the time-dependent mean demand level (t). Moreover, the time-dependency of (t) allows to capture daily or seasonal patterns in the demand further motivating the choice of (7) as demand process.
From a mathematical point of view, this process has some nice analytical properties. For example, we can derive the solution to Equation (7) explicitly via the Itô formula. It reads as Moreover, it is possible to derive its distribution explicitly as For further details on the demand process, we refer the reader to Reference 28.

Chance constraints
Requiring demand satisfaction for every realization of the demand process might be too restrictive as, in some cases, it might lead to an infeasible optimization problem. 37(p5) One possibility to overcome problems of infeasibility and to reduce the average undersupply is to introduce an undersupply penalty term in the cost function. The effect of an undersupply penalty on the optimal output has been analyzed in Reference 38. A comparison of different types of undersupply can be found in Reference 39. Another approach is to guarantee with a certain probability that there is no undersupply within a 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 possibility is to require at each point in time that the probability of a demand satisfaction is at least equal to one minus a given risk level (see (10a)). This results in a so called single chance constraint (SCC). Another possibility is a joint chance constraint (JCC), that is, we require that the probability of a demand satisfaction is at least equal to one minus a given risk level on a whole interval simultaneously (see (10b)).

Objective function and stochastic optimal control problem
Having formulated all the optimization constraints, we now address the objective function. We aim 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 (C1). Another component are tracking type costs that arise from a mismatch between the externally given demand (7) and our supply S u realized at the demand vertex v d . We measure the tracking type costs in terms of the expected quadratic deviation between the demand and supply (C2). Particularly accounting for negative mismatches, we introduce undersupply costs as a third cost component (C3). If a sale of excess supply is possible, the revenue can be included into the cost function as a fourth component (R).
Putting the deterministic costs C1, the tracking costs C2, the undersupply costs C3, and the excess revenue C4 together, we obtain a cost function of the following type: where All together, we come up with the stochastic optimal control (SOC) problem )dt subject to (1), (7), and (10).
To determine the optimal control, we need to specify measurability assumptions on the control by defining the space of admissible controls  ad .
Results for this control method and two other control methods with an objective function of pure tracking type (w 1 = w 3 = w 4 = 0, w 2 = 1) for the linear advection equation without imposing CCs can be found in Reference 28.

DETERMINISTIC REFORMULATION OF THE STOCHASTIC PROBLEM
Having set up the SOC problem (12), the question of how to solve this minimization problem naturally arises. One way is to trace back the SOC problem (12) to a deterministic setting, in which we can apply well-known methods from deterministic PDE-constrained optimization such as adjoint calculus. In order to do so, in Section 3.1, we analytically treat both types of CCs presented in Subsection 2.2. The deterministic expression of the objective function introduced in Section 2.3 is derived in Section 3.2.

Reformulation of chance constraints
The reformulation of the CC heavily depends on the type of CC. Whereas the SCC (10a) can be reformulated via quantiles of a normal distribution, the reformulation of the JCC (10b) is not obvious. Therefore, we need to treat the different types of CCs separately.

Single chance constraint
For the SCCs (10a), we use a quantile-based reformulation as mentioned in Reference [13, p580]. By using the known distribution (9) of the OUP, this results in the deterministic state constraints where

Joint chance constraint
JCCs (10b) are mathematically by far more involved than SCCs (see Reference [4, p1213]). No longer considering the constraint pointwise in time, we now have to deal with a joint probability distribution. Unfortunately, we can no longer make use of the deterministic reformulation as state constraint as for (10a).
As the integration of the JCC into the optimization framework is a core issue to tackle the SOC problem (12), we need to come up with a different approach. As in Reference 27, we now use a nondeterministic reformulation of the JCC as a first passage time problem. We denote by for Y t 0 = y t 0 < S u(t 0 ) the first passage time of the OUP, and obtain the equivalent formulation of the JCC (10b) as a first passage time problem of the form For further details on first passage times, we refer the reader to Reference [40, chapter 4]. We will see that reformulation (15) enables to include constraint (10b) into our SOC framework. However, this is not obvious. The drawback is that we have to deal with the distribution of the first passage time of the OUP with time-dependent mean demand level for a time-dependent absorbing 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 the distribution of the first passage time of a Brownian motion (see References 41,42). A closed-form solution is only known for a few special cases.
One reason making the present case particularly hard is that there is no pure diffusion equation any more that would allow to apply the formula presented in the paper. 43 So far, to the best of our knowledge, no closed-form solution of the first passage time density in our case is known. Several semianalytical approaches exist: In Reference 44, they derive an integral representation of the first-passage time of an inhomogeneous OUP with an arbitrary continuous time-dependent barrier and extend their results to continuous Markov processes. The idea of exploiting transformations among Gauss-Markov processes to relate the problem to the known first passage time density of a Wiener process. However, as stated in Reference 45, the transformation of the OUP to the Brownian motion entails exponentially large times. In an alternative approach presented in Reference 45, they consider this broader class of real continuous Gauss-Markov process with continuous mean and covariance functions. The latter approach is the one that we tailor to our time-dependent OUP.
In a first step, we show that the approach in Reference 45 is applicable to our case of the time-dependent OUP, 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 (7) is indeed a Gauss-Markov process. (7 ) is a Gauss-Markov process.

Lemma 1. The OUP given by Equation
Proof. We first verify that the OUP is a Gauss process, that is, that for any integer n ≥ 1 and times 0 ≤ t 1 has a joint normal distribution. From (9), we know that, for any n ≥ 1, and arbitrary k ∈ {1, … , n}, Y t k is normally distributed. From Reference [46, p259], we know that, to verify the joint normal distribution of (7) satisfy the assumptions for existence and uniqueness of a strong solution in Reference [47, thm 3.1]. Hence, Reference [47, thm 3.9] is applicable, which states that the corresponding SDE, in our case (7), is a Markov process (see Reference [48, def 4.6 We are now in the Gauss-Markov setting of Reference 45. To prepare the numerical computation of the first passage time density, we need to calculate some basic characteristics of our OUP. We recall, that Y t is normally distributed with mean and variance given by The probability density function of the OUP starting at time t 0 in y t 0 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 stochastic integral term (8) and the deterministic part can be neglected for its calculation. Moreover, note that E[ The covariance function (16) can be decomposed in Also, the functions h1, and h 2 are elements of C 1 ([0, T]), and their derivatives are

Numerical calculation of the first passage time density
In Reference 45, it is shown in a first step that the first passage time density for a C 1 -barrier satisfies a nonsingular Volterra second-kind integral equation. In a second step, this equation is iteratively solved by a repeated Simpson's rule yielding an approximation to the desired first passage time density in discretized form. Below, we state theorem 3.1. of Reference 45 adapted to our OUP (7) and our notation.
Thereby, the function Ψ is defined via ⋅ p y,s (S(t), t).
We adopt the notational short cuts introduced in Reference [45, p466]: We introduce a grid t The iterative procedure based on the repeated Simpson's rule to obtain an approximationg(t k ) of the first passage time density g(t k ) reads as follows: where the weights are specified by The iterative procedure has been proven to converge in Reference 45.
Theorem 2. We shall be given the above discretization where Δt is the discretization step. The first passage time density obtained by the iterative procedure (18) converges to the true first passage time density as the step size tends to zero, that is, To use the iteratively approximated first passage time density to obtain the risk level corresponding to S(t), we apply Algorithm 1. It enables to integrate the JCC (10b) in our optimization procedure to solve the SOC problem (12). In our case the optimal inflow control will take the role of the C 1 -boundary S(t). As this control is a result of the optimization procedure, the calculation of the first passage time density needs to be repeated in every optimization iteration. It is therefore worthwhile to mention that the above introduced iterative procedure is well suited for computational efficiency. This is because the algorithm only requires the characteristics of the OUP in terms of its initial data (t 0 , y 0 ), its mean m OUP (t), its variance v OUP (t), its covariance decomposition in terms of the functions h 1 and h 2 as well as a prespecified boundary S(t) and a chosen discretization step Δt. No Monte Carlo (MC) methods, and no high-dimension integral computations are involved and no particular software packages are necessary (see Reference 45).

Algorithm 1. Algorithm to calculate the risk level corresponding to S(t)
Require: OUP characteristics: mean m OUP (t), variance v OUP (t), covariance decomposition in terms of h 1 (t), h 2 (t), and its derivatives h ′ 1 (t), h ′ 2 (t); boundary S(t); discretization step Δt, final active time T CC of CC; risk level Ensure: First passage time density g, and cumulative distribution function G 1: Define parameters of OUP. 2: Choose time discretization Δt.  In Table 1, we show the risk of hitting the boundary S(t) based on a MC simulation (risk-MC) and based on the evaluation of the cumulative distribution function G from Algorithm 1 (risk-fptd), and calculate the differences (diff). We observe that our Algorithm 1 gives rather precise results already for large step sizes. This is beneficial when using it within Note the rare event character of undersupply. This is even more pronounced in real-world settings where most 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, for example, in Reference 27 in the context of power flow reliability, where the probability of an outage is very small, might lead to more accurate risk estimates even for larger step sizes. However, even with the plain MC method, we observe that the MC-risk and the risk-fptd values approach up to a precision in the range of 10 3 indicating the correct functioning of the algorithm.

Reformulation of the cost function
We benefit from an analytical treatment of the cost function that has been set up in Reference 28 for the tracking type costs C2 of (11). We use the slightly more general formulation of it allowing for arbitrary t 0 instead of only t 0 = 0: It remains to extend the approach for the undersupply cost component C3 and the excess supply revenue component R.
To this end, we make use of the truncated normal distribution. The following definition is taken from Reference [49, p.81] and adapted to our notation.

Definition 1.
Let X be a random variable on a probability space (Ω, , P). We say that X follows a doubly truncated normal distribution with lower and upper truncation points a and b, respectively (X ∼  b a ( , 2 )) if its probability density function is given by where is the density, and Φ the cumulative distribution function of a standard normally distributed random variable. We denote by and 2 the mean and the variance of the nontruncated normal distribution. We call the distribution singly truncated from above respectively from below if a is replaced by −∞ respectively b is replaced by ∞.
Proposition 1 (cf Reference [49, p81]). Let X ∼  b a ( , 2 ). Then, the expected value of X reads as We use Equation (20) denoting the expectation of a truncated normally distributed random variable to reformulate the expectations in C3 and R in (11). For C3, we obtain where Y t denotes the density, and Φ Y t denotes 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 equation results from formula (20) with a = S u(t) , and b = ∞.
In a similar manner, by setting a = −∞, and b = S u(t) in (20), for the expectation in R, we have With Equations (19), (21), and (22), we have a completely deterministic reformulation of our cost function (11) at hand and denote it by OF detReform (Y t , t 0 , y t 0 , S u(t) ). This reformulation of the cost function together with the reformulation of the CC as state constraint (14) allows us to drop the OUP (7) as constraint in the optimization problem (12).
We are left with a completely deterministic PDE-constrained optimization problem, that is, (12) without the OUP (7) as constraint and with (11) replaced by OF detReform (Y t , t 0 , y t 0 , S u(t) ), and (10) replaced by (14). Hence, we are free to apply a suitable method from deterministic PDE-constrained optimization of our choice to solve the problem. Here, we use a first-discretize-then-optimize approach, and numerically calculate the discrete optimal solution based on discrete adjoints. All together, we come up with the deterministic reformulation of the stochastic optimal control (SOC) problem given by (14), and/or (15).
Note that the JCC (15) is still a probabilistic constraint, which can, however, be handled in a deterministic way by applying Algorithm 1.

NUMERICAL RESULTS
We apply deterministic discrete adjoint calculus (see, eg, References 50-52) to solve the deterministically reformulated SOC problem (23). Numerically, this is implemented in ANACONDA, a modularized simulation and optimization environment for hyperbolic balance laws on networks that allows the inclusion of state constraints, which has been developed in Reference  To be more precise about the applied schemes, we use the IBOX scheme (3) to discretize the energy system (1), and the IPOPT solver 53 as well as the DONLP2 method 54,55 within ANACONDA for the optimization procedure. We start this section with a validation of our numerical routine in Subsection 4.1 for a special case of (1), for which we are able to derive an analytical solution of the SOC problem in case of an SCC. In Subsections 4.2 and 4.3, we apply our numerical routine in case of three different parameter settings. The choices in the energy system (1), the uncertain demand (7), and the CCs (10) stated in Table 2 correspond to the settings for the Telegrapher's equation (Tele) in Section 4.2 and the ones for the gas-to-power system (GtP) in Section 4.3.

Validation via scalar linear advection with source term
For a numerical validation of our optimization routine, we consider the special case of the linear advection on one edge with velocity ∈ R, and nonzero, linear source term s . As we only consider the dynamics on one edge, we omit the dependence on the edge e. The latter equation results 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, where id denotes the identity function, and by setting the flux function to f ( ) = . It reads as This IBVP has the advantage that we are able to derive an analytical solution to the corresponding SOC problem (23) with SCC (14) 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 of the IBVP (24). Therefore The solution on I 0 is determined by the initial condition 0 (x), and the solution on I b is determined by the inflow control u(t) acting as a left BC. For the solution of (24) on I b , we can change the role of x, and t as we deal with an empty system at the beginning. Then, we can use the explicit solution of the IVP for the linear advection equation with source term (see Reference 56). All together, the solution of (24) reads as The solution of (24) is illustrated graphically in Figure 2  For the particular case of (24), we are able to derive an analytical expression for the optimal control for the particular cost function OF detReform with w 1 = w 3 = w 4 = 0, and w 2 = 1. (24 ), and the uncertain demand by the OUP (7 ) with existing second moment. Furthermore, set w 1 = w 3 = w 4 = 0, and w 2 = 1 in (11). Then solving (12) with (10) being an active SCC (10a) on the whole interval [0, T] yields the optimal control

Theorem 3. Let the supply dynamics (1 ) be given by
Proof. From Reference 28, we know that the optimal supply without imposing a CC in the linear advection setting (24) with s = 0 for CM1 is given by the conditional expectation S u(t) = E[Y t+ 1 |Y 0 ]. As the source term, does not influence the characteristics of the linear advection (see Equation (25)), we again have a constant time delay between inflow and outflow. That enables us to still consider the solution pointwise in time. Due to the analytical solution (25), the optimal supply is always reachable as long as S u(t) ∈ [u min , u max ], where u min , and u max are lower and upper bounds on the control (u min = 0, and u max = +∞ in (12)). Since we have for the chance constraint level

Numerical validation
We now validate our numerical procedure by considering the particular SOC problem given by where we set w 1 = w 3 = w 4 = 0, and w 2 = 1 in (11). We consider the interval [0, 1]. For this setting, we derived the analytical optimal control in Theorem 3 such that we can compare our numerical solution against the analytical one. We activate the SCC within I CC = [0. 6,1]. The transport velocity is = 4 and we consider a rate for the source term of s = −0.1. In the numerical implementation, we have included a nonnegativity constraint u(t) ≥ 0 on the inflow control. Note however that this constraint does not affect our test case (see Figure 3). Therefore, a direct comparison of our numerical solution against the analytical one is possible.
In our validation, we used the parameters y 0 = 1, (t) = 1 + 2 sin(8 t), = 3 and = 0.1 for the OUP. In Figure 3, we observe that the numerical optimal control is close to the analytical control (26). This finding suggests the correct functioning of our routine. Note that the small spike at the activation time of the SCC in the numerical optimal control is inherited from the jump in the analytical solution.

Linear system of hyperbolic balance laws: Telegrapher's equation
Having seen that our numerical routine provides very good approximations to the analytical solution in the case of the linear advection with source term, in this section, we now consider the Telegrapher's equation on a network to model power flow along transmission lines as it has already been done in References 32,57,58. This model is generic in the sense that a steady-state consideration of the Telegrapher's equation leads to the well-known power flow equations. We apply our routine to the Telegrapher's equation on a network here. The corresponding energy system is obtained by specifying the quantities in the IBVP (1) according to the values in Table 2 according to the parameter setting Tele 4.2. In the following assertions, we stick to the conventional notation of the Telegrapher's equation and denote (q1e, q2e) T = (U e , I e ) T , where U e is the voltage, and I e is the current on edge e. We deal with a system of linear hyperbolic balance laws with linear source terms, where we assume the same constant parameters R, L, C, G for resistance, inductance, capacitance, and conductance on each edge. We Before we present our numerical results for the Telegrapher's equation, we first address the question of well-posedness of the system on one edge, and in a second step the well-posedness of the system in the network case as well as the existence of an optimal control. Note that the system (28a) is diagonalizable. Hence, we can rewrite it in characteristic variables denoted by = ( + , − ) T (see Reference 32). In case of lossless transmission (R = G = 0), we trace the system back to a decoupled system of two classical linear advection equations by splitting the dynamics into left-and right-traveling waves with characteristic speeds − and + . In the lossless case for one edge normed to a length of 1, the IBVP (28) reads as Note that voltage U, and current I can be expressed in terms of the characteristic variables as , and I(x, t) = + + − . As we deal with an empty system at the beginning, we can change the role of x, and t, and obtain for + > 0, and − < 0 0, t)), whose solution is By interchanging the role of x and t again (empty system at the beginning), the existence of a unique weak solution of (28) with nonzero source term is ensured as a special case of the IVP for the system of hyperbolic balance laws with dissipative source term studied in Reference 59. For further details, we refer the reader to Reference 56.
As we consider the dynamics given by the Telgrapher's equation (28a) on each edge of our network, we need to impose suitable coupling conditions at the inner nodes v ∈  int =  ⧵ ( in ∪  out ). Therefore, we define the set of all edges connected to node v ∈  as We define the coupling function as for each vertex v ∈  int . Thereby, m c denotes the number of coupling conditions. Now, we simply state the coupling conditions as c v = 0 for all v ∈  int . Remark 1. The coupling conditions with respect to U e v , e v ∈  v , state the equality of the voltages at the corresponding boundaries of all edges with common node v. The coupling condition for I e at node v models the conservation of the flow of current.
We are interested in the existence of an optimal inflow control u(t) for the system (28). Similar problems have been tackled in Reference 60. There, they derive the well-posedness of a system of nonlinear hyperbolic balance laws on a network with edges represented by the interval [0, ∞), and the existence of an optimal control minimizing a given cost functional under certain conditions in the context of gas networks and open canals.
As our flow function f e ( ) = Λ is linear with − < 0 < + , it satisfies the assumption (F) in Reference 60. Moreover, our source term s( ) = B is a constant function implying its Lipschitz property and its bounded total variation. Hence, condition (G)from Reference 60 is also fulfilled. This gives us the well-posedness of our IBVP (28) on a network on the positive half-line (see Reference 60, theorem 2.3).
We now turn our attention to the full SOC problem (23) in consideration. An existence result of an optimal control in a gas network setting in the presence of state constraints has been derived in Reference 61, where the state constraints enter the optimization problem via a barrier method. Note that the SCC (10a) for our SOC problem can be incorporated as a state constraint. The existence of a minimizer in the presence of a JCC (10b) is more involved as a reformulation as state constraint is no longer possible.
Therefore, we focus on the numerical analysis of the SOC problem (12) with particular focus on the influence of different types of chance constraints.

Telegrapher's equation with single chance constraint
We start with a single chance constraint (10a) in the Telegrapher's setting on the network depicted in Figure 1 with parameters specified in Table 2 (Tele 4.2), and consider the cost function (11) with w 1 = w 3 = w 4 = 0, and w 2 = 1. We consider the network topology from 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 follows the course of the optimally available current without SCC (black dotted line). For t ∈ I CC , the optimal available current matches the given CC-Level (turquoise solid line). This goes along with our 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 with the theoretical result in Theorem 3.

Nonlinear system of hyperbolic balance laws: Gas-to-power system
The energy transition phase comes with a lot of challenges in its realization. Aiming for a high percentage of energy provided from renewable energy sources, one has to somehow cope with the large volatility in the energy generation. One possibility to react to fluctuations to still ensure a stable demand satisfaction are gas turbines. The possibility of gas-to-power, that is, the withdrawal of gas from the gas system and its transformation to power, seems to be a promising approach. A gas turbine can be booted to full performance within several minutes. Thus, a gas turbine might be appropriate to overcome short-term bottlenecks in energy generation. Therefore, we finally focus on the modeling of gas-to-power in the context of uncertain demands. We particularly pay attention how the withdrawal of gas affects the behavior of the gas network. Hyperbolic balance laws can capture complex phenomena inherent to gas flow such as shocks or rarefaction waves. We refer the reader to Reference 29 for more details. To account for the latter, we do not make a steady-state assumption here (see also Reference 62), which clearly complicates the control task.
Using a time-dependent model for the gas flow, one major mathematical challenge is that its governing equations are no longer a linear system of hyperbolic balance laws. The nonlinear gas transport can be described by the so-called isentropic Euler equations. They can be obtained from Equation (1) by putting in the parameters of GtP_a 4.3 of Table 2: Thereby, q1e is denoted by e and describes the density, q 2 e is identified with q and gives the flow, g a given source term chosen as in Reference 63, and p e is the pressure on edge e, which is specified by the pressure law p e ( e ) = d 2 ⋅ for all edges e (with = 1 and d = 340 m s in the examples below). Note that the exponent = 1 means that we deal with the isothermal Euler equations, a special case of the isentropic Euler equations. However, the analysis is not limited to this choice and results for various pressure laws can for example be found in Reference 62.
To ensure the well-posedness of the system, we also need coupling conditions at the nodes. As in Reference 63, we use pressure equality and mass conservation at all nodes at all times (Kirchhoff-type-coupling), that is, for all v ∈  and for all t ∈ [0, T], we have

q̃e(t).
To couple the gas network to the power system, we adapt the deterministic coupled gas-to-power-system described in Reference 63, and extend it by an uncertain power demand given by the OUP (7). The adapted 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 matching this demand Y t best possible by the power S u(t) provided by the gas-to-power turbine, whereas the 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 missing powerỸ t needs to be covered by an external power source.
The gas consumption to generate the power is described by a quadratic function (S u(t) ) = a 0 + a 1 S u(t) + a 2 (S u(t) ) 2 .
As the amount of gas withdrawn from the network is our control, we have u(t) = (S u(t) ).
Pressure bounds play an important role in gas networks. In our case, we add a lower pressure bound at node v 5 , which is in the first scenario below, appearing as an additional algebraic constraint in the SOC (23). There is a pressure drop in the gas network caused by the gas withdrawal for the conversion to power. It can be compensated by a compressor station. The modeling of the compressor station is taken from Reference 63: the compressor is modeled as edge e 1 with time-independent in-and outgoing pressure p v 0 and p v 1 , and flux values q v 0 , and q v 1 through the nodes v 0 and v 1 . We assume that the compressor is run via an external power source only linked to the scenario in Figure 5 via the cost component C1 in (11) accounting for operator costs of the compressor. This entails the flux equality q v 1 = q v 0 . Note that the operating costs increase if the ratio p v 0 increases. The compressor can be controlled by the gas network operator. Therefore, a second control variable, that is, the control of the compressor station u compr (t), is added to the SOC problem (23) in terms of It appears in the operating costs C1 compr in the cost component C1 of OF detReform . Moreover, costs occur for the gas consumption to satisfy the power demand. These costs complete the deterministic costs C1 in (11) as For our purpose, it is important to note the deterministic nature of the compressor costs C1 compr and the costs for the gas consumption within the objective function (11).
Due to the active element in terms of the compressor station, we need to adapt two of the above coupling conditions at v 1 and v 2 as Equation (31aa) describes the gas coupling accounting for the possible pressure increase at v 1 . The mass conservation at v 2 except for the gas withdrawn for the gas-to-power conversion is modeled by Equation (31ab).
It might well be the case that the power demand cannot completely be served by the gas network due 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 difference 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 withdrawn from the network, we consider the parameter setting GtP_a in Table 2. In the cost function (11) of the SOC (23), we consider a pure tracking type functional by setting w 1 = w 3 = w 4 = 0, and w 2 = 1.
Our numerics refer to the setting schematically represented in Figure 5. The gas network is a subgrid of the GasLib-40 network, which approximates a segment of the low-calorific gas network located in the Rhine-Main-Ruhr area in Germany. 64 An extension of the gas network by a compressor station in a purely deterministic setting has already been analyzed in Reference 63. Here, we present results including the compressor station in the presence of an uncertain demand stream as well as an SCC. To the best of our knowledge, this has not been considered before. We choose a discretization of dt = 900 seconds and dx ≈ 1 km, slightly adapted to the individual length of the edges. The specifications of the intervals [a e , b e ] for the edges are due to the network specifications in Reference 64 and can be found in Reference [63, table 1] for the subgrid considered here.
In Figure 6, we see that the optimal amount of gas-to-power conversion (purple solid line) follows well the course of the expected value of the OUP (7) (blue dashed line), and from t = 6 on matches the given SCC level (turquoise solid line) until T = 12. This again coincides with the theoretical result for the linear advection with source term in Theorem 3.
In Figure 7, the necessity of the compressor station as an active element to keep the lower pressure bound is illustrated. There is a pressure increase at the compressor station with time-shift to ensure the lower pressure bound at v 5 , which is set to p v 5 (t) ≡ 43 (bar) here. This behavior has already been observed and investigated in a deterministic setting in Reference 63. As we would expect, this increase is particularly pronounced while the SCC is active. The lower bound on the supply induces the need for a higher power level. This entails a higher gas withdrawal inducing a pressure drop in the gas network. The latter drop needs to be compensated by the compressor station to 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 (11) being able to depict real costs. To do 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 1 = w 3 = 10 −4 , w 2 = 0, and w 4 = 10 −6 . Since the cost component C2 cannot directly be interpreted as real costs, we set w 2 = 0. To handle the terms in the cost function, we use the deterministic reformulation of the additional cost components from Subsection 3.2. Furthermore, we impose a JCC on the interval I CC = [0, 4]. We consider a discretization of dt = 60 seconds and dx ≈ 1 km, again slightly adapted to the individual length of the edges, applied to the same subgrid of the GasLib-40 network 64 as above. Note that this coarse discretization is possible due to the chosen numerical scheme, the IBOX scheme (3). The coarse time discretization is possible due to the properties of the IBOX scheme introduced in Subsection 2.1. The control grid differs from the one in the simulation. The amount of gas-to-power conversion can be adapted every 15 minutes. We work with a lower pressure bound of p v 5 (t) ≡ 42.5. The remaining parameters are stated as GtP_b in Table 2. In Figure 8, we compare the optimal amount of gas withdrawal 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 gray scale shows the pointwise quantile levels of the OUP and the blue dashed line indicates the mean of the OUP.
Our numerical results reveal that the 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: the JCC acts pathwise and thus represents a more restrictive constraint. Note that a path that once violated the constraint cannot be counted as save path, that is, a path below the imposed CC, at any later 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 supply with a Monte Carlo (MC) investigation of the obtained optimal supply levels for the JCC. We introduce two time-dependent sets of paths distinguishing those paths that at least once hit the prescribed supply boundary S(t) until time t from those that stay below this level until time t, that is, Note 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. The confidence intervals of the OUP are shown in gray scale. Note that they are calculated pointwise in time corresponding to an SCC. We investigate the running maximum of the save paths (green dotted line) and over all paths (light blue solid line), which we define via the evaluations of the numerical approximationŶ j of the OUP obtained for M = 10 3 MC samples. This leads to r(j) = max where N Δt + 1 is the number of time grid points. Note that the running maximum r save (j) is only well defined if the index set is nonempty. Furthermore, we depict the instances of the first passage 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), whereas the running maximum r hit (t) evolves slightly above S(t). This can be interpreted in the following way: we minimize the expected quadratic deviation between our supply and the demand at the market as a major component in our cost function. At the same time, a guarantee that no undersupply occurs with a 95% probability is given. Hence, it appears logical that we control our energy system in a way being close to this lower probabilistic undersupply bound to avoid costly oversupply occurrences.
With respect to the range where we expect the demand to evolve within (confidence intervals in gray scale), our optimal supply with JCC appears by far too large. However, looking at the save distance d save (j) = S(t j ) − r save (j) in Figure 9B, we observe primarily values below 0.1 indicating a very limited buffer explaining the comparably large optimal supply.

CONCLUSION
In this work, we analyzed the optimal inflow control in hyperbolic energy systems subject to uncertain demand and particularly focused on quantifying the related uncertainty. By imposing chance 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 a quantile-based reformulation of the SCC enabling us to include the SCC as a state constraint into the optimization, the JCC case was more involved. We came up with a first passage time reformulation of the JCC, which we solved by using an algorithm set up in the general context of Gauss-Markov processes. 45 To show that our numerical procedure can be applied to real-world phenomena, we took a real-gas network from the GasLib-40 library 64 and calculated the optimal amount of gas withdrawn from the network to cover an uncertain power demand stream described by an OUP in the presence of a JCC.
An interesting aspect for further research is the theoretical existence of an optimal control for the above analyzed hyperbolic energy systems in the presence of a JCC.