An Optimization-Based State Estimation Framework for Large-Scale

Subscriber access provided by Kaohsiung Medical University

Process Systems Engineering

An Optimization-Based State Estimation Framework for Large-Scale Natural Gas Networks Jordan Jalving, and Victor M. Zavala Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b04124 • Publication Date (Web): 13 Mar 2018 Downloaded from http://pubs.acs.org on March 13, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

An Optimization-Based State Estimation Framework for Large-Scale Natural Gas Networks Jordan Jalving and Victor M. Zavala∗ Department of Chemical and Biological Engineering University of Wisconsin-Madison, 1415 Engineering Dr, Madison, WI 53706, USA E-mail: [email protected] Abstract We propose an optimization-based state estimation framework to track internal spacetime flow and pressure profiles of gas networks during dynamic transients. We find that the estimation problem is ill-posed (due to the infinite-dimensional nature of the states) and that this leads to instability of the estimator when short estimation horizons are used. To circumvent this issue, we propose moving horizon strategies that incorporate prior information. In particular, we propose a strategy that initializes the prior using steady-state information and compare its performance against a strategy that does not initialize the prior. We find that both strategies are capable of tracking the state profiles but that superior performance is obtained with steady-state prior initialization. We also find that, under the proposed framework, pressure sensor information at junctions is sufficient to track the state profiles. We also derive approximate transport models and show that some of these can be used to achieve significant computational speed-ups without sacrificing estimation performance. We show that the estimator can be easily implemented in the graph-based modeling framework Plasmo.jl and use a multi-pipeline network study to demonstrate the developments.

ACS Paragon Plus 1 Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Keywords: large-scale; optimization; state estimation; natural gas; networks; prior

1

Motivation and Background

A tighter interdependence between natural gas and electrical power infrastructures is motivating the development of new dynamic modeling, optimization, and control technologies. 1 New fracking technologies have led to abundant and cheaper natural gas supplies and has promoted significant investment in gas-fired power plants. Currently, gas-fired plants account for the majority of new generating capacity. 2 Gas-fired plants have smaller capacities than nuclear and coal counterparts but have much faster response times. From a power grid perspective, such dynamic flexibility is critical to manage intermittent wind and solar power sources and contingencies. On the other hand, frequent ramping of gas-fired power plants triggers complex spatio-temporal transient phenomena in gas transmission systems. These complex interactions present new risks and challenges in the operation of gas infrastructures. In particular, reliable gas delivery relies on the availability of line-pack, which is gas stored inside the network pipeline segments. 3 To build up and use line-pack, pipeline operators need to carefully synchronize compressors, gas injections, and withdrawals to ensure that sufficient line-pack is available at different locations and that pressure conditions are met at all delivery points. Recent literature has been concerned with modeling of gas transport dynamics with the ultimate goal of developing sophisticated optimal control and line-pack management strategies. 4–6 Such schemes rely on accurate estimates of the state of the gas network (pressure and flow pipeline profiles). State information is also essential from a power grid perspective because it provides a mechanism to determine if sufficient line-pack is available for gas-fired power plants. Internal pressure and flow profiles can be estimated reliably during steadystate conditions (e.g., flow is constant across the pipeline and pressure can be estimated using simple pressure drop calculations). Tracking internal pipeline profiles during dynamic transients, however, is significantly more challenging as complex spatiotemporal profiles are triggered by gas injections and withdrawals. To enable this, it is necessary to use more sophisticated state estimation techniques. Surprisingly, the literature on state estimation for


Page 2 of 33

Page 3 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


natural gas networks is quite sparse. Most existing studies have been concerned with obtaining state estimates to perform gas leak detection. Among these studies, estimation methods used include adaptive particle filters, 7 transfer functions, 8,9 and Kalman filters. 10–12 The drawbacks of these approaches when contrasted with optimization-based approaches such as moving-horizon estimation is well-documented and understood, 13,14 particularly in the context of large-scale systems. 15 To the best of our knowledge, optimization-based estimation approaches for natural gas networks have not been studied before. We attribute this to the computational complexity of the associated models, which comprise large sets of hyperbolic and nonlinear partial differential equations (PDEs). In this work, we study the problem of state estimation in natural gas networks using an optimization-based approach. Our goal is to estimate internal pipeline profiles during dynamic transients using flow and pressure sensor information at the junction nodes. The proposed framework incorporates detailed physical models (which are highly nonlinear PDEs) and complex constraint sets. We use a spatiotemporal discretization scheme to cast the estimation problem as a large-scale and sparse nonlinear programming problem (NLP). We find that the estimation problem is inherently ill-posed (due to the infinite-dimensional nature of the states) and that this leads to instability of the estimator when short horizons are used. To circumvent this issue, we propose moving horizon strategies and propose a strategy to obtain initial prior information using an easily-computable steady-state. We find that the estimator is capable of closely tracking the state during dynamic transients. We also find that, under the proposed framework, junction pressure sensor information is sufficient to track the states. We also derive different approximate transport models and show that a Weymouth approximation can achieve significant computational speed-ups without sacrificing much estimation performance. We use the graph-based modeling framework Plasmo.jl to implement the estimator and use a case study for a multi-pipeline network to demonstrate our developments. The paper is structured as follows. Section 2 provides an overview of the transport equations describing spatiotemporal dynamics in gas networks. In this section, we also provide model approximations that are used to improve computational tractability. In Section 3, we develop the estimation formulation and in Section 4 we develop the computational strategy



to solve this formulation. In Section 5 we provide the case study.

2

Gas Network Modeling

In this section, we present the governing equations that describe mass and momentum conservation in gas pipelines, equations for network interconnections (junctions), and expressions to compute compressor power and line-pack. Model derivations as well as nomenclature and units for all variables and parameters can be found in the supporting information.

