Prostate cancer is a worldwide problem. While the role of caveolin-1 has been extensively studied, little is known about the role of caveolin-2 (CAV2) in prostate cancer. Up-regulation of CAV2 in androgen independent PC3 cells compared to normal prostate cell line and androgen dependent prostate cancer cell lines has been observed. Recent studies suggest that up-regulation of CAV2 plays an important role in androgen independent prostate cancer. This study investigates whether CAV2 is important in mediating the aggressive phenotypes seen in androgen independent prostate cancer cells. The androgen independent prostate cancer cell line, PC3 was used that has been shown to express CAV2, and CAV2 knock down was performed using siRNA system. Changes to cell number, migration and invasion were assessed after knocking down CAV2. Our results showed that down-regulating CAV2 resulted in reduced cell numbers, migration and invasion in PC3 cells. This preliminary study suggests that CAV2 may act to promote malignant behavior in an androgen independent prostate cancer cell line. Further studies are required to fully elucidate the role of CAV2 in androgen independent prostate cancer.
1.
Introduction
Modeling, simulation and optimization of critical infrastructure systems such as traffic, electricity, water or natural gas networks play an increasingly important role in our society. Many of these problems involve aspects of dynamic energy transportation and distribution in networks. The optimization with respect to efficiency, robustness, or environmental performance requires the use of high-resolution dynamical models in form of time-dependent differential equations. As a particular case, we consider transportation and distribution of natural gas in a network of pipelines, where the high-resolution models refer to the one-dimensional Euler gas equations coupled with further dynamics representing active or passive network elements such as compressors, resistors and control valves.
Detailed models for gas flow in pipeline networks are well established, see e.g., [4,9,17]. Similar maturity has been achieved for time-simulation methods, see e.g., [1,16,29], the analysis of stationary states [22], the time-continuous optimal control of compressors using adjoint techniques [26,37], and feedback stabilization based on classical solutions [15]. If the discrete nature of control valves being open or closed is to be taken into account, then the optimization needs to deal with mixed integer and continuous type variables simultaneously. Similar decisions are often to be taken into account also in other critical infrastructure systems.
In industrial practice the treatment of switched systems is often carried out by first applying a full space-time discretization of the systems and then using a mixed-integer nonlinear programming tool that incorporates the switches as extra variables to be optimized, see e.g., [7,5,47]. Another approach is to restrict the optimization to the stationary case where special purpose techniques can be successfully applied [44]. A temporal expansion of these techniques on a full network currently seems to be out of reach, see, e.g., [23,36]. We proceed in yet another way via a model switching approach, see [24]. Since networks already provide a natural spatial partition by the edges representing pipes, it seems natural to minimize a global model error by selecting either one of several models in a dynamic model hierarchy [16,18] or a stationary model hierarchy [37] for each pipe as a function of time. To do this, it is important to identify regions in the time expanded network problem where stationary models still provide a reasonable approximation in the sense that the global error remains small. A particular difficulty that the model switching problem for gas networks has in common with most of the other mentioned applications is the high-resolution model being a partial differential equation (typically a system of balance laws). While the computation of optimal switching for ordinary differential equations and differential-algebraic equations is theoretically and numerically well studied, see [11,19,27,28,33,38,51,52,53,54], the corresponding theory involving certain types of partial differential equations is still under development [45,46].
Our main contribution in this paper is to provide a theoretical and numerical framework for solving the model switching problem using the example of gas networks. We show that the problem can be casted in the sense of switching among a family of abstract evolution equations on an appropriate Banach space. This allows us to use adjoint based gradient representations for switching time and mode sequence variations recently developed in [45] to characterize locally optimal solutions for the model switching problem by introducing the notion of first order stationarity. This, in turn, motivates a two stage gradient descent approach conceptually introduced and analyzed in [2] for the optimization of switching sequences in the context of ordinary differential equations for numerical solutions of the model switching problem. We provide results for a proof-of-concept implementation for a gas network comprising 340 km of pipes on a 30 min time horizon.
Our approach relies on a semigroup property for networked transport systems. For related results without nodal control, see [6,30,41], and boundary conditions of a delay-differential-type are considered in [48,49]. In [20,21], classical solutions for a system separating a semigroup equation from a linear nodal control condition are derived. The recent work [31] considers mild solutions if the nodal control is semilinear. In contrast, we derive a semigroup formulation for the entire system, allowing for semilinear dynamics both on the edges and the boundary control, based on a characteristic decomposition of the system.
The paper is organized a follows. In Section 2 we provide a more detailed description of the common gas network models. In particular, we briefly discuss a semilinear simplification of the Euler gas equations as the most detailed model on a pipe and consider the corresponding stationary solutions. Moreover, we introduce the model switching problem. In Section 3, we show that the equations on a network coupling these models for the pipes, along with appropriate coupling conditions on the nodes allowing further network elements such as valves and compressors, is well-posed in the sense of being equivalently represented by a nonlinear perturbation of a strongly continuous semigroup in a suitable Sobolev space. In Section 4 we apply the well-posedness result to define an appropriate system of adjoint equations and to derive a stationarity concept based on switching-time gradients and mode-insertion gradients as a first order optimality condition for model optimality. Further, we present details of a conceptual algorithm alternating between switching-time optimization and mode-insertion in order to compute a solution of the model switching problem in the sense of our stationarity concept. In Section 5 we present a numerical study. In Section 6 we discuss applications and directions of future work.
2.
Gas network modeling
2.1. A model hierarchy for single pipes
The motion of a compressible non-viscous gas in long high-pressure pipelines is described by the one-dimensional isothermal Euler equations. They are given by a system of nonlinear hyperbolic partial differential equations (PDEs) and consist of the continuity equation and the balance of moments (see, e.g., [9,34,35,50])
where ϱ denotes the density in kgm3, v the velocity in ms, P the pressure of the gas in kgms2, g the gravitational constant and h′ the slope of the pipe. Furthermore θ=λ2D, where λ is the friction coefficient of the pipe, and D is the diameter of the pipe. The conserved, respectively balanced, quantities of the system are the density ϱ and the flux q=ϱv. Pressure and density are related by the constitutive law for a real gas
where Z(P,T0) is the gas compressibility factor at constant temperature T0 and Rs is the specific gas constant. For an ideal gas one has Z(P)≡1. This system of equations yields a hierarchy of simplified models for pipe dynamics [9,17,25,42].
For instationary situations, we consider the semilinear PDE model
with a constant speed of sound c=√P/ϱ obtained under the assumption of a constant gas compressibility factor Z(P,T0)≡ˉZ. This model neglects inertia effects and assumes small velocities |v|≪c. For natural gas, indeed one typically has |v|≤10ms and c≈340ms, see e.g., [17]. This PDE exhibits two simple characteristics with speeds λ1=−c and λ2=c.
For stationary situations, we consider the model
obtained from (2) with ∂tϱ=∂tq=0. Here, the flux q is constant in space and time q(t,x)=ˉq and ϱ(t,x)=ˉϱ(x) with ˉϱ being a solution of the momentum equation in (3), which is a Bernoulli-equation. A solution of the momentum equation can therefore be obtained algebraically
Analogous considerations apply also for the case of non-isothermal flow, see e.g. [9,17]. More generally, similar model hierarchies also exist for other infrastructure systems such as water distribution networks [24].
2.2. Networks with pipes, valves, and compressors
For m,n∈N we consider a network of pipes that we model by a metric graph G=(V,E) with nodes V=(v1,…,vm) and edges E=(e1,…,en)⊆V×V. For each edge e∈E, call e(1) the left node and e(2) the right node of e. We demand the incident nodes of every edge to be different, so e(1)≠e(2) for any e∈E and thus self-loops are not allowed. On the other hand, if v∈V is any node, then we define
The number |δv| then is called the degree of node v∈V.
With each edge ej∈E of such a network, we associate a pipe model from the hierarchy in Section 2.1 and a given pipe length Lj>0. Furthermore, depending on the role of each node in the network, we impose appropriate coupling conditions for the gas density and flow at the boundary of pipes corresponding to edges being incident to that, see [13]. To this end, we define for v∈V and ej∈δv
For each node v∈V, we then impose a transmission condition for the density and a balance equation for the fluxes at the node. The transmission condition states that the density variables ϱj weighted by given factors α∈(0,∞)m×2 coincide for all incident edges e∈δv and can be expressed as
The nodal balance equation for a given outflow function qv:[0,T]→R is similar to a classical Kirchhoff condition for the fluxes qj and can be written as
The choice of α corresponds to the nodal types in the network. Prototypically, we will consider the following cases:
Junctions: nodes v such that qv≡0 and αkx(v,ek)=1 for all ek∈δv.
Boundary nodes: nodes v such that αkx(v,ek)=1 for all ek∈δv, but qv≢0.
We also refer to v as an entry node, if qv<0, or an exit node, if qv>0.
Compressors: nodes v with qv≡0 and |δ+v|=|δ−v|=1. A description established via the characteristic diagram based on measured specific changes in adiabatic enthalpy Had of the compression process yields the model
where κ is a compressor specific constant, ˉZ is the gas compressibility factor that is assumed to be constant and Had is within flow dependent and compressor specific bounds obtained from the characteristic diagram. In consistency with the pipe models, we assume that Had is given by a known reference ˉHad. Then we get
with a compressor specific factor
This yields αk1=1 and αl0=ˉα. The compression ratio of centrifugal and turbo compressors mostly lies within the range ˉα∈[1,2], with ˉα≈1.3 being a typical value, see [13], [39,Chapter 4] and [12,Chapter 5.1]. For realistic long distance gas networks, compressors are set along the pipes in a distance of around 100-300 km, see [12,Chapter 5.3].
We note that there are other conceivable ways to model gas pipe junctions, including geometry conditions, leading to possibly nonlinear coupling conditions, see e.g., [3,9,40]. Furthermore, gas networks typically involve additional network components such as valves, resistors, gas coolers, etc. For appropriate coupling conditions we refer to [44].
2.3. The model switching problem
We now discuss switching between the two models (2) and (3) in order to efficiently resolve the dynamics of the gas flow in a network. The idea is that, with the exception of locally high fluctuation, in realistic scenarios, the solution to (2) is on big parts of the network close to the stationary model (3). In these regions we thus can freeze the solution with an acceptable loss in accuracy to save computational effort. By comparison with the solution fully simulated with (2) we then can set up a cost functional measuring the deviation of the partially frozen solution in some appropriate norm. Adding a performance function measuring the time steps where the costly model (2) is calculated, we can set up the optimization problem of weighting the accuracy against the computational effort. Solving this problem enables us to identify a time-dependent model selection for the simulation of gas dynamics on networks that can be used in further examinations to get a cheaply solvable system.
In order to apply adjoint-based gradient methods to solve this problem efficiently, we consider (3) as a limit of appropriately perturbed instationary semilinear dynamics. To this end, we note that, after a linear transformation, system (2) can be written for t∈[0,T] and x∈[0,L], L,T>0, as partially diagonalized system
with the characteristic variable w, see Section 3 for details. Let us assume for the moment that g:R2→R2 is globally Lipschitz-continuous. The same transformation applied to (3), together with given boundary values w1∈C([0,T],R2) yields the ODE-system
To define an appropriate hyperbolic system to approximate (6), for ε>0 and any initial values w0∈L1([0,L],R2) we set up the auxiliary system
Referring to [8,Chapter 3.4] for the definition of a broad solution to system (7), we get the following result:
Lemma 2.1. For any ε>0 there is a unique solution ˉw∈C([0,T]×[0,L],R2) to (6) and a unique broad solution w:[0,T]×[0,L]→R2 to (7). The restriction of w to the set
is a continuous function. Furthermore, for any t0∈(0,T) the function w converges uniformly on [t0,T]×[0,L] for ε↘0 to ˉw, i.e.
Proof. The global existence and uniqueness of ˉw follows from standard ODE-theory. For (t,x)∈Ω, the broad solution w can be written as
Its existence and the continuity on Ω is shown in [8,Theorem 3.1,Theorem 3.4]. Given a fixed t0∈(0,T), we obviously can choose ε small enough such that [t0,T]×[0,L]⊆Ω. We then have
for ε↘0 and all (t,x)∈[t0,T]×[0,L]. The convergence is uniform, since both w12 and g2(w) are continuous and defined on a compact set, thus uniformly continuous. The same argumentation holds for w1.
Retransforming (7) back to the original variables ϱ and q, we get
By Lemma 2.1, this system yields an approximation to model (3).
Now let G=(V,E) be a network, where the type of each node v∈V is given by the parameters α and qv as in Section 2.2. On each edge ej of length Lj we have an initial gas density ϱj0 and gas flow qj0. For an increasing sequence of switching times τ=(τk)k=0,…,N+1⊆[0,∞) and a finite sequence μ=(μk)k=1,…,N, where μk∈{0,1}n denotes for each edge ej∈E, whether model (2) (μk(j)=1) or (8) (μk(j)=0) is used on the time interval (τk−1,τk), consider the PDE-system
where x∈[0,Lj] for j=1,…,n and t∈(τk−1,τk) for k=1,…,N+1. Here,
for a fixed ˉε>0 with ˉε≪1 and
denote the right-hand sides in (8) and (2), respectively. Let ϱj(0,x)=ϱj0(x), qj(0,x)=qj0(x) for x∈[0,Lj] be given initial states. Moreover, at any time point t∈[τ0,τN+1], density and flux additionally satisfy the node coupling conditions
for each v∈V, explained in Section 2.2.
Denote by zd=(ϱ1d,q1d,…,ϱnd,qnd)⊤ the reference solution to (9), (11) for the choice N=1, μ=1 and τ=(0,T), which corresponds to the fine model (2) being fully solved on the complete network and the existence of which will be proven in Section 3. For any other z=(ϱ1,q1,…,ϱn,qn)⊤ then define the cost functional
with γ1,…,γ4≥0, where the first term measures the deviation of z from zd, the second term penalizes using the fine model (2) and the third term penalizes the number of switching time points. Note that, since longer pipes mean more computational effort when using the fine model, the lengths Lj of the pipes enter into the cost as well. For later reference, we set
then J=J1+J2.
The challenge now is to choose the sequences μ and τ such that, with z being the corresponding solution to (9), (11), the cost functional J is minimized. Hence our objective is to solve the minimization problem
3.
Semigroup formulation on networks
In this section, we will set up an abstract formulation of the PDE-system in (9) for τ=(0,T) and any fixed choice of modes per edge and prove that the linear part generates a C0-semigroup. By nonlinear perturbation theory and induction, this yields a well-posedness result of (9) for any finite switching sequence.
For each j∈{1,…,n}, we now consider the initial boundary value problem for zj=(ϱj,qj)⊤ on pipe ej∈E given by
and the coupling conditions
for each node v∈V. Here, Lj>0, zj0:[0,Lj]→R2, fj∈C1(R2,R2) is globally Lipschitz-continuous and
for each j∈{1,…,n}. Moreover, α∈(0,∞)m×2 and qv0:[0,∞)→R is a given outflow for each v∈V. Introduce the space
the vectors
the operators
Moreover, define the operator
on the domain
of all vectors of absolutely continuous functions that satisfy the boundary conditions (16) and (17). Note that, though the inflow functions qv is only evaluated at the origin, it also is shifted in time due to the transport-type evolution represented by operator B - this together enforces the originally time-dependent coupling condition (17). We then can reformulate system (15)-(17) as the abstract initial-value problem
Now we can state the following result:
Theorem 3.1. The operator
is the infinitesimal generator of a strongly continuous semigroup on Z.
Proof. Note that diag(A,B) is a densely defined, closed operator on Z with a nonempty resolvent set (for instance 0∈ρ(diag(A,B))). Following [43,Chapter 4,Theorem 1.3], to prove the claimed semigroup property, it thus suffices to show that the homogeneous system
has a unique solution z∈C1([0,T],Z) for any T>0 and every choice z0∈D(diag(A,B)). So let
be given and first assume T≤ˉT:=min{εjLj2cj|j=1,…,n}>0. We recognise in the operator B a transport equation with velocity −1 for the variables qv1,…,qvm with the well-known unique solution
Next, for each note v∈V we can set
and
Since ϱ1,q1,…,ϱn,qn,qv1,…,qvm are absolutely continuous, so is ϱv for each v∈V. For each edge ej∈E, j=1,…,n, construct the functions ϱj,qj:[0,T]×[0,1]→R as follows: for t∈(0,T] set
Substituting the coupling conditions stated in (19), one can indeed show that, with these definitions,
We skip these calculations here for brevity. Next, set
for (t,x)∈(0,T]×(0,Lj), where
The above considerations show that (ϱj∗,qj∗) is continuous and piecewise absolutely continuous, thus absolutely continuous everywhere. To see that these functions in fact solve system (21), first note that the initial condition and the coupling condition (16) are satisfied by construction. Moreover, for each v∈V we have
Observe that ϱj and qj are continuous at the boundaries x=0 and x=Lj and that
for a.e. (t,x)∈(0,T]×(0,Lj), thus the above construction indeed yields a solution to (21) for T≤ˉT with z(t)∈D for t∈[0,T]. For the case T>ˉT note that the same steps can be repeated successively to expand the solution for the times [ˉT,2ˉT],[2ˉT,3ˉT] and so on. To prove uniqueness of the solution, consider the case z0=0 of initial data and introduce the energy
Obviously, we have E(t)=0 if and only if ϱj(t,.)=qj(t,.)=0 for all j=1,…,n, E(t)≥0 for all t≥0 and E(0)=0. Furthermore,
since ϱv≡0 for all v∈V. But then E≡0, hence z≡0. This concludes the proof for the semigroup property.
We now can apply standard arguments from semigroup theory for further discussion of (20): by [43,Chapter 6,Theorem 1.2 and 1.5], if the inhomogeneity f is continuously differentiable and globally Lipschitz-continuous, (20) has a unique classical solution z. Though these assumptions do not apply to the friction term f stated in (2), we can use the following technical consideration to still get a unique solution: choose any ˉϱ,ˉq>0 and define the modified right-hand side
The function ˜f is continuously differentiable and globally Lipschitz-continuous as a function mapping L2([0,L],R) onto itself for any L>0. If we replace f by the ˜f in (9), we then have a system fitting the assumptions made above for each fixed k=1,…,N and thus we have a unique solution. If we choose ˉϱ sufficiently small and ˉq sufficiently large, then the numerical simulation with realistic data shows that both the boundaries ˉϱ and ˉq will in fact never be reached by ϱ and q, respectively. In this case, however, the solution coincides with the unique solution to the original problem (9). The same argumentation holds for the right-hand side in (8). In summary, this proves the necessary results on well-posedness of the system (20) that we need to examine the optimization problem (14).
4.
Optimality conditions and a solution method
By the results developed in Section 3, we find that the optimization problem (14) is of the form
where Aμk is a strongly continuous semigroup, fμk is a semilinear perturbation and gμk,μk+1=id is a transition map that can be chosen trivially here for each k=1,…,N. The ordering cone T(0,T) is for fixed T>0 defined by
In general, problem (22) does not have a unique global minimum. Instead, we will define a concept of stationarity fitted to the hybrid nature of the problem and set up an algorithm capable of finding such stationary points.
First, however, we apply on system (22) some results from the preliminary work in [45]: by Theorem 3.1 and the subsequent remarks, [45,Lemma 2] can be applied, yielding a control-to-state-map (μ,τ)↦z(μ,τ) that can be substituted into J to define the reduced cost functional Φ(μ,τ)=J(μ,τ,z(μ,τ)). We can now apply [45,Theorem 8] to show that Φ is continuously differentiable with respect to the switching time τk. Recalling the shortened notation introduced in (13), we find that, while J2 can be differentiated directly with respect to (μ,τ), the term J1 depends on the solution z=z(μ,τ). Again by [45,Theorem 8], we can state gradient formulae for the derivatives of J1 based on the adjoint equation
Here, fμz:Z→L(Z) and lz:Z→Z∗ denote the derivatives of the friction term fμ and l as in (13) with respect to z. Substituting the definitions made in (18) yields that in fact γ≡0 and p=(pj)j=1,…,n again can be partitioned into edgewise defined functions pj=(pj1,pj2)⊤ for j=1,…,n satisfying
where
are the derivatives of f0 and f1 in (10), as well as the coupling conditions
By [45,Lemma 6], system (23) (and thus (24), (25), (26) as well) have a unique mild solution and, applying [45,Theorem 8], the switching time gradient is given by
If the mode sequence μ is fixed, then we can conclude that a switching sequence τ=(τk)k=0,…,N+1 is a KKT-point of the minimization problem (14), if the following holds: For n∈{1,…,N} set a(τ,n)=min{m∈{0,…,n}|τm=τn} and b(τ,n)=max{m∈{n,…,N+1}|τm=τn}, then
Similarly, by [45,Theorem 10], the sensitivity of the cost function with respect to introducing a new mode μ′ on an infinitesimal time interval at the point τ′ can be represented by the mode insertion gradient given by
where μk is the original mode at time τ′. In summary, a switching sequence (μ,τ) is called stationary, if τ is a KKT-point for μ fixed and if ∂Φ∂μ′(τ′)≥0 for all modes μ′ and all times τ′∈[0,T]. Any global minimum of the problem (14) then is stationary.
In order to compute such stationary switching signals, we consider a conceptual algorithm originally proposed for optimal mode scheduling in hybrid ODE-dynamical systems [2]. The main idea is a two phase approach as follows: the algorithm is initialized with a switching sequence (μ0,τ0), for instance μ0=1 and τ0=(0,T) where no switching occurs and the system is solved by keeping the mode constant at 1.
In a first phase, the positions of available switching time points are optimized, while conserving their order, by using a projected-gradient method with Armijo step size. To this end a projection P onto the ordering cone T(0,T) is used. In a second phase, a new mode μ′ is inserted at a time point τ′ where ∂Φ∂μ′(τ′)<0. If no such point exists, the switching sequence is stationary in the above sense, otherwise the algorithm continues with the first phase again.
A more precise description of the procedure is given in Algorithm 1. Note that Algorithm 1 does not necessarily terminate with the global solution to the minimization problem (22), but with a stationary point in the above sense. A convergence analysis for this algorithm can be found in [2]. A possible implementation of the projection P can be found in [19].
5.
Numerical study
As a proof of concept we consider a gas network outlined in Figure 1b. At the boundary node N1 gas is supplied to the network whereas nodes N2 and N3 are customer nodes where gas can possibly be taken out of the network. This can be seen as a simplified example of a big regional gas pipeline network with a local distribution subnetwork. In the particular scenario we are looking at, there is gas transported from the supply N1 to node N2 to satisfy a given demand while N3 is inactive. All pipes are assumed to be horizontal (h′=0), the compressor is assumed to be running at a constant compression factor ˉα, compare (5), and at initial time the network is assumed to be stationary with constant densities ϱ0 on the outer circle (thus αϱ0 on the inner circle) and flux q0 everywhere, moreover the outflows at the nodes N1 and N2 are outlined in Figure 1c. See the table in Figure 1a for specific values for those and other constants.
Due to the almost decoupled inner and outer circle, it can be expected that the numerical solution highly varies on the outer circle connecting N1 and N2 but is near to constant on the inner circle. This is confirmed by the simulation of the full model, see Figure 2 for a snapshot of the simulation at the time t=900 s. We therefore can suspect that the simulation on the inner circle can be widely frozen with only some small losses in the accuracy of the solution. Indeed, letting ˉz be the distinguished solution to (9), (11) resulting from freezing the solution completely on the inner edges 6 to 10, our simulations show that the L2-errors of the density and the flux relative to the respective maximum values in zd is less than 1% for ϱ and less than 6.5% for q; see Figure 3d. Moreover, compared to a full simulation zd about half of the computation time in the sense of J2 defined in (13) is saved.
In our implementation the system (9), (11) is solved via splitting of the hyperbolic part and the friction term. For the simulation of the hyperbolic part we use the 2-step-Richtmyer-method with artificial viscosity, see [32,Chapter 18.1]. Given a system matrix A and a discretization zk=(zk1,…,zkn)⊤ with spatial step size Δx of the solution at time point tk, it computes the discretized values zk+1 at time tk+1=tk+Δt by the explicit finite-volume-scheme
for j=1,…,n, where ε∼0.05 introduces an additional, minor but stabilizing smoothing to the solution. The discretization is chosen in a way such that the cells zk1 and zkn are centered at the boundary points x=0 and x=1, respectively. In order to handle the boundary values and the inner cells simultaneously with the same scheme, we add appropriate ghost cells on each end of the spatial domain. The appropriate ghost cell values (ze)k0,(ze)kn+1 for each e∈E are derived from the coupling conditions. We mention without further proof that these can be solved explicitly for each node v∈V by first setting the weighted mean value
and then using the ghost cell values
Here, in each time step, we incorporate only those edges where the solution is actively calculated. In the special case if v∈V is a compressor, we only have one ingoing edge e+∈δ+v and one outgoing edge e−∈δ−v and instead have to set
For the friction term, we add an explicit Runge-Kutta-step using the classical Runge-Kutta-scheme of order 2 or midpoint rule, see [10,Chapter 23]. The same methods are used on the same discretization grid backwards in time to solve the adjoint system stated in (24), (25), (26).
To evaluate the gradient formulae (27) and (29), we use the trapezoidal rule over the spatial grid, see [14,Chapter 2.1]. The expression [(Aμk)j−(Aμk−1)j]zj(τk,x) occurring in both (27) and (29) represents the difference of the time derivatives of the solution depending on which mode μk is switched to. In the numerical scheme, this is realized by calculating a time step of the solution for each choice of μk and then substracting the forward difference quotients.
The fixed temporal discretization grid for the actual solution is supplemented by a grid of switching time points that may vary due to the superordinated optimization where the data needed for the gradient formulae is calculated. For our study we start the optimization with the fully frozen solution, where in no time step the model (2) is actually calculated, and iteratively insert regions of active numerical solving wherever the gradients indicate a major loss of accuracy. Applying the projected-gradient method with sequential mode insertion described in Algorithm 1 yields the results shown in Figure 3a. We observe that it is in fact almost unnecessary to calculate the fine model on the edges 6 to 10 of the inner circle and that the algorithm indeed approximates the distinguished solution ˉz as expected. Note that Algorithm 1 does not remove switching points during its process, because doing so might lead to the iteration being caught in an infinite loop, where the same switching points are added and removed repeatedly. This way, however, there remain scattered short intervals where the mode is switched twice almost instantly, which can be interpreted as the algorithm eliminating the respective intermediate modes. We suggest a post-processing step to remove this scattering, for instance by applying the following filtering rule with d∈N a fixed parameter: if within d time steps after a switching point in the numerical solution the mode is switched again, then both switching time points are removed. Relative to the initial cost, we get for the optimized solution a cost reduction of 96.4% with 65 switching points, for d=5 shown in Figure 3b a cost reduction of 96.3% with 19 switching points and for d=10 shown in Figure 3c a cost reduction of 95.1% with 7 switching points.
6.
Conclusion
We have presented an application of the theory of switching systems to a model hierarchy for the dynamics of gas in a pipeline network. A semigroup formulation was given for the model on gas networks including time-dependent outflow at each node as well as a linearized model for compressors that allowed us to prove the existence of unique solutions. Using adjoint-based gradient representations for switching abstract evolution systems, we implemented a two stage projected-gradient descent method for the optimization of switching between different levels of accuracy in the hierarchy. As an example, we optimized the simulation of a small supply network by freezing the calculation on edges whenever the numerical error made is small compared to the computational costs.
Our prototypical approach can be applied in a similar fashion to realistic industrial networks. The technique can also be extended to identify a model switching strategy for a reduced model using a range of different parameter configurations such as representing different boundary flow scenarios, compressor settings and valve positions. The resulting reduced simulation model can then be used for a time expanded mixed-integer optimization technique based on full discretization with a minimum of variables or within a bilevel optimisation method which optimises a cost functional on the outer level with optimal efficiency on the lower level as proposed, for example, in [24]. Further directions include a combination of this method in a network of submodel hierarchies such as those used in [37] for the optimisation of stationary models. Moreover, it remains an open question, whether the chosen numerical method lead to a consistent discretization of the optimality system. Well-balanced finite-volume schemes also seem to be an promising alternative approach in this context.
Acknowledgments
The authors thank the Deutsche Forschungsgemeinschaft for their support within the projects A03 and B03 in the Sonderforschungsbereich/-Transregio 154 Mathematical Modelling, Simulation and Optimization using the Example of Gas Networks. Furthermore, the authors are indebted to the anonymous reviewers who considerably improved the exposition of this paper.