- Email: [email protected]
A joint IMEX-MOR approach for Water Networks
8th Vienna Vienna International International Conference Conference on on Mathematical Mathematical Modelling Modelling 8th 8th Vienna International Conference on Mathematical Modelling February -- 20, University of 8th Vienna18 Conference on Mathematical Available onlineModelling at Vienna, www.sciencedirect.com February 18International 20, 2015. 2015. Vienna Vienna University of Technology, Technology, Vienna, February Austria February 18 18 -- 20, 20, 2015. 2015. Vienna Vienna University University of of Technology, Technology, Vienna, Vienna, Austria Austria Austria
ScienceDirect
IFAC-PapersOnLine 48-1 (2015) 260–261
A joint A joint IMEX-MOR IMEX-MOR approach approach for for Water Water Networks Networks A A joint joint IMEX-MOR IMEX-MOR approach approach for for Water Water Networks Networks ∗∗ Lennart Jansen ∗∗ ∗∗ Sara Grundel Sara Grundel ∗ Lennart Jansen ∗∗ Sara Grundel Sara Grundel ∗ Lennart Lennart Jansen Jansen ∗∗ ∗∗ Max-Planck Institute for Dynamics of Complex Technical Systems, Institute for Dynamics of Complex Technical Systems, ∗∗ Max-Planck Max-Planck Institute Dynamics of Max-PlanckGermany Institute for for Dynamics of Complex Complex Technical Technical Systems, Systems, Magdeburg, Germany (e-mail: [email protected]). Magdeburg, (e-mail: [email protected]). ∗∗ Magdeburg, Germany (e-mail: ∗∗ Heinrich Magdeburg, Germany (e-mail: [email protected]). [email protected]). Heine Univerist¨ tt D¨ uusseldorf, Germany (e-mail: Heinrich Heine Univerist¨ D¨ sseldorf, Germany (e-mail: ∗∗ ∗∗ Heinrich Univerist¨ Heinrich Heine Heine Univerist¨aatt D¨ D¨uusseldorf, sseldorf, Germany Germany (e-mail: (e-mail: [email protected]) [email protected]) [email protected]) [email protected])
Abstract: Modeling and Simulation of fluids large network is aa rather problem. We Abstract: Modeling and Simulation of fluids in in large network is rather challenging challenging problem. We Abstract: Modeling and Simulation of in large network is challenging problem. Abstract: Modeling and Simulation of fluids fluids in large Order network is aa rather rather challenging problem. We We provide an approach combining techniques in Model Reduction (MOR) and implicit-explicit provide an approach combining techniques in Model Order Reduction (MOR) and implicit-explicit provide an an approachto combining combining techniques in Model Model Order Reduction Reduction (MOR) (MOR) and and implicit-explicit implicit-explicit provide approach techniques in Order (IMEX) integration create efficient and stable simulations. (IMEX) integration to create efficient and stable simulations. (IMEX) (IMEX) integration integration to to create create efficient efficient and and stable stable simulations. simulations. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Model Reduction, Integration, Differential Equations, Networks, Directed Graphs Keywords: Model Reduction, Integration, Differential Equations, Networks, Directed Graphs Keywords: Model Model Reduction, Reduction, Integration, Integration, Differential Differential Equations, Equations, Networks, Networks, Directed Directed Graphs Graphs Keywords: 1. INTRODUCTION A demand node has demand . Thus, the + 1. INTRODUCTION A demand node has aaa given given demand qqqss ::: III −→ −→ R R the + .. Thus, 1. INTRODUCTION INTRODUCTION A demand node given R Thus, the ss : I −→ 1. A demandbetween node has has given demand demand −→towards R+ the difference between theaamount amount of water waterqflowing flowing towards a node, node, + . Thus, difference the of a difference between the amount of water flowing towards a node, difference between the amount of water flowing towards a node, and the amount of water flowing away from the node has to be The simulation of fluids within aa large network of pipes poses and the amount of water flowing away from the node has to be The simulation of fluids within large network of pipes poses and the amount of water flowing away from the node has to be The simulation simulation of fluids fluids within aa large large network of pipes pipes poses and the amount of water flowing away from the node has to be q The of within network of poses several mathematical challenges. Typically after spatial disqq several mathematical challenges. Typically after spatial disseveral mathematical challenges. Typically after spatial disq several mathematical challenges. Typically after spatial discretization the resulting mathematical system is a nonlinear difcretization the resulting mathematical system is nonlinear difcretization the resulting resulting mathematical system is is aaaare nonlinear difcretization the mathematical system nonlinear differential algebraic system. Standard techniques often slow qqiiiin − qqiiiout = q ∑ ∑ ferential algebraic system. Standard techniques are often slow − qss (t) (t) ∑ ∑ out = in i ferential algebraic system. Standard techniques are often slow q − = i∈I i∈I ∑ ∑ out in ferential algebraic system. StandardWe techniques are often slow due to the stiffness of the equation. will show aa several step − i∈I∑out qqiout =q qss (t) (t) i∈I ∑ in qin out due to the stiffness of the equation. We will show several step in i∈I i∈I out in due to to the the stiffness of the equation. equation. We will show show several step i∈Iin i∈Iout due stiffness of the We will aa several step process on how to improve on the timing. A first and major process on how to improve on the timing. A first and major process on how how to improve on the the timing. A first first and and major process on to improve on timing. A major qqiiiout being incoming and outgoing flow of edges with qqiiiin and step in order to achieve stable and fast simulators for these and being the the incoming and outgoing flow of edges with step in order to achieve stable and fast simulators for these out in i iout being and qqthe the incoming and outgoing flow of edges with q step in order to achieve stable and fast simulators for these and being the incoming and outgoing flow of edges with q in step in order to achieve stable and fast simulators for these connected to demand node, respectively. It is possible that problems is what we call the decoupling step. In that step we out in connected to the demand node, respectively. It is possible that problems is what we call the decoupling step. In that step we connected to the demand node, respectively. It is possible that problems is what we call the decoupling step. In that step we connected to the demand node, respectively. It is possible that problems ismodel what we call the as decoupling step. In 11that stepThis we qqss (t) (t) = = 00 for for all all tt ∈ ∈ I. I. are able to the system aa discrete index DAE. are able to model the system as discrete index DAE. This qss (t) (t) = = 00 for for all all tt ∈ ∈ I. I. are able able to model model thedue system aschoosen a discrete discrete index 11 DAE. DAE. This q are to the system as a index This step is only possible to the discretization we use. step is only possible due to the choosen discretization we use. step is isweonly only possible due to to of theModel choosen discretization we use. step possible due the choosen discretization we use. Next use aa combination Order Reduction (MOR) Next we use combination of Model Order Reduction (MOR) 2.2 Edge Elements Edge Elements Next we use a combination of Model Order Reduction (MOR) Next we use a combination Modelscale Orderindex Reduction (MOR) 2.2 methods in order to create aaofsmaller 11 Differential 2.2 Edge Elements methods in order to create smaller scale index Differential 2.2 Edge Elements methods in order to create a smaller scale index 1 Differential methods in order to create a smaller scale index 1 Differential Algebraic Equation. And last but not least we use an implicitAlgebraic Equation. And last but not least we use an implicitFirst we will discuss the pipe model. The behavior of aa pipe Algebraic Equation. And last last but not nottoleast least wethe usetime-step an implicitimplicitFirst we will discuss the pipe model. The behavior of Algebraic Equation. And but we use an explicit (IMEX) integration method reduce for explicit (IMEX) integration method to reduce the time-step for First we will discuss the pipe model. The behavior of aa pipe pipe First we will discuss the pipe model. The behavior of pipe is described by a continuity equation and an equation describexplicit (IMEX) integration method to reduce the time-step for is described by a continuity equation and an equation describexplicit (IMEX) integration method to reduce the time-step for the stiff nonlinear differential equation. We will only present is described by a continuity equation and an equation describthe stiff nonlinear differential equation. We will only present is described by a continuity equation and an equation describing the movement inside a pipe. We consider circular pipes the stiff nonlinear differential equation. We will only present ing the movement inside a pipe. We consider circular pipes the stiff nonlinear differential equation. We willreservoirs only present aa simplified network here which includes pipes, and π simplified network here which includes pipes, reservoirs and 2 ing the movement inside We circular pipes πconsider ing the movement inside aa pipe. pipe.A We consider circular�. D The with diameter D, cross-section = a simplified simplified network hereThis which includes pipes, reservoirs and D2 and and length length �.pipes The with diameter D, cross-section A = π aso network here which includes pipes, reservoirs and called demand nodes. system will actually result in an π44 D22 and length �. The so called demand nodes. This system will actually result in an with diameter D, cross-section A = D and length �. The with diameter D, cross-section A = independent variables are space x ∈ [0, �] := Ω ⊂ R and time 4 so called demand nodes. This system will actually result in an independent variables are space x ∈ [0, 4 �] := Ω ⊂ R and time so called demand nodes. This system will actually result in an ODE which simplifies the discussion. ODE which simplifies the discussion. independent variables are space xx ∈ [0, �] := Ω ⊂ R and time independent variables are space ∈ [0, �] := Ω ⊂ R and time t ∈ [t , T ] := I ⊂ R. The time dependent variables are the mass ODE which simplifies the discussion. 0 tt ∈ [t0 ,, T ] := I ⊂ R. The time dependent variables are the mass ODE which simplifies the discussion. ∈ R. time tflow ∈ [t [t00m T:: ]]Ω :=× ⊂−→ R. The The time dependent variables are the mass flow m, T Ω:= ×IIII⊂ −→ R and and thedependent pressure variables p :: Ω Ω× × II are −→the R.mass The R the pressure p −→ R. The 2. MODELING flow m :: Ω Ω× × −→ Rdepend and the theonpressure pressure Ω× × II −→ −→ R. pipe The 2. MODELING flow m II −→ and ppproperties :: Ω R. The parameters a, ρ and ccR material of the parameters a, ρ and depend on material properties of the pipe 2. MODELING 2. MODELING parameters a, ρ and ccangle depend on material properties of the pipe parameters a, ρ and depend on material properties of the pipe and the gas. α is the of the pipe and g is the gravitational and the gas. α is the angle of the and g is the gravitational and the α the of the pipe pipe and the It is common to define aa connected and directed simple graph and the gas. gas. α is isthis, the angle angle of the pipe and gg is is the gravitational gravitational It is common to define connected and directed simple graph constant. With we get following partial differential constant. With this, we get the following partial differential It is common to define a connected and directed simple graph It is common to define a connected and directed simpleaagraph constant. With this, we get the following partial differential G = G(V, E) representing the pipe network. This allows more constant. With this, we get the following partial differential G = G(V, E) representing the pipe network. This allows more equation, which describes the behavior in pipes equation, which describes the behavior in pipes G = G(V, E) representing the pipe network. This allows a more G = G(V,representation E) representing the pipe network. This allows aV more equation, which describes the behavior in pipes compact of the model equations. The set are equation, which describes the behavior in pipes compact representation of the model equations. The set V are 2 2 m compact representation of the theand model equations. The set V are are ∂∂∂ ppp (x,t) + aaa22 ∂∂∂ m compact representation of model equations. The set V the nodes and E are the edges we will describe the different the nodes and E are the edges and we will describe the different (x,t) = 0 + = 0 ∂∂tp (x,t) aA ∂∂m m the nodes nodes and E Eelements are the the edges edges and we will will describe describe the the different different xx (x,t) the and are and we node and edge in the following. (x,t) + (x,t) = 0 ∂t A ∂ node and edge elements in the following. (x,t) + (x,t) = 0 (1) (1) ∂t A ∂ x node and edge elements in the following. ∂ pp m ∂t A ∂x node and edge elements in the following. (1) ∂ ∂∂∂ m (1) (x,t) + A (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = ∂∂ pp (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = 0 m (x,t) + A 0 ∂ m ∂t ∂ x 2.1 Node Elements (x,t) + A (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = 0 ∂t 2.1 Node Elements (x,t) + A ∂∂ xx (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = 0 ∂t 2.1 Node Node Elements Elements ∂t ∂x 2.1 Hyperbolic PDE pipe: Reservoirs are water sources with unlimited capacity. Thus, we Hyperbolic PDE pipe: Reservoirs are water sources with unlimited capacity. Thus, we Hyperbolic PDE pipe: Reservoirs are water sources with unlimited capacity. Thus, we Hyperbolic PDE pipe: Reservoirs are water sources with pressure unlimitedppcapacity. Thus, we Furthermore no assume that they have aa constant ss .. Furthermore assume that they have constant pressure no Further edge edge elements elements are are valves valves and and pumps, pumps, which which we we also also assume that that they holds have aaatconstant constant pressure pss ..arbitrary Furthermore no Further Furthermore no assume they have pressure p balance equation a reservoir, since an amount balance equation holds at aa reservoir, since an arbitrary amount Further edge elements are Further edgeextended elementsabstract. are valves valves and and pumps, pumps, which which we we also also omit in this balance equation holds at reservoir, since an arbitrary amount omit in this extended abstract. balance equation holds at a reservoir, since an arbitrary amount of water may leave or enter the reservoir. of water may leave or enter the reservoir. omit in this extended abstract. omit in this extended abstract. of water may leave or enter the reservoir. of water may leave or enter the reservoir. 2.3 Network Model 2.3 Network Model p = p 2.3 Network Model p = p ss 2.3 Network Model p = p p = pss From now on we aa network nn p many pipes and From now on we consider consider network with with pipes and now on consider with nn ppp many many pipes and From nowdemand on we we nodes consider network with many pipes and ndd many In contrast contrast to to reservoirs, reservoirs, tanks tanks have have limited limited capacity. capacity. Never Never the the From many demand nodes andaannnetwork many reservoirs. For each pipe rs n In and many reservoirs. For each pipe rs ni,ddwe In contrast to reservoirs, reservoirs, tanks have have limited capacity. Never the many demand nodes and n many reservoirs. For each pipe rs n In contrast to tanks limited capacity. Never the many demand nodes and n many reservoirs. For each pipe less, pressure can inor decrease even though the tank is full or get a flow m : Ω × I −→ R, Ω = [x , x ] ⊂ R, I = [t , T ]] ii : Ωii × I −→rsR, Ωii = [xL R 0 less, pressure can inor decrease even though the tank is full or i, we get a flow m , x ] ⊂ R, I = [t i i Li , xRi ] ⊂ R, I = [t0 ,, T less, pressure can inor decrease even though the tank is full or i, we get aa flow m : Ω × I −→ R, Ω = [x T i i i L R 0 less, pressure can inor decrease even though the tank is full or i, we get flow m : Ω × I −→ R, Ω = [x , x ] ⊂ R, I = [t T ]]ii i i empty respectively. We will not talk about tanks in more details and a pressure p : Ω ×I −→ R both depending on space x ∈ i −→ R bothi depending Li Ri on space x ∈ 0,Ω empty respectively. We will not talk about tanks in more details and aa pressure ppii ::i Ω Ω ii ×I n empty respectively. We will not talk about tanks in more details and pressure Ω ×I −→ R both depending on space x ∈ Ω d empty respectively. We will not talk about tanks in more details and atime pressure pii :The Ωii ×I −→variables R both depending space Ωii here. tt ∈ node are R here. and time ∈ I. I. The node variables are pppdd ::: IIIon−→ −→ Rnnnxdd ∈and and here. and here. and time time tt ∈ ∈ I. I. The The node node variables variables are are pdd : I −→ −→ R R d and and
2405-8963 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright © 2015, IFAC 260 Copyright 2015,responsibility IFAC 260Control. Peer review© under of International Federation of Automatic Copyright © 2015, IFAC 260 Copyright © 2015, IFAC 260 10.1016/j.ifacol.2015.05.201
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria
Sara Grundel et al. / IFAC-PapersOnLine 48-1 (2015) 260–261
prs : I −→ Rnrs , the pressures at demand nodes and reservoirs. We define mL (t) = (mi (xLi ,t))i∈A pi ,mR (t) = (mi (xRi ,t))i∈A pi . pL (t) = (pi (xLi ,t))i∈A pi ,pR (t) = (pi (xRi ,t))i∈A pi . mL being the vector with all pipe flows at their tail-node mR the flow vector at their head-nodes and similarly for pL pR . Note, that the node pressures coincide with the headtail pressures of the pipes. We call the vector of demand reservoir pressures by p.
and and and and
Last we define the following incidence matrices 1 if reservoir node j is head of pipe i n p ×nrs rs ∈ R (A ) = Ars R R ij 0 else −1 if reservoir node j is tail of pipe i n p ×nrs rs Ars (AL )i j = L ∈R 0 else 1 if demand node j is head of pipe i AR ∈ Rn p ×nd ,(AR )i j = 0 else −1 if demand node j is tail of flow i AL ∈ Rn p ×nd ,(AL )i j = 0 else We can combine them all and get the full incidence matrix A
n p ×+drs +nd rs A := (Ars R + AL AR + AL ) ∈ R With the help of these matrices we can write the spatial discretized system of equations as (2) p˙R + Dα (mL − mR ) = 0
m˙ L + Dβ AT p + γ + G(mL )mL = 0 (3) (4) AR mR + AL mL = qset (5) AR pd = pR AL pd = pL (6) p = p (7) Ars R R rs rs AL prs = pL (8) (9) prs = pset To obtain these equation it is crucial to chose a suitable spatial discretization. In particular the time derivative of the pressure is evaluated at the right end of the pipe and the time derivative of the flux at the left end. The size of the first two equation is the number of pipes, equation (4) the number of junction, (5,6,7,8) number of pipes and the last equation number of tanks. Dα is a diagonal matrix containing αi = a2i /Ai /�i on the diagonal and Dβ similar with βi = Ai /�i . γ is a vector with γi = ρi Ai g sin αi and G(mL ) is a diagonal matrix function such that G(mL )i = ci miL . The matrix A is the incidence matrix as described above. The vector qset has an entry for every demand node showing the given demand at that particular node given by qs (t) and similar is pset the vector of the given pressures ps at the reservoirs. In the modeling of the graph it is crucial to pick the direction of the edges such that every demand node has a right end of a pipe. This is possible for any topology as long as one of the nodes in the graph is not a demand node which means in our case it has to be a reservoir. We furthermore want all pipes the end in a reservoir to end in a left node there. 3. DECOUPLING By selecting a matrix Aselect which picks one pipe for each node that has a right end in that given node we can rewrite the system of equation in the variables pd and mL .
261
261
p˙d + Aselect Dα (−Cqset +CAmL ) = 0 pd T + γ + G(mL )mL = 0 m˙ L + Dβ A pset
We will explain in detail how we create the matrix C which is the crucial part in this decoupling process. This resulting ODE is of size nd + n p and has the general structure x˙ = T x + g(x,t) = f (x,t),
(10)
where the matrix T is given by 0 Aselect Dal phaCA T= , Dβ (AR + AL )T 0 and the vector x is combined by pd and mL . 4. IMEX In order to solve this stiff and nonlinear ODE we make use of implicit-explicit (IMEX) integration methods. This allow us to deal with the stiffness in an efficient way while not having to solve large-scale nonlinear problems. First order methods are of the flavor xn+1 − xn (11) = (1 − γ)T xn + γT xn+1 + g(xn ,t) h for γ ∈ [0, 1] and time setp h, which leads to the iteration
xn+1 = (1 − hγT )−1 (xn + (1 − γ)T xn + hg(xn ). We study convergence properties of that by analyzing the matrix T and the function g as well as the differences for several values of γ. If γ = 0 we get explicit Euler and if γ = 1 we get a combination of implicit Euler for the linear part an explicit Euler for the nonlinear part. We will also show the differences within this class of methods as well as the difference to second order methods following Ascher et al. (1995). 5. MODEL ORDER REDUCTION On the resulting ODE (10) we use the Model Order Reduction techniques Proper Orthogonal Decomposition (POD) together with Discrete Empirical Interpolation (DEIM) by Chaturantabut and Sorensen (2010). POD is a projection-based method where we find a projection matrix W such that the solution of (10) x ≈ W xˆ for xˆ in a lower dimensional space. The resulting low-dimensional ODE is then given by ˆ xˆ˙ = W T TW xˆ +W T g(W x,t). DEIM is then used to create a truely low-dimensional function approximating W T g(W x,t). ˆ 6. CONCLUSIONS The combination of POD-DEIM with the IMEX integration results in a significant speedup of simulation time. REFERENCES Ascher, U.M., Ruuth, S.J., and Wetton, B.T. (1995). Implicitexplicit methods for time-dependent partial differential equations. SIAM Journal on Numerical Analysis, 32(3), 797–823. Chaturantabut, S. and Sorensen, D.C. (2010). Nonlinear model reduction via discrete empirical interpolation. SIAM J. Sci. Comput., 32(5), 2737–2764. doi:10.1137/090766498.
ScienceDirect
IFAC-PapersOnLine 48-1 (2015) 260–261
A joint A joint IMEX-MOR IMEX-MOR approach approach for for Water Water Networks Networks A A joint joint IMEX-MOR IMEX-MOR approach approach for for Water Water Networks Networks ∗∗ Lennart Jansen ∗∗ ∗∗ Sara Grundel Sara Grundel ∗ Lennart Jansen ∗∗ Sara Grundel Sara Grundel ∗ Lennart Lennart Jansen Jansen ∗∗ ∗∗ Max-Planck Institute for Dynamics of Complex Technical Systems, Institute for Dynamics of Complex Technical Systems, ∗∗ Max-Planck Max-Planck Institute Dynamics of Max-PlanckGermany Institute for for Dynamics of Complex Complex Technical Technical Systems, Systems, Magdeburg, Germany (e-mail: [email protected]). Magdeburg, (e-mail: [email protected]). ∗∗ Magdeburg, Germany (e-mail: ∗∗ Heinrich Magdeburg, Germany (e-mail: [email protected]). [email protected]). Heine Univerist¨ tt D¨ uusseldorf, Germany (e-mail: Heinrich Heine Univerist¨ D¨ sseldorf, Germany (e-mail: ∗∗ ∗∗ Heinrich Univerist¨ Heinrich Heine Heine Univerist¨aatt D¨ D¨uusseldorf, sseldorf, Germany Germany (e-mail: (e-mail: [email protected]) [email protected]) [email protected]) [email protected])
Abstract: Modeling and Simulation of fluids large network is aa rather problem. We Abstract: Modeling and Simulation of fluids in in large network is rather challenging challenging problem. We Abstract: Modeling and Simulation of in large network is challenging problem. Abstract: Modeling and Simulation of fluids fluids in large Order network is aa rather rather challenging problem. We We provide an approach combining techniques in Model Reduction (MOR) and implicit-explicit provide an approach combining techniques in Model Order Reduction (MOR) and implicit-explicit provide an an approachto combining combining techniques in Model Model Order Reduction Reduction (MOR) (MOR) and and implicit-explicit implicit-explicit provide approach techniques in Order (IMEX) integration create efficient and stable simulations. (IMEX) integration to create efficient and stable simulations. (IMEX) (IMEX) integration integration to to create create efficient efficient and and stable stable simulations. simulations. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Model Reduction, Integration, Differential Equations, Networks, Directed Graphs Keywords: Model Reduction, Integration, Differential Equations, Networks, Directed Graphs Keywords: Model Model Reduction, Reduction, Integration, Integration, Differential Differential Equations, Equations, Networks, Networks, Directed Directed Graphs Graphs Keywords: 1. INTRODUCTION A demand node has demand . Thus, the + 1. INTRODUCTION A demand node has aaa given given demand qqqss ::: III −→ −→ R R the + .. Thus, 1. INTRODUCTION INTRODUCTION A demand node given R Thus, the ss : I −→ 1. A demandbetween node has has given demand demand −→towards R+ the difference between theaamount amount of water waterqflowing flowing towards a node, node, + . Thus, difference the of a difference between the amount of water flowing towards a node, difference between the amount of water flowing towards a node, and the amount of water flowing away from the node has to be The simulation of fluids within aa large network of pipes poses and the amount of water flowing away from the node has to be The simulation of fluids within large network of pipes poses and the amount of water flowing away from the node has to be The simulation simulation of fluids fluids within aa large large network of pipes pipes poses and the amount of water flowing away from the node has to be q The of within network of poses several mathematical challenges. Typically after spatial disqq several mathematical challenges. Typically after spatial disseveral mathematical challenges. Typically after spatial disq several mathematical challenges. Typically after spatial discretization the resulting mathematical system is a nonlinear difcretization the resulting mathematical system is nonlinear difcretization the resulting resulting mathematical system is is aaaare nonlinear difcretization the mathematical system nonlinear differential algebraic system. Standard techniques often slow qqiiiin − qqiiiout = q ∑ ∑ ferential algebraic system. Standard techniques are often slow − qss (t) (t) ∑ ∑ out = in i ferential algebraic system. Standard techniques are often slow q − = i∈I i∈I ∑ ∑ out in ferential algebraic system. StandardWe techniques are often slow due to the stiffness of the equation. will show aa several step − i∈I∑out qqiout =q qss (t) (t) i∈I ∑ in qin out due to the stiffness of the equation. We will show several step in i∈I i∈I out in due to to the the stiffness of the equation. equation. We will show show several step i∈Iin i∈Iout due stiffness of the We will aa several step process on how to improve on the timing. A first and major process on how to improve on the timing. A first and major process on how how to improve on the the timing. A first first and and major process on to improve on timing. A major qqiiiout being incoming and outgoing flow of edges with qqiiiin and step in order to achieve stable and fast simulators for these and being the the incoming and outgoing flow of edges with step in order to achieve stable and fast simulators for these out in i iout being and qqthe the incoming and outgoing flow of edges with q step in order to achieve stable and fast simulators for these and being the incoming and outgoing flow of edges with q in step in order to achieve stable and fast simulators for these connected to demand node, respectively. It is possible that problems is what we call the decoupling step. In that step we out in connected to the demand node, respectively. It is possible that problems is what we call the decoupling step. In that step we connected to the demand node, respectively. It is possible that problems is what we call the decoupling step. In that step we connected to the demand node, respectively. It is possible that problems ismodel what we call the as decoupling step. In 11that stepThis we qqss (t) (t) = = 00 for for all all tt ∈ ∈ I. I. are able to the system aa discrete index DAE. are able to model the system as discrete index DAE. This qss (t) (t) = = 00 for for all all tt ∈ ∈ I. I. are able able to model model thedue system aschoosen a discrete discrete index 11 DAE. DAE. This q are to the system as a index This step is only possible to the discretization we use. step is only possible due to the choosen discretization we use. step is isweonly only possible due to to of theModel choosen discretization we use. step possible due the choosen discretization we use. Next use aa combination Order Reduction (MOR) Next we use combination of Model Order Reduction (MOR) 2.2 Edge Elements Edge Elements Next we use a combination of Model Order Reduction (MOR) Next we use a combination Modelscale Orderindex Reduction (MOR) 2.2 methods in order to create aaofsmaller 11 Differential 2.2 Edge Elements methods in order to create smaller scale index Differential 2.2 Edge Elements methods in order to create a smaller scale index 1 Differential methods in order to create a smaller scale index 1 Differential Algebraic Equation. And last but not least we use an implicitAlgebraic Equation. And last but not least we use an implicitFirst we will discuss the pipe model. The behavior of aa pipe Algebraic Equation. And last last but not nottoleast least wethe usetime-step an implicitimplicitFirst we will discuss the pipe model. The behavior of Algebraic Equation. And but we use an explicit (IMEX) integration method reduce for explicit (IMEX) integration method to reduce the time-step for First we will discuss the pipe model. The behavior of aa pipe pipe First we will discuss the pipe model. The behavior of pipe is described by a continuity equation and an equation describexplicit (IMEX) integration method to reduce the time-step for is described by a continuity equation and an equation describexplicit (IMEX) integration method to reduce the time-step for the stiff nonlinear differential equation. We will only present is described by a continuity equation and an equation describthe stiff nonlinear differential equation. We will only present is described by a continuity equation and an equation describing the movement inside a pipe. We consider circular pipes the stiff nonlinear differential equation. We will only present ing the movement inside a pipe. We consider circular pipes the stiff nonlinear differential equation. We willreservoirs only present aa simplified network here which includes pipes, and π simplified network here which includes pipes, reservoirs and 2 ing the movement inside We circular pipes πconsider ing the movement inside aa pipe. pipe.A We consider circular�. D The with diameter D, cross-section = a simplified simplified network hereThis which includes pipes, reservoirs and D2 and and length length �.pipes The with diameter D, cross-section A = π aso network here which includes pipes, reservoirs and called demand nodes. system will actually result in an π44 D22 and length �. The so called demand nodes. This system will actually result in an with diameter D, cross-section A = D and length �. The with diameter D, cross-section A = independent variables are space x ∈ [0, �] := Ω ⊂ R and time 4 so called demand nodes. This system will actually result in an independent variables are space x ∈ [0, 4 �] := Ω ⊂ R and time so called demand nodes. This system will actually result in an ODE which simplifies the discussion. ODE which simplifies the discussion. independent variables are space xx ∈ [0, �] := Ω ⊂ R and time independent variables are space ∈ [0, �] := Ω ⊂ R and time t ∈ [t , T ] := I ⊂ R. The time dependent variables are the mass ODE which simplifies the discussion. 0 tt ∈ [t0 ,, T ] := I ⊂ R. The time dependent variables are the mass ODE which simplifies the discussion. ∈ R. time tflow ∈ [t [t00m T:: ]]Ω :=× ⊂−→ R. The The time dependent variables are the mass flow m, T Ω:= ×IIII⊂ −→ R and and thedependent pressure variables p :: Ω Ω× × II are −→the R.mass The R the pressure p −→ R. The 2. MODELING flow m :: Ω Ω× × −→ Rdepend and the theonpressure pressure Ω× × II −→ −→ R. pipe The 2. MODELING flow m II −→ and ppproperties :: Ω R. The parameters a, ρ and ccR material of the parameters a, ρ and depend on material properties of the pipe 2. MODELING 2. MODELING parameters a, ρ and ccangle depend on material properties of the pipe parameters a, ρ and depend on material properties of the pipe and the gas. α is the of the pipe and g is the gravitational and the gas. α is the angle of the and g is the gravitational and the α the of the pipe pipe and the It is common to define aa connected and directed simple graph and the gas. gas. α is isthis, the angle angle of the pipe and gg is is the gravitational gravitational It is common to define connected and directed simple graph constant. With we get following partial differential constant. With this, we get the following partial differential It is common to define a connected and directed simple graph It is common to define a connected and directed simpleaagraph constant. With this, we get the following partial differential G = G(V, E) representing the pipe network. This allows more constant. With this, we get the following partial differential G = G(V, E) representing the pipe network. This allows more equation, which describes the behavior in pipes equation, which describes the behavior in pipes G = G(V, E) representing the pipe network. This allows a more G = G(V,representation E) representing the pipe network. This allows aV more equation, which describes the behavior in pipes compact of the model equations. The set are equation, which describes the behavior in pipes compact representation of the model equations. The set V are 2 2 m compact representation of the theand model equations. The set V are are ∂∂∂ ppp (x,t) + aaa22 ∂∂∂ m compact representation of model equations. The set V the nodes and E are the edges we will describe the different the nodes and E are the edges and we will describe the different (x,t) = 0 + = 0 ∂∂tp (x,t) aA ∂∂m m the nodes nodes and E Eelements are the the edges edges and we will will describe describe the the different different xx (x,t) the and are and we node and edge in the following. (x,t) + (x,t) = 0 ∂t A ∂ node and edge elements in the following. (x,t) + (x,t) = 0 (1) (1) ∂t A ∂ x node and edge elements in the following. ∂ pp m ∂t A ∂x node and edge elements in the following. (1) ∂ ∂∂∂ m (1) (x,t) + A (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = ∂∂ pp (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = 0 m (x,t) + A 0 ∂ m ∂t ∂ x 2.1 Node Elements (x,t) + A (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = 0 ∂t 2.1 Node Elements (x,t) + A ∂∂ xx (x,t) + ρAg sin α + c|m(x,t)|m(x,t) = 0 ∂t 2.1 Node Node Elements Elements ∂t ∂x 2.1 Hyperbolic PDE pipe: Reservoirs are water sources with unlimited capacity. Thus, we Hyperbolic PDE pipe: Reservoirs are water sources with unlimited capacity. Thus, we Hyperbolic PDE pipe: Reservoirs are water sources with unlimited capacity. Thus, we Hyperbolic PDE pipe: Reservoirs are water sources with pressure unlimitedppcapacity. Thus, we Furthermore no assume that they have aa constant ss .. Furthermore assume that they have constant pressure no Further edge edge elements elements are are valves valves and and pumps, pumps, which which we we also also assume that that they holds have aaatconstant constant pressure pss ..arbitrary Furthermore no Further Furthermore no assume they have pressure p balance equation a reservoir, since an amount balance equation holds at aa reservoir, since an arbitrary amount Further edge elements are Further edgeextended elementsabstract. are valves valves and and pumps, pumps, which which we we also also omit in this balance equation holds at reservoir, since an arbitrary amount omit in this extended abstract. balance equation holds at a reservoir, since an arbitrary amount of water may leave or enter the reservoir. of water may leave or enter the reservoir. omit in this extended abstract. omit in this extended abstract. of water may leave or enter the reservoir. of water may leave or enter the reservoir. 2.3 Network Model 2.3 Network Model p = p 2.3 Network Model p = p ss 2.3 Network Model p = p p = pss From now on we aa network nn p many pipes and From now on we consider consider network with with pipes and now on consider with nn ppp many many pipes and From nowdemand on we we nodes consider network with many pipes and ndd many In contrast contrast to to reservoirs, reservoirs, tanks tanks have have limited limited capacity. capacity. Never Never the the From many demand nodes andaannnetwork many reservoirs. For each pipe rs n In and many reservoirs. For each pipe rs ni,ddwe In contrast to reservoirs, reservoirs, tanks have have limited capacity. Never the many demand nodes and n many reservoirs. For each pipe rs n In contrast to tanks limited capacity. Never the many demand nodes and n many reservoirs. For each pipe less, pressure can inor decrease even though the tank is full or get a flow m : Ω × I −→ R, Ω = [x , x ] ⊂ R, I = [t , T ]] ii : Ωii × I −→rsR, Ωii = [xL R 0 less, pressure can inor decrease even though the tank is full or i, we get a flow m , x ] ⊂ R, I = [t i i Li , xRi ] ⊂ R, I = [t0 ,, T less, pressure can inor decrease even though the tank is full or i, we get aa flow m : Ω × I −→ R, Ω = [x T i i i L R 0 less, pressure can inor decrease even though the tank is full or i, we get flow m : Ω × I −→ R, Ω = [x , x ] ⊂ R, I = [t T ]]ii i i empty respectively. We will not talk about tanks in more details and a pressure p : Ω ×I −→ R both depending on space x ∈ i −→ R bothi depending Li Ri on space x ∈ 0,Ω empty respectively. We will not talk about tanks in more details and aa pressure ppii ::i Ω Ω ii ×I n empty respectively. We will not talk about tanks in more details and pressure Ω ×I −→ R both depending on space x ∈ Ω d empty respectively. We will not talk about tanks in more details and atime pressure pii :The Ωii ×I −→variables R both depending space Ωii here. tt ∈ node are R here. and time ∈ I. I. The node variables are pppdd ::: IIIon−→ −→ Rnnnxdd ∈and and here. and here. and time time tt ∈ ∈ I. I. The The node node variables variables are are pdd : I −→ −→ R R d and and
2405-8963 © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Copyright © 2015, IFAC 260 Copyright 2015,responsibility IFAC 260Control. Peer review© under of International Federation of Automatic Copyright © 2015, IFAC 260 Copyright © 2015, IFAC 260 10.1016/j.ifacol.2015.05.201
MATHMOD 2015 February 18 - 20, 2015. Vienna, Austria
Sara Grundel et al. / IFAC-PapersOnLine 48-1 (2015) 260–261
prs : I −→ Rnrs , the pressures at demand nodes and reservoirs. We define mL (t) = (mi (xLi ,t))i∈A pi ,mR (t) = (mi (xRi ,t))i∈A pi . pL (t) = (pi (xLi ,t))i∈A pi ,pR (t) = (pi (xRi ,t))i∈A pi . mL being the vector with all pipe flows at their tail-node mR the flow vector at their head-nodes and similarly for pL pR . Note, that the node pressures coincide with the headtail pressures of the pipes. We call the vector of demand reservoir pressures by p.
and and and and
Last we define the following incidence matrices 1 if reservoir node j is head of pipe i n p ×nrs rs ∈ R (A ) = Ars R R ij 0 else −1 if reservoir node j is tail of pipe i n p ×nrs rs Ars (AL )i j = L ∈R 0 else 1 if demand node j is head of pipe i AR ∈ Rn p ×nd ,(AR )i j = 0 else −1 if demand node j is tail of flow i AL ∈ Rn p ×nd ,(AL )i j = 0 else We can combine them all and get the full incidence matrix A
n p ×+drs +nd rs A := (Ars R + AL AR + AL ) ∈ R With the help of these matrices we can write the spatial discretized system of equations as (2) p˙R + Dα (mL − mR ) = 0
m˙ L + Dβ AT p + γ + G(mL )mL = 0 (3) (4) AR mR + AL mL = qset (5) AR pd = pR AL pd = pL (6) p = p (7) Ars R R rs rs AL prs = pL (8) (9) prs = pset To obtain these equation it is crucial to chose a suitable spatial discretization. In particular the time derivative of the pressure is evaluated at the right end of the pipe and the time derivative of the flux at the left end. The size of the first two equation is the number of pipes, equation (4) the number of junction, (5,6,7,8) number of pipes and the last equation number of tanks. Dα is a diagonal matrix containing αi = a2i /Ai /�i on the diagonal and Dβ similar with βi = Ai /�i . γ is a vector with γi = ρi Ai g sin αi and G(mL ) is a diagonal matrix function such that G(mL )i = ci miL . The matrix A is the incidence matrix as described above. The vector qset has an entry for every demand node showing the given demand at that particular node given by qs (t) and similar is pset the vector of the given pressures ps at the reservoirs. In the modeling of the graph it is crucial to pick the direction of the edges such that every demand node has a right end of a pipe. This is possible for any topology as long as one of the nodes in the graph is not a demand node which means in our case it has to be a reservoir. We furthermore want all pipes the end in a reservoir to end in a left node there. 3. DECOUPLING By selecting a matrix Aselect which picks one pipe for each node that has a right end in that given node we can rewrite the system of equation in the variables pd and mL .
261
261
p˙d + Aselect Dα (−Cqset +CAmL ) = 0 pd T + γ + G(mL )mL = 0 m˙ L + Dβ A pset
We will explain in detail how we create the matrix C which is the crucial part in this decoupling process. This resulting ODE is of size nd + n p and has the general structure x˙ = T x + g(x,t) = f (x,t),
(10)
where the matrix T is given by 0 Aselect Dal phaCA T= , Dβ (AR + AL )T 0 and the vector x is combined by pd and mL . 4. IMEX In order to solve this stiff and nonlinear ODE we make use of implicit-explicit (IMEX) integration methods. This allow us to deal with the stiffness in an efficient way while not having to solve large-scale nonlinear problems. First order methods are of the flavor xn+1 − xn (11) = (1 − γ)T xn + γT xn+1 + g(xn ,t) h for γ ∈ [0, 1] and time setp h, which leads to the iteration
xn+1 = (1 − hγT )−1 (xn + (1 − γ)T xn + hg(xn ). We study convergence properties of that by analyzing the matrix T and the function g as well as the differences for several values of γ. If γ = 0 we get explicit Euler and if γ = 1 we get a combination of implicit Euler for the linear part an explicit Euler for the nonlinear part. We will also show the differences within this class of methods as well as the difference to second order methods following Ascher et al. (1995). 5. MODEL ORDER REDUCTION On the resulting ODE (10) we use the Model Order Reduction techniques Proper Orthogonal Decomposition (POD) together with Discrete Empirical Interpolation (DEIM) by Chaturantabut and Sorensen (2010). POD is a projection-based method where we find a projection matrix W such that the solution of (10) x ≈ W xˆ for xˆ in a lower dimensional space. The resulting low-dimensional ODE is then given by ˆ xˆ˙ = W T TW xˆ +W T g(W x,t). DEIM is then used to create a truely low-dimensional function approximating W T g(W x,t). ˆ 6. CONCLUSIONS The combination of POD-DEIM with the IMEX integration results in a significant speedup of simulation time. REFERENCES Ascher, U.M., Ruuth, S.J., and Wetton, B.T. (1995). Implicitexplicit methods for time-dependent partial differential equations. SIAM Journal on Numerical Analysis, 32(3), 797–823. Chaturantabut, S. and Sorensen, D.C. (2010). Nonlinear model reduction via discrete empirical interpolation. SIAM J. Sci. Comput., 32(5), 2737–2764. doi:10.1137/090766498.