2.1 Network Topology and Sets Gas networks are typically modeled as sets of pipeline segments (links) that are coupled together via junctions (nodes). We consider a network with a set L of links and a set N of nodes. Junctions may also include a gas supply or demand for which we define the sets S ⊆ N and D ⊆ N , respectively. For each node n ∈ N , we define its set of inlet and outlet out links Lin n := {`|rec(`) = n} and Ln := {`|snd(`) = n}. Here, rec(`) is the receiving node

of link ` and snd(`) is the sending node of link `. We define the subsets La ⊆ L and Lp ⊆ L to represent active and passive links, respectively. Active links have compressor stations at their inlet while passive ones do not. Figure 1 illustrates our notation using a simple segment that comprises a supply (at node n1 ), a demand (at node n3 ), a passive link (`1 ), and and active link (`2 ).

Figure 1: Sketch of pipeline segment along with model nomenclature.


Page 4 of 33

Page 5 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2.2

Link Equations

We assume an isothermal flow through a horizontal pipeline segment ` ∈ L. The isothermal mass and momentum conservation equations are presented in 16 and take the form: ∂ρ` (t, x) ∂(ρ` (t, x)v` (t, x)) + =0 ∂t ∂x ∂(ρ` (t, x)v` (t, x)) ∂(ρ` (t, x)v` (t, x)2 + p` (t, x)) λ` + =− ρ` (t, x)v` (t, x) |v` (t, x)| ∂t ∂x 2D`

(2.1a) (2.1b)

These equations are also known as the Euler equations. Here, the link diameter is D` and the friction coefficient is λ` . The states of each link vary in space and time and are given by the gas density ρ` (t, x), gas velocity v` (t, x), and gas pressure p` (t, x). From these fundamental quantities we can derive expressions for transversal area A` , volumetric flow q` (t, x), and mass flow rate f` (t, x): 1 A` = πD`2 4

(2.2a)

q` (t, x) = v` (t, x)A`

(2.2b)

f` (t, x) = ρ` (t, x)v` (t, x)A` .

(2.2c)

For an ideal gas, the density and pressure are related by the gas speed of sound c (2.3): p` (t, x) = c2 ρ` (t, x)

(2.3)

Typically, gas pipeline operators monitor pressure and mass flow rate since they can be directly measured. To represent the transport equations in terms of these states, we use relations (2.2) and (2.3) to cast (2.1) into (2.4). The details of the reformulation can be found in the supporting information. The final form of the transport equations is: c2 ∂f` (t, x) ∂p` (t, x) + = 0, ` ∈ L (2.4a) ∂t A` ∂x ∂f` (t, x) 2c2 f` (t, x) ∂f` (t, x) + ∂t A` p` (t, x) ∂x c2 f` (t, x)2 ∂p` (t, x) ∂p` (t, x) 8c2 λ` A` f` (t, x) − + A = − |f` (t, x)| , ` ∈ L (2.4b) ` A` p` (t, x)2 ∂x ∂x π 2 D`5 p` (t, x)



Page 6 of 33

We can improve the numerical behavior of formulation (2.4) by scaling pressure and flow variables. We define constants c1,` , c2,` and c3,` (see supporting information) and cast the Euler equations in the scaled form: ∂p` (t, x) ∂f` (t, x) = −c1,` , `∈L ∂t ∂x ∂f` (t, x) f` (t, x) ∂f` (t, x) f` (t, x)2 ∂p` (t, x) = −2 · c1,` + c1,` ∂t p` (t, x) ∂x p` (t, x)2 ∂x ∂p` (t, x) f` (t, x) − c2,` − c3,` |f` (t, x)| , ` ∈ L. ∂x p` (t, x)

2.3

(2.5a)

(2.5b)

Compressor Equations

We define the pressures at nodes as θn (t), n ∈ N and define the boost pressures for the active links ` ∈ La as ∆θ` (t), which is the additional pressure added to the link inlet by the compressor. We can thus define the boundary conditions for the link pressures as:

p` (L` , t) = θrec(`) (t), ` ∈ L

(2.6a)

p` (0, t) = θsnd(`) (t), ` ∈ Lp

(2.6b)

p` (0, t) = θsnd(`) (t) + ∆θ` (t), ` ∈ La .

(2.6c)

The total power consumed by each compressor as part of each active link is then given by: 

θsnd(`) (t) + ∆θ` (t) P` (t) = cp · T · f` (0, t)  θsnd(`) (t)

γ−1 γ

 − 1 , ` ∈ La

(2.7)

where cp , T , and γ are constant parameters defined in the supporting information.

2.4

Node Equations

For convenience, we introduce extra variables for each link fìn (t) and fòut (t), that represent the inlet and outlet flow rates for each link. With this, we can express the node balances as: X `∈Lrec n

fòut (t) −

X `∈Lsnd n

fìn (t) +

X i∈Sn

gi (t) −

X j∈Dn


dj (t) = 0, n ∈ N

(2.8)

Page 7 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


We note that the supply flows si (t), i ∈ S and the demand flows dj (t), j ∈ D are algebraic states of the system (if the pressures of the corresponding nodes are specified). 5 The flow variables can also be used to define boundary conditions for the link flows as:

f` (t, L` ) = fòut (t)

(2.9a)

f` (t, 0) = fìn (t).

(2.9b)

In this formulation we assume the direction of the flow is constant and given. Strategies for modeling flow reversals are proposed in 17 and more recently. 18

2.5 Line-Pack Line-pack refers to the actual inventory of gas stored within a pipeline segment at any time. This quantity is important from an operational stand-point, since it provides flexibility to account for unforeseen situations. The amount of linepack in any pipeline segment m` can be calculated as: L`

Z m` (t) = A`

ρ` (t, x)dx.

(2.10)

0

By using the pressure law (2.3) and the scaling constant c1,` , we can express the line-pack in each link in terms of the gas pressure as:

m` (t) =

1 c1,`

Z

L`

p` (t, x)dx.

(2.11)

0

2.6 Approximate Models The transport equations are highly nonlinear and computationally challenging to solve, particularly when large networks and long time horizons need to be considered. To deal with this issue, we consider different approximate models proposed in literature. 19

2.6.1 Weymouth Approximation The Weymouth approximation 16 is derived by noticing that the momentum term ∂x (ρv 2 ) in (2.1b) has a negligible effect compared to the rest of the terms. By dropping this term from ACS Paragon Plus 7 Environment


Page 8 of 33

the Euler equations we obtain: ∂p` (t, x) ∂f` (t, x) + c1,` =0 ∂t ∂x ∂f` (t, x) ∂p` (t, x) f` (t, x) + c2,` + c3,` |f` (t, x)| = 0. ∂t ∂x p` (t, x)

(2.12a) (2.12b)

This is a coupled set of wave equations for flow and pressure with a nonlinear source term. This system is easier to solve than the Euler equations but it is still a challenging hyperbolic equation.

2.6.2 Adiabatic Approximation The adiabatic approximation 20 can also be derived from the isothermal Euler equations. We obtain this form by neglecting the momentum ∂x (ρv 2 ) term and the partial derivative ∂t (ρv). This gives the following formulation: ∂p` (t, x) ∂f` (t, x) + c1,` =0 ∂t ∂x f` (t, x) ∂p` (t, x) + c3,` |f` (t, x)| = 0. c2,` ∂x p` (t, x)

(2.13a) (2.13b)

2.6.3 Steady-State Model The steady-state version of the Euler conditions is given by 16 : ∂f` (t, x) =0 ∂x

(2.14a)

∂p` (t, x) f (t, x) + c3,` |f` (t, x)| = 0 ∂x p` (t, x)

(2.14b)

c1,` c2,`

This form is known as the quasi-static approximation and is decoupled in time. This form of the Euler equations is mostly used as an approximation to capture long-term (planning) behavior and can be used as a proxy to capture situations under which the gas network is under “calm” conditions.


Page 9 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3 State Estimation Framework In this section, we derive a state estimation framework to infer pressure and flow profiles in real-time using junction sensor data.

3.1

Space-Time Discretization

To enable computational implementation, we discretize the Euler equations by using an implicit Euler scheme in time and a finite difference scheme in space. We define the discrete set of times as T and the discrete set of spatial points X . The discrete version of the Euler equations (2.5) is given by: p`,t+1,k − p`,t,k c2 f`,t+1,k+1 − f`,t+1,k =− , t ∈ T ,k ∈ X ∆t A` ∆x` f`,t+1,k − f`,t,k 2c2 f`,t+1,k f`,t+1,k+1 − f`,t+1,k =− ∆t A` p`,t+1,k ∆x` +

(3.15a)

2 c2 f`,t+1,k p`,t+1,k+1 − p`,t+1,k A` p2`,t+1,k ∆x`

− A`

p`,t+1,k+1 − p`,t+1,k 8c2 λ` A` f`,t+1,k |f`,t+1,k | − 2 5 , ∆x` p`,t+1,k π D`

t ∈ T ,k ∈ X. (3.15b)

We use the compact notation fL,T ,X to denote the entire set of flows for all links L, times T , and spatial points X . We use a similar notation for the rest of the variables such as the in and pressure profiles pL,T ,X , junction pressures θN ,T , as well as input/output flows fL,T out . fL,T

3.2 Sensor Measurements Typical sensor measurements available in gas networks are the pressures at nodal junctions and flows directed in and out of the junctions. 21 In this work, we assume that only pressure measurements are available. We will demonstrate that such information is sufficient to infer the internal state of the pipelines. For simplicity of implementation, we also assume sensor measurement noise is sufficiently filtered. An example of applying output filtering to an optimization-based estimator can be found in. 15 We consider sensor measurements disACS Paragon Plus 9 Environment


tributed over a historical (finite) time horizon consisting of N sampling times of equal length δ = tj − tj−1 . At sampling time tj , we define the past time window Tj = {tj−N , tj−N +1 , ..., tj }. For convenience, we also define the set T¯j =Tj \ {tj }. Using a history of measurements over Tj we seek to estimate the current state at time tj . Specifically, at sampling time tj , we estim . In addition, we assume that mate the current state using the pressure sensor history θN ,Tj

measurements for the control (input) variables are available; these include the compressor m boosts ∆θL (this can be inferred from suction and discharge pressures). We use ITmj to a ,Tj

denote all the sensor information over time window Tj . With this information we estimate e the dynamic states (pe`,Tj ,X , f`,T ) as well as the algebraic states (seS,Tj and deD,Tj ). We use j ,X e ITej = (pe`,tj ,X , f`,t , seS,Tj , deD,Tj ) to denote the estimated state trajectory. j ,X


Page 10 of 33

Page 11 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.3

Estimation Problem

The states are estimated at sampling time tj by solving the following estimation problem: XX

min

XX m 2 (p`,tj−N ,k − p`,k )2 + (f`,tj−N ,k − f `,k )2 + (θn,t − θn,t )

s.t.

(3.16a)

t∈Tj n∈N

`∈L k∈X

p`,t+1,k − p`,t,k f`,t+1,k+1 − f`,t+1,k = −c1,` , ` ∈ L, t ∈ T j , k ∈ X ∆t ∆x` f`,t+1,k − f`,t,k 2c2 f`,t+1,k f`,t+1,k+1 − f`,t+1,k =− ∆t A` p`,t+1,k ∆x`

(3.16b) (3.16c)

2

+

c2 f`,t+1,k p`,t+1,k+1 − p`,t+1,k A` p2`,t+1,k ∆x`

p`,t+1,k+1 − p`,t+1,k 8c2 λ` A` f`,t+1,k |f`,t+1,k | − 2 5 , t ∈ T j, k ∈ X ∆x` p`,t+1,k π D` X X X X f`,L,t − f`,0,t + gi,t − dj,t = 0, t ∈ Tj , n ∈ N − A`

`∈Lrec n

`∈Lsnd n

i∈Sn

(3.16d) (3.16e) (3.16f)

j∈Dn

p`,t,L = θrec(`),t , ` ∈ L, t ∈ Tj

(3.16g)

p`,t,0 = θsnd(`),t , ` ∈ Lp , t ∈ Tj

(3.16h)

m , ` ∈ La , t ∈ Tj p`,t,0 = θsnd(`),t + ∆θ`,t

(3.16i)

out , ` ∈ L, t ∈ Tj f`,t,L = f`,t

(3.16j)

in f`,t,0 = f`,t , ` ∈ L, t ∈ Tj   γ−1 m γ θsnd(`),t + ∆θ`,t P`,t = cp · T · f`,t,0  − 1 , ` ∈ La , t ∈ Tj θsnd(`),t

m`,t =

1 X c1,`

p`,t,k ∆x` ,

` ∈ L, t ∈ Tj .

(3.16k) (3.16l) (3.16m)

k∈X

The first term in the objective (3.16a) is known as the prior term, which summarizes past information through the prior information p`,X and f `,X . We use I¯ = (p`,X , f `,X ) to denote this prior information. The second term in the objective are standard squared error terms, which penalize deviations of the model outputs from the corresponding sensor measurements. We use equal weighted objective terms to simplify the presentation, but note that different weights applied to the individual objective terms can be used to tune estimator performance. We also omit state noise objective terms (model-plant mismatch) which can be used to improve estimator performance and refer to 15 as an example of including this kind of term. The constraints of the estimator problem are the discretized Euler equations deACS Paragon Plus 11 Environment


Page 12 of 33

fined over the time domain Tj . A key benefit of optimization-based estimation approaches is that they can accommodate constraints directly. This provides a mechanism to add prior information. In our formulation we add physical bounds for all variables of the form:

θnL ≤ θn,t ≤ θnU , n ∈ N , t ∈ Tj

(3.17a)

P`L ≤ P`,t ≤ PÙ , ` ∈ La , t ∈ Tj

(3.17b)

f`L ≤ f`,t,k ≤ fÙ , ` ∈ L, k ∈ X , t ∈ Tj

(3.17c)

U pL ` ≤ p`,t,k ≤ p` , ` ∈ L, k ∈ X , t ∈ Tj .

(3.17d)

¯ I) ˆ to express the estimator model in compact form. We use the notation ITej ← M(ITmj , I, Here, Iˆ is an initial guess (warm-start) for the estimated states.

3.4

Prior Information

The initial axial profiles p`,tj−N ,X and f`,tj−N ,X in problem (3.16) are free variables that are estimated. As a result, the estimation problem has a large number of degrees of freedom, given by the number of link segments times the cardinality of the set X (the mesh resolution). Specifically, the number of degrees of freedom is 2 · |L| · |X |. The number of degrees of freedom might make the problem inherently ill-posed, in the sense that sensor information might be insufficient to infer the initial axial profiles uniquely. This problem becomes exacerbated when the length of the estimation horizon |Tj | is short. This presents a practical problem because short horizons are often desired from a computational performance standpoint (to make the problem more tractable). A mechanism to mitigate ill-posedness is to use regularization in the form of prior information. Prior information is also critical to ensure that the estimator converges to the actual states. 22,23 In other words, in the absence of prior information (or if a poor prior is used), the estimator is not guaranteed to converge. Obtaining initial prior information in complex physical models, however, can be a challenging task. In this work, we propose to obtain an initial prior by calculating a steady-state solution for the transport equations. Consider an initial horizon T0 with sensor information history ITm0 . We use the first measurement in the history Itm0−N to solve the following set of nonlinear ACS Paragon Plus 12 Environment

Page 13 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


equations: f`,k+1 − f`,k , ` ∈ L, k ∈ X ∆x` 2c2 f`,k f`,k+1 − f`,k 0=− A` p`,k ∆x`

0 = −c1,`

2 c2 f`,k p`,k+1 − p`,k + A` p2`,k ∆x`

p`,k+1 − p`,k 8c2 λ` A` f`,k |f`,k | − 2 5 , ` ∈ L, k ∈ X ∆x` p`,k π D` X X X X f`,0 + gi − dj = 0, n ∈ N f`,L − `∈Lsnd n

i∈Sn

(3.18b) (3.18c)

− A`

`∈Lrec n

(3.18a)

(3.18d) (3.18e)

j∈Dn

p`,L = θrec(`) , ` ∈ L

(3.18f)

p`,0 = θsnd(`) , ` ∈ Lp

(3.18g)

m , ` ∈ La p`,0 = θsnd(`) + ∆θ`,t 0−N

(3.18h)

f`,L = fòut , ` ∈ L

(3.18i)

f`,0 = fìn , ` ∈ L

(3.18j)

m θn = θn,t , n ∈ S × D. 0−N

(3.18k)

The solution of this problem gives us the prior information I p . The steady-state formulation is motivated by the degree of freedom analysis presented in, 5 which concludes that the steady-state of the gas network can be fully defined by fixing supply and demand pressures and boost pressures. The above formulation also highlights that the state of natural gas networks can be easily computed when the system is at rest. During dynamic transients, however, the states must be inferred using state estimation techniques. At time t0 (with j = 0), we use the steady-state prior to solve the estimator problem ¯ I). ˆ We use the solution of the estimation problem to update the prior as ITe0 ← M(ITm0 , I, I¯ ← Itej−N +1 . In other words, we use the model state prediction for time tj−N +1 to update e the prior as p¯`,X ← pe`,tj−N +1 ,X and f¯`,X ← f`,t . The prior information is updated in j−N +1 ,X

this form at all subsequent sampling times j > 0. In our analysis we contrast the performance of the estimator using the initial steady-state prior to that of an estimator in which the prior is not initialized due to lack of knowledge (the prior is updated at subsequent sampling times time using estimated states). We also ACS Paragon Plus 13 Environment


consider an estimator in which no prior information is used at all (the prior term in (3.16a) is removed at all sampling times). With this last strategy we seek to evaluate if sensor measurements are sufficient to track the state profiles.

3.5

Summary of Estimation Algorithm

The state estimator is implemented over a receding (moving) horizon of fixed length. We recursively shift the sensor data window and adjust the prior accordingly using the state model prediction. This provides a mechanism to avoid accumulating sensor data over an infinitely long horizon. The estimator computational workflow is presented in Algorithm 1.

Algorithm 1 State Estimation ¯ and warm-start I. ˆ Given N , prior information I, Set sampling time j ← 0 and time window Tj ← {tj−N : tj }. loop Obtain sensor information ITmj from system. ¯ I). ˆ Solve estimator problem ITej ← M(ITmj , I, Update warm-start information Iˆ ← ITej . Update sampling time j ← j + 1 and time window Tj ← {tj−N : tj }. Update prior information I¯ ← Itej−N . end loop

4 Computational Implementation In this section we discuss computational aspects associated to the implementation of the estimator.

4.1

Graph-Based Modeling

Developing gas network models is a challenging task because of the need to incorporate complex sets of equations. Moreover, state estimation applications require the repetitive solution of optimization problems, which creates complex computational workflows where large amounts of input data and solution data is dynamically updated. To facilitate the implementation of state estimation applications in gas networks, we use Plasmo.jl. 24


Page 14 of 33

Page 15 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Plasmo.jl is a graph-based modeling package written in Julia 25 and built upon the algebraic modeling language JuMP. 26 Plasmo.jl facilitates the expression of complex models by modularization, as done in other modeling languages such as Pyomo. 27 In particular, Plasmo.jl enables the creation of component models that are independent of the graph structure of the problem (e.g., the gas network topology). For instance, pipeline segments can be represented as individual edges with associated component models containing the transport equations. Similarly, junctions are represented by individual nodes with associated component models containing the junction conservation equations (see figure 2). Such component models can be implemented without having any knowledge of the network topology. Once the network topology is defined, the component models are attached to the nodes and edges. This modeling paradigm facilitates the model re-use (one can create models for different network topologies by re-using component models). Moreover, given a fixed network topology, one can create models of different resolution by simply swapping component models (e.g., Euler equations can be exchanged by the Weymouth or quasi-static approximation). The modular modeling abstraction provided by Plasmo.jl also facilitates data management because the input and output data of the component models remains modularized. Modularization also facilitates warm-starting because it is possible to retain the component model solutions and the underlying algebraic structure and use it in subsequent problem instances.

Figure 2: Modeling in Plasmo.jl using component models for nodes and edges.



Page 16 of 33

By implementing the state estimator in the programming language Julia, one can use high-level data and function constructs to facilitate the implementation of complex applications, such as moving horizon estimation. In Figure 3 we present the high-level front end needed to implement the estimator.

4.2 Optimization Solver The estimator problem is a large-scale and sparse nonlinear programming problem (NLP) that can be expressed in the abstract form:

min f (w)

(4.19)

s.t. c(w) = 0

(4.20)

w ≥ 0.

(4.21)

We use the interior-point solver Ipopt to solve this problem. 28 This strategy allows us to handle highly nonlinear models as well as large numbers of constraints and degrees of freedom. Moreover, the linear algebra operations performed in Ipopt can help us identify if sensor information is sufficient to infer the states reliably. 15 This is based on the observation that a minimizer of (4.19) is unique when the Karush-Kush-Tucker (KKT) matrix has a correct inertia.

5

Case Study

In this section, we provide computational results to demonstrate the performance of the proposed state estimation framework. Animated results depicting estimator performance can be found at https://github.com/zavalab/JuliaBox/tree/master/MHEGasNetwork.

5.1 Multi-pipeline system We consider a system of 13 connected pipelines in series, as depicted in Figure 4. The network includes a gas supply at one end, a time-varying demand at the other end, and 12 compressor stations. We run a simulation of this network with typical operating conditions ACS Paragon Plus 16 Environment

Page 17 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3: Implementation of moving horizon estimator in Plasmo.jl 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

using Plasmo, Ipopt, JuMP include("data_set.jl") # defines measurement info include("model_functions.jl") # defines ”solve ss problem” and ”create estimation problem” time_grid = Grid(0,T,Nt) space_grid = Grid(0,L,Nx) #Create Graph Topology graph = PlasmoGraph() for l in Links node = add_node(graph) end for n in Nodes from_node = get_node(graph,n) to_node = get_node(graph,n+1) edge = add_edge(graph,from_node,to_node) end for j = 1:T # loop over sampling times # set prior information if j == 1 #use steady-state p_prior,f_prior = solve_ss_problem(graph,measurement_info[j]) else #update using state estimate k = horizon_points[1] p_prior = px_estimate[k] f_prior = fx_estimate[k] end # create estimation problem status=create_estimation_problem(graph,measurement_info[j]) # update warm-start information if i > 1 status=set_solution(state_estimate,graph) end # solve estimation problem status=solve(graph) state_estimate = get_solution(graph) end

using the Euler equations 2.1. This simulation is for a 24 hour horizon with a time discretization of 15-minute intervals. The control policy used for this simulation is obtained by solving an optimal control problem that tracks a sudden demand withdrawal. 5 Figure 5 summarizes the simulation results. The network is at steady-state at t1 (given by a constant flow profile) and experiences a sudden gas withdrawal at the demand node at time t10 that induces a complex dynamic transient that triggers a forward propagating wave. At time t40 , the demand returns to its original state, which results in a backward ACS Paragon Plus 17 Environment


∆θℓ2

pℓ2 (Lℓ2 ) = θn3

∆θℓ3

∆θℓ11

n2

n1 ℓ1

n3

n11

s

n13

n12

ℓ2 θ n2

d

fℓ1 (0) = fℓin 11

θ n1

ℓ11

pℓ2 (0) = θn2 + ∆θℓ2

fℓ1 (Lℓ1 ) =

ℓ12

θn13

fℓout 11

Figure 4: Sketch of multi-pipeline system. propagating wave. Towards the end of the time horizon, the system is pressurized to return the total line-pack to its original value. This produces the final dynamic flow profile at t96

Flow [scmhr* 104 ]

Demand [scmhr* 104 ]

(which is not a steady-state).

Pressure [bar]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 33

46 44 42 40 50 48 46 44 42 40

t10

t1

0

t40

5

10

tfinal

15

Time [hr]

20

tfinal

t40 t10

t1

0

50

100 Distance [km]

150

0

50

100 Distance [km]

150

60 50 40

Figure 5: Simulated transient. Demand step (top), flow profiles (middle), and pressure profiles (bottom).

5.2

Implementation Details

For our implementation we use the Weymouth formulation (2.12) and a space-time discretization mesh with Nx = 50 and Nt = 96. We study the performance of the estimator both with and without prior information for horizon lengths of N =2,5, and 10. Each opACS Paragon Plus 18 Environment

Page 19 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


timization problem is solved with Ipopt using the MA57 linear solver of Harwell. We use default algorithmic parameters and set the maximum number of iterations to 3,000. All results were obtained on an Intel(R) Core(TM) i7-6600U 2.60 GHz CPU. We assume that in the nominal case, the available system measurements are the pressures at each node and that metering stations are not necessarily installed at junctions. We therefore discard the flow measurements from the problem formulation but note that they are straight-forward to include if available.

5.3

Estimator Performance

We examine the performance for the state estimator. We analyze convergence under different prior update strategies and the effect of using model approximations.

5.3.1 Convergence Figure 6 shows the relative errors for both pressure and flow at each time step in the algorithm. We define scaled errors for both pressure and flow estimates at each sampling time tj as:

p,tj f,tj

X X p`,k,tj − pe`,k,tj = p`,k,tj `∈L k∈X e X X f`,k,tj − f`,k,t j = . f`,k,tj

(5.22a) (5.22b)

`∈L k∈X

Here, p`,k,tj and f`,k,tj are the true states obtained from the simulation at current sampling time tj . For the estimator that does not use prior information we can see that, for N = 2, there is a considerable error in the pressure and flow estimates at each sample time and the estimates do not converge. Moreover, about half of the estimator problems failed to solve with Ipopt over the complete simulation horizon, resulting in the sporadic error profile. By increasing the horizon length, the estimation error improves significantly, but the estimator is still prone to convergence issues. For the case N = 5, the estimation problem failed to solve at t10 . This behavior is a clear manifestation of the ill-posedness of the problem. No failures



were observed for a long horizon of N = 10. 102

N=2 N=5 N = 10

101

100

100

N=2 N=5 N = 10

f, tj

101

10

1

10

2

10

3

10

4

Estimate Error =

p, tj

102

Estimate Error =

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 33

25

50

Time Step

75

10

1

10

2

10

3

10

4

25

50

Time Step

75

Figure 6: Time evolution of estimator error for pressure (left) and flow (right) without prior information. By including prior information, the estimator error converges for every horizon length (even for a short horizon of N = 2). In Figure 7 we present the estimation errors with no prior information (I), with no initial prior information (II), and with an initial steady-state prior (III) for N = 2. We recall that in the case of (I), the prior term in (3.16a) is dropped at all sampling times while in (II) it is only dropped in the first sampling time (assumes that no information is initially available but this is updated at later sampling times using the state prediction). It is clear that the error levels of (I) are drastically improved by using prior information (II and III) and that the steady-state prior (III) provides the best performance. This highlights that, if the pipeline system is at rest at the beginning of the horizon, using an initial steadystate prior will be sufficient to track the state during the dynamic transient. This result has important practical implications because most real pipeline systems are indeed returned to rest conditions nearly every day. Remarkably, even if no initial prior is used, estimator (II) is also able to track the states and the errors decrease to the same levels of (III) for the flow profiles. The error level for the pressure profiles do not decrease to the same levels but remain fairly small. Unfortunately, under approach (II) we run the risk that the estimation ACS Paragon Plus 20 Environment

Page 21 of 33

problem at the first sampling time cannot be solved due to ill-posedness. Consequently, we conclude that the steady-state approach provides a more reliable way to initialize the prior. 102

102

102

p, tj f, tj

Estimate Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


p, tj f, tj

p, tj f, tj

101

101

101

100

100

100

10

1

10

1

10

1

10

2

10

2

10

2

10

3

10

3

10

3

10

4

10

4

10

4

25

50

Timestep

75

25

50

Time Step

75

25

50

Time Step

75

Figure 7: Estimation error with N = 2 for no prior information (I) (left), no initial prior information (II) (middle), and initial steady-state prior (III) info (right). Figure 8 compares the flow profile from the original simulation, the flow profiles estimated with (II), and the flow profiles estimated with (III) for N = 2. Notice the initial error profile of (II) corresponds with its estimate error in Figure 7 (obtained with I). While there is considerable error at t1 (the black line) the error converges by time t10 (the red line). In contrast, the flow profiles corresponding to the steady-state prior are visually indistinguishable with the true values. Figure 9 compares the no-prior (I) and steady-state prior (III) estimates with the true state values at sampling time t45 with N = 2 and N = 10. Here, it is obvious that (I) introduces high instability in the flow estimates when the horizon is short. Moreover, the estimation error persists up to time t45 , which highlights that prior information (as done in II and III) must be used to ensure convergence.



Flow [scmhr* 104 ]

50 48

tfinal

46

t10

42

t1 0

50

Distance [km]

Flow [scmhr* 104 ]

50 48

100

150

tfinal

t1

46

t40

44

t10

42 40

0

50

Distance [km]

50

Flow [scmhr* 104 ]

t40

44 40

100

48

150

tfinal

46

t40

44

t10

42

t1

40

0

50

Distance [km]

100

150

Figure 8: Comparison of simulated flow profile (top) and reconstructed state with N = 2 using no-initial-prior (II) (middle) and the steady-state prior (III) (bottom) 60

50

45

True N=2 N = 10

52.5

Flow [scmhr* 104 ]

Flow [scmhr* 104 ]

55.0

True N=2 N = 10

55

50.0 47.5 45.0 42.5

40

0

50

100

Distance [km]

150

40.0

True N=2 N = 10

50

40

0

50

100

Distance [km]

150

True N=2 N = 10

60

Pressure[bar]

60

Pressure[bar]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 33

50

40 0

50

100

Distance [km]

150

0

50

100

Distance [km]

150

Figure 9: Time snapshot of estimated and true profiles with N = 2 and N = 10 at sampling time t45 . No prior information (I) (left) and steady-state prior (III) (right). 5.3.2 Line-Pack Estimates The system line-pack is an important estimate with regard to system resilience, since it provides flexibility to absorb fluctuations in demand and equipment failures. Figure 10 shows the time evolution of the true and estimated line-pack in each pipeline segment. The results ACS Paragon Plus 22 Environment

Page 23 of 33

show that the estimator can accurately track the line-pack distribution in the system during dynamic transients. To give an idea on the value of this information, at a commercial price of $8.77 per thousand cubic feet (https://www.eia.gov/dnav/ng/ng_pri_sum_dcu_ nus_m.htm 29 ), the total value of natural gas stored in the first pipeline segment (300 km) is approximately $168, 000.

300

[kg] Line-pack 1000

[kg] Line-pack 1000

150

200

tj = 1 tj = 10 tj = 40 tj = 96

100

100 1 2 3 4 5 6 7 8 9 10 11 12 13 300 100 100 100 100 100 100 100 100 100 100 100 300 km km km km km km km km km km km km km

2 100 km

Pipeline Segment

3 100 km

4 100 km

5 100 km

Pipeline Segment

6 100 km

7 100 km

150 300

[kg] Line-pack 1000

[kg] Line-pack 1000

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


200

tj = 1 tj = 10 tj = 40 tj = 96

100

100 1 2 3 4 5 6 7 8 9 10 11 12 13 300 100 100 100 100 100 100 100 100 100 100 100 300 km km km km km km km km km km km km km

Pipeline Segment

2 100 km

3 100 km

4 100 km

5 100 km

Pipeline Segment

6 100 km

7 100 km

Figure 10: Time evolution of true line-pack axial profile (top left) with magnified region (top right). Time evolution of estimated line-pack dynamics (bottom left) with corresponding magnified region (bottom right) 5.3.3 Model Approximations We now analyze the computational performance of estimator (III). In typical applications, the state estimates are used to generate planning and control decisions and/or to quickly detect abnormal events (e.g., unstable flows inside the pipelines). Consequently, the estimator problem must be solved at a high frequency. Figure 11 shows computational times using a horizon length of N = 10. We note that the first problem takes over a minute to complete because the initial warm-state information is not as reliable. After this, the estimator problem solves in less than five seconds due to the availability of improved warm-start information. Figure 11 also shows the time-average solution times for different horizon lengths. For the longest horizon, the estimator problem solves in around 40 seconds. We also highlight that the use prior information induces computational stability. Without prior ACS Paragon Plus 23 Environment


information, the estimation problem required hundreds to thousands of iterations to solve, whereas with prior information the problem solved in less than 30 iterations.

CPU Time [sec]

60

40

20

0

0

20

40

40

Average CPU Time [sec]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 33

Time Step

60

80

30

20

10

0

N=2

N=5

N = 10

Horizon Length

N = 20

N = 30

Figure 11: Computational performance for N = 10 (top) and for increasing N (bottom). With both our accuracy and computational goals in mind, we also quantify the effect of using approximate models in the estimator formulation (III). Figure 12 compares the estimation error for the different approximations relative to the true state of the full-resolution Euler equations (2.1). In particular, we consider the Weymouth approximation (2.12) and the adiabatic approximation (2.13) with N = 10. As can be seen, the Weymouth approximation achieves reasonably accurate estimates and the adiabatic approximation shows significant deviations towards the initial and final sampling times. This is because the adiabatic approximation neglects the advection term ∂t (ρv) in the momentum balance, which is in fact significant during strong transients. We thus attribute the observed error to the strong compressor ramping used during the depletion and replenishment of the line-pack at the beginning and end of the planning horizon (see Figure 10). We compare computational times for the three different model approximations in Figure 13. The Weymouth approximation reduces the computational time by one half. The adiabatic approximation requires the least time but the time improvements over the Weymouth approximation are marginal. We can thus conclude that the Weymouth approximation provides a promising model reduction technique from a state estimation stand-point. ACS Paragon Plus 24 Environment

Page 25 of 33

100

100

100

Estimate Error

p, tj f, tj

p, tj f, tj

p, tj f, tj

10

1

10

1

10

1

10

2

10

2

10

2

10

3

10

3

10

3

10

4

10

4

10

4

10

5

10

5

10

5

25

50

Time Step

75

25

50

Time Step

75

25

50

Time Step

75

Figure 12: Error profiles for estimator with different model approximations for N = 10. Adiabatic (left), Weymouth (middle), and Euler (right).

Euler Weymouth Adiabatic

100

CPU Time [sec]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


50

0

t=1

t=5

t = 10

t = 40

Time Step

t = 60

t = 90

average

Figure 13: Computational performance for model approximations.



5.3.4 Highly Transient Systems We finally demonstrate the estimator performance using a highly dynamic demand scenario for the 13 pipeline system as depicted in figure 14. The estimator is additionally started during the transient phase of operation (at t20 in the simulation) to evaluate its performance in non-steady state conditions. Once again, the estimator error converges for every horizon length (even for the short horizon N = 2) by including the steady-state prior information. Figure 15 presents the estimation error with no prior information (I), with no initial prior information (II), and with an initial steady-state prior (III) for N = 2 for the highly dynamic case. As was the case in figure 6 the error for (I) is drastically improved using prior information from (II) or (III). In contrast to the previous case, formulation (II) shows better convergence than (III). This makes intuitive sense; formulation (II) is able to find a better initial arrival cost for the transient system versus using the steady-state prior. However, the approach using (II) still runs the risk of failing due to ill-posedness at the first sampling time. We therefore again conclude that the steady-state approach provides a more reliable way to initialize the prior. We further see from figure 16 that the flow estimate error could be considered acceptable for long operation (for N = 2).

6

Conclusions and Future Work

We have presented a state estimation framework to track internal space-time flow and pressure profiles of gas networks during dynamic transients. Our results indicate that prior information is key to infer the states using the limited sensor information and horizons. We demonstrate that prior information can be obtained by solving the steady-state Euler equations and that this information is sufficient to track the states during dynamic transients. We also analyze the performance of estimators that use different model approximations. We find that a Weymouth approximation can reduce solution times by half without sacrificing accuracy. We also demonstrate that graph-based modeling tools can greatly facilitate the implementation of complex state estimation applications. As part of future work, we are interested in solving larger scale problems by using decomposition techniques. Model ACS Paragon Plus 26 Environment

Page 26 of 33

Flow [scmhr* 104 ] Pressure [bar]

90 80 70 60 50 40 t1 0 90 80 70 60 50 40 30 55 50 45 40 35

t10

t20 (Estimator starts)

5

10

15

Time [hr]

tfinal

t70

20 t20 t70 t10 t1 tfinal

0

50

100 Distance [km]

150

0

50

100 Distance [km]

150

Figure 14: Highly dynamic simulated transient. Demand profile (top), flow profiles (middle), pressure profiles (bottom).

102

102

102

p, tj f, tj

Estimate Error

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Demand [scmhr* 104 ]

Page 27 of 33

p, tj f, tj

p, tj f, tj

101

101

101

100

100

100

10

1

10

1

10

1

10

2

10

2

10

2

10

3

10

3

10

3

10 420

40

60

Timestep

80

10 420

40

60

Time Step

80

10 420

40

60

Time Step

80

Figure 15: Highly dynamic estimation error with N = 2 for no prior information (I) (left), no initial prior information (II) (middle), and initial steady-state prior (III) info (right).



t20 t30 t40 t60 tfinal

90

80

70

Flow [scmhr* 104 ]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 33

60

50

40

30 0

50

Axial Position

100

150

Figure 16: Snapshots of simulated transient (thick lines) and corresponding estimated flow profiles (thin lines) with N = 2. extensions are also required to capture non-ideal gas and non-isothermal behavior 30 over large geographical regions and the incorporation of models for other assets such as valves and pressure regulators. 31 We are also interested in exploring the use of state estimation to perform leak detection tasks and to estimate process noise resulting from model mismatch. 7 The proposed optimization-based framework can also be extended to estimate other important parameters such as friction factors and to correct strong errors in input and output measurements. 32

Supporting Information Model nomenclature; set definitions; variables and units; conversion of euler equations to flow and pressure form.

Author Information ∗

Tel.: (608) 890-2111

Email:[email protected]


Page 29 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ORCID Jordan Jalving: 0000-0002-2299-0119 Victor Zavala: 0000-0002-5744-7378 Notes The authors declare no competing financial interest

Acknowledgments This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, under Contract No. DE-AC02-06CH11357 as well as the DOE Office of Electricity Delivery and Energy Reliability’s Advanced Grid Research and Development program (AGR&D). DOE’s Advanced Grid Modeling program under AGR&D directs the research and development of advanced modeling, computational, and control technologies to improve the reliability, resiliency, security, and flexibility of the nation’s electricity system.

References (1) Massachusetts Institute of Technology, Growing Concerns, Possible Solutions: The Interdependency of Natural Gas and Electricity Systems. 2013, (2) Lyons, C.; Litra, G. Gas-Power Interdependence: Knock-on effects of the dash to gas; 2013. (3) Carter, R. G.; Dupont, T. F.; Rachford Jr, H. H. Pack management and transient optimization of natural gas transmission lines. IGRC Conference, Vancouver. 2004. (4) Zlotnik, A.; Chertkov, M.; Backhaus, S. Optimal control of transient flow in natural gas networks. arXiv 2015, 4563–4570. (5) Zavala, V. M. Stochastic optimal control model for natural gas networks. Computers & Chemical Engineering 2014, 64, 103–113. (6) Chiang, N. Y.; Zavala, V. M. Large-scale optimal control of interconnected natural gas and electrical transmission systems. Applied Energy 2016, 168, 226–235. ACS Paragon Plus 29 Environment


(7) Liu, M.; Zang, S.; Zhou, D. Fast leak detection and location of gas pipelines based on an adaptive particle filter. International Journal of Applied Mathematics and Computer Science 2005, 15, 541–550. (8) Reddy, H. P.; Narasimhan, S.; Bhallamudi, S. M. Simulation and State Estimation of Transient Flow in Gas Pipeline Networks Using a Transfer Function Model. Industrial & Engineering Chemistry Research 2006, 45, 3853–3863. (9) Reddy, H. P.; Narasimhan, S.; Bhallamudi, S. M.; Bairagi, S. Leak detection in gas pipeline networks using an efficient state estimator. Part-I: Theory and simulations. Computers and Chemical Engineering 2011, 35, 651–661. (10) Ahmadian Behrooz, H.; Boozarjomehry, R. B. Modeling and state estimation for gas transmission networks. Journal of Natural Gas Science and Engineering 2015, 22, 551–570. (11) Ahmadian Behrooz, H.; Boozarjomehry, R. B. Distributed and decentralized state estimation in gas networks as distributed parameter systems. ISA Transactions 2015, 58, 552–566. (12) He, Y.; Li, S.; Zheng, Y. Distributed state estimation for leak detection in water supply networks. IEEE/CAA Journal of Automatica Sinica 2017, PP, 1–9. (13) Rawlings, J. B.; Bakshi, B. R. Particle filtering and moving horizon estimation. Computers and Chemical Engineering 2006, 30, 1529–1541. (14) Rawlings, J. B.; Ji, L. Optimization-based state estimation: Current status and some new results. Journal of Process Control 2012, 22, 1439–1444. (15) Zavala, V. M.; Biegler, L. T. Optimization-based strategies for the operation of lowdensity polyethylene tubular reactors: Moving horizon estimation. Computers & Chemical Engineering 2009, 33, 379–390. (16) Osiadacz, A. Simulation of transient gas flows in networks. International Journal for Numerical Methods in Fluids 1984, 4, 13–24.


Page 30 of 33

Page 31 of 33 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(17) T. Baumrucker, B.; Biegler, L. MPEC strategies for cost optimization of pipeline operations. 2010, 34, 900–913. (18) Hahn, M.; Leyffer, S.; Zavala, V. M. Mixed-Integer PDE-Constrained Optimal Control of Gas Networks. 2017, (19) Herty, M. Multi scale modeling and nodal control for gas transportation networks. 2015, 4585–4590. (20) Herty, M.; Mohring, J.; Sachers, V. A new model for gas flow in pipe networks. Mathematical Methods in the Applied Sciences 2010, 33, 845–855. (21) Folga, S. Natural Gas Pipeline Technology Overview; 2008; Vol. 25; pp 1163–306. (22) Rao, C. V.; Rawlings, J. B.; Mayne, D. Q. Constrained state estimation for nonlinear discrete-time systems: Stability and moving horizon approximations. IEEE transactions on automatic control 2003, 48, 246–258. (23) Alessandri, A.; Baglietto, M.; Battistelli, G.; Zavala, V. Advances in moving horizon estimation for nonlinear systems. Decision and Control (CDC), 2010 49th IEEE Conference on. 2010; pp 5681–5688. (24) Jalving, J.; Abhyankar, S.; Kim, K.; Hereld, M.; Zavala, V. M. A graph-based computational framework for simulation and optimisation of coupled infrastructure networks. IET Generation, Transmission & Distribution 2017, 1–14. (25) Bezanson, J.; Edelman, A.; Karpinski, S.; Shah, V. B. Julia: A Fresh Approach to Numerical Computing. SIAM Review 2017, 59, 65–98. (26) Lubin, M.; Dunning, I. Computing in Operations Research Using Julia. INFORMS Journal on Computing 2015, 27, 238–248. ´ (27) Santamar´ıa, F. L.; Gomez, J. M. Framework in PYOMO for the assessment and implementation of (as) NMPC controllers. Computers & Chemical Engineering 2016, 92, 93–111. (28) Wächter, A.; Biegler, L. T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. 2006, 106, 25–57. ACS Paragon Plus 31 Environment


(29) Natural Gas Prices,U.S. Energy Information Association. https://www.eia.gov/ dnav/ng/ng_pri_sum_dcu_nus_m.htm, 2017. (30) Chaczykowski, M. Transient flow in natural gas pipeline - The effect of pipeline thermal model. Applied Mathematical Modelling 2010, 34, 1051–1067. (31) Rami, E. G.; Jean-Jacques, B.; Bruno, D.; François, M. Modelling of a pressure regulator. International Journal of Pressure Vessels and Piping 2007, 84, 234–243. ´ (32) Nicholson, B.; Lopez-Negrete, R.; Biegler, L. T. On-line state estimation of nonlinear dynamic systems with gross errors. Computers & Chemical Engineering 2014, 70, 149– 159.


Page 32 of 33

Page 33 of 33

Simulated Gas Network

Pipeline Junction Pressure Measurements

Optimization-based estimator infers dynamic profiles Estimated Dynamics

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 17: For Table of Contents Only


An Optimization-Based State Estimation Framework for Large-Scale

Recommend Documents