# Linear System Analysis and Optimal Control of Natural Gas Dynamics in
Pipeline Networks

Luke S. Baker1, Sachin Shivakumar2, Dieter Armbruster1, Rodrigo B. Platte1,
and Anatoly Zlotnik3 *This study was partially supported by the National
Science Foundation grants No. 1739990 and 1935453 and by the U.S. Department
of Energy’s Advanced Grid Modeling (AGM) project “Dynamical Modeling,
Estimation, and Optimal Control of Electrical Grid-Natural Gas Transmission
Systems”. Research conducted at Los Alamos National Laboratory is done under
the auspices of the National Nuclear Security Administration of the U.S.
Department of Energy under Contract No. 89233218CNA000001.1L. Baker, D.
Armbruster, and R. Platte are with the School of Mathematical and Statistical
Sciences at Arizona State University, Tempe, Arizona, 85281;
<EMAIL_ADDRESS>Sachin Shivakumar is with School for
Engineering of Matter, Transport and Energy at Arizona State University,
Tempe, Arizona, 85281<EMAIL_ADDRESS>3Anatoly Zlotnik is in the Applied
Mathematics & Plasma Physics Group at Los Alamos National Laboratory, Los
Alamos, New Mexico, 87545<EMAIL_ADDRESS>
###### Abstract

We derive a linear system of ordinary differential equations (ODEs) to
approximate the dynamics of natural gas in pipeline networks. Although a
closed-form expression of the eigenvalues of the state matrix does not
generally exist, the poles of an irrational transfer function corresponding to
the linearized partial differential equations are used to approximate the
eigenvalues of the ODE system. Our analysis qualitatively demonstrates that
the eigenvalues of the state matrix of the entire network system are “pipeline
separable” in the sense that the eigenvalues are dominated by the individual
pipeline parameters and not the incidence connectivity of the network graph.
The linear system is used as the dynamic constraints of a linear optimal
control problem (OCP) to design the control actions of compressor units to
minimize the energy that they expend. The motivation of this work is to reduce
the computational complexity of optimizing gas dynamics in large networks to
meet the unpredictable and highly variable demand from electric generators.
The linear and corresponding nonlinear OCPs are discretized in time to obtain
linear and nonlinear optimization problems, which are demonstrated on a test
network to illustrate the validity of linear programming. Moreover, an
analytical bound on the error between the solutions of the linear and
nonlinear flow dynamics is presented using Lyapunov functions and verified
computationally by plotting the error against the size of the flow variation
around the steady-state solution.

##  I Introduction

Natural gas is the primary source of energy used to generate electricity in
the United States, but substantial investments are being made to transition
the U.S. economy from fossil fuels such as natural gas and coal to cleaner and
more sustainable resources. The U.S. Energy Information Administration
projects a rapid increase in the installation of renewable energy resources
over the next 30 years as coal transitions into retirement [1]. However, the
reliability of renewable energy sources, such as wind and sun, is a persistent
issue because the energy output from these sources is generally unpredictable
and highly variable throughout the course of a day. During hours of peak
electricity demand, electric generators balance temporal shortages of
renewable energy by using natural gas-fired power plants. Thus, there is a
temporal variability in gas withdrawal from associated pipelines that creates
highly transient flows throughout the natural gas infrastructure. Pipeline
system operators must consistently moderate pressure, mass flows, and
compressor activity as part of intra-day planning and operations.

The transmission of natural gas through networks of pipelines has been studied
in steady-state [2, 3, 4, 5, 6] and transient operations [7, 8, 9]. A major
goal is to design compressor activity to minimize the expenditure of
compressor energy [10] or maximize the economic value generated by the
pipeline system [11], while satisfying the physics of gas flow throughout the
network in such a way that engineering limits on pressure and compression are
satisfied. In steady-state operation, the flow of gas in the network is
balanced, so that the totality of inflows from processing plants and supply
stations equates to the outflows from withdrawal stations. Steady-state
pipeline flows are described using simple time-invariant algebraic equations
that relate pressure drop in the direction of flow to mass flow along each
pipeline. In transient operation, the computational complexity increases
significantly because the flow in each pipeline is governed by a system of
nonlinear partial differential equations (PDEs) [12, 13].

Optimization of gas network dynamics is usually implemented in a digital
computer by using discretization methods in space and time to approximate the
continuous PDEs with algebraic equations. Since these optimization problems
are large, nonlinear, and nonconvex with nonlinear algebraic constraints,
model reduction methods have been proposed to reduce the complexity in both
steady-state and transient operations [14, 15, 16]. Green’s functions are used
in [17] to form a linear optimization program, which is applied on a test
network to demonstrate a reduction of optimization time by two orders of
magnitude in comparison to the nonlinear optimization program. A mixed-integer
linear programming approach based on piecewise linearization of the nonlinear
terms was applied in [18] for the optimal control of transient flows. A
similar linearization approach for the coupled gas network and electric power
grid is presented in [19] for steady-state gas flows. These linearization
methods usually require additional discretization points to accurately
interpolate the nonlinear terms with piecewise linear segments, which could
significantly increase the size of the network.

The key contribution of our study is the development of a linear state
equation and linear program (LP) that approximate the nonlinear state equation
and nonlinear program (NLP) within the same state space, i.e., without the
addition of interpolating piecewise segments. The linear system is derived in
a way that can be utilized for the design of feedback controllers and the
investigation of asymptotic stability and transient responses. We begin by
defining the PDE-constrained optimal control problem (OCP) and discretize it
in space using a finite volume method to form a nonlinear control system of
ordinary differential equations (ODEs) that is subsequently written in matrix-
vector form using the edge incidence matrix of the network [20]. This
discretization approach was used to formulate NLPs in previous studies [9,
21], with applications in gas reservoir storage [22] and coordinated
scheduling of natural gas and electric power infrastructures [23]. We derive
the linear system and LP from the nonlinear system and NLP using linearization
techniques around a steady-state solution.

The rest of this paper is organized as follows. In Section II, the system of
nonlinear PDEs that dictate the flow of gas in the network is presented and
used to formulate the PDE-constrained OCP. In Section III, the PDE system is
discretized in space to form a nonlinear ODE state equation and ODE-
constrained OCP. The linear state equation and linear ODE-constrained OCP are
derived in Section IV. The eigenvalues of the linear state matrix are compared
to the poles of an irrational transfer function in Section V to gain insight
into the transient behavior of the network system. The nonlinear and linear
programs are implemented in Section VI using Euler’s approximation for the
time derivatives that arise in the state equations. Section VII demonstrates
the performance of the two programs on a test network that was used in a
previous study [24]. The error between the solutions of the linear and
nonlinear systems is analyzed computationally and analytically in Section
VIII. Concluding remarks are made in Section IX and a proof of the error bound
is presented in the Appendix.

##  II Network Flow Control Formulation

A gas pipeline network is modeled as a connected and directed graph
$(\mathcal{E},\mathcal{V})$ consisting of edges $\mathcal{E}=\\{1,\dots,E\\}$
and nodes $\mathcal{V}=\\{1,2,\dots,V\\}$, where $E$ and $V$ denote the
numbers of edges and nodes of the graph. Edges represent pipelines and nodes
represent junctions or stations where gas can be injected into or withdrawn
from the network. It is assumed that the edges and nodes are ordered within
their sets according to their integer labels. The symbol $k$ is reserved for
indexing edges in $\mathcal{E}$ and the symbols $i$ and $j$ are reserved for
indexing nodes in $\mathcal{V}$. Supply nodes
$\mathcal{V}_{s}\subset\mathcal{V}$ and withdrawal nodes
$\mathcal{V}_{w}\subset\mathcal{V}$ are assumed to be disjoint sets that
partition $\mathcal{V}$, i.e.,
$\mathcal{V}_{s}\cup\mathcal{V}_{w}=\mathcal{V}$ and
$\mathcal{V}_{s}\cap\mathcal{V}_{w}=\emptyset$. It is assumed that supply
nodes are ordered in $\mathcal{V}$ before withdrawal nodes so that $i<j$ for
all $i\in\mathcal{V}_{s}$ and $j\in\mathcal{V}_{w}$. The graph is directed by
assigning a positive flow direction along each edge. Mass flux and velocity
have positive values along an edge if gas physically flows along this edge in
the prescribed positive direction of the graph. We assume that gas flows in
the positive orientation of the graph. The notation $k:i\mapsto j$ means that
edge $k\in\mathcal{E}$ is directed from node $i\in\mathcal{V}$ to node
$j\in\mathcal{V}$. For each node $j\in\mathcal{V}$, we define (potentially
empty) incoming and outgoing sets of pipelines by
${}_{\mapsto}j=\\{k\in\mathcal{E}|k:i\mapsto j\\}$ and
$j_{\mapsto}=\\{k\in\mathcal{E}|k:j\mapsto i\\}$, respectively.

For each pipe $k\in\mathcal{E}$, the flow variables are density
$\rho_{k}(t,x)$ and mass flux $\varphi_{k}(t,x)$ for $t\in[0,T]$ and
$x\in[0,\ell_{k}]$, where $T$ denotes the time horizon and $\ell_{k}$ denotes
the length of the pipe. Assuming that the pipe is horizontal, the flow is
isothermal, and the transients do not excite shock waves, then the flow
through edge $k\in\mathcal{E}$ may be governed by the semilinear hyperbolic
PDE system [8, 25]

$\displaystyle\partial_{t}\rho_{k}+\partial_{x}\varphi_{k}$ $\displaystyle=$
$\displaystyle 0,$ (1)
$\displaystyle\partial_{t}\varphi+\sigma^{2}\partial_{x}\rho_{k}$
$\displaystyle=$
$\displaystyle-\frac{\lambda_{k}}{2D_{k}}\frac{\varphi_{k}|\varphi_{k}|}{\rho_{k}},\quad$
(2)

where $D_{k}$ and $\lambda_{k}$ are the diameter and the friction factor of
the pipe, respectively, and $\sigma$ is the speed of sound through natural
gas. The compressibility factor of the gas is assumed to be constant so that
the equation of state is ideally given by $p_{k}=\sigma^{2}\rho_{k}$, where
$p_{k}$ is the pressure. These assumptions are often adopted for analyzing
flows in natural gas systems [26, 27].

Forces of friction between the interior wall of a pipe and gas flowing through
the pipe cause pressure to decrease in the direction of flow. Compressor
stations are strategically installed throughout the network to increase the
pressure of gas at these locations to within limits required for
transportation. We define $\mathcal{C}\subset\mathcal{E}$ to be the set of
edges $k\in\mathcal{E}$ that are adjacent to a compressor station and we
simply speak of the compressor station $k\in\mathcal{C}$. We assume that each
compressor station $k\in\mathcal{C}$ is located at the entrance of its
adjacent edge $k\in\mathcal{E}$ with respect to the positive flow direction.
The action of the compressor station $k\in\mathcal{C}$ is modeled with the
control input $\mu_{k}(t)$. In particular, the pressure of gas leaving the
compressor unit $(k:i\mapsto j)\in\mathcal{C}$ and entering the pipeline
$(k:i\mapsto j)\in\mathcal{E}$ is $\mu_{k}(t)$ times larger than the pressure
of gas leaving node $i\in\mathcal{V}$ and entering compressor
$k\in\mathcal{C}$. To simplify notation, we define $\mu_{k}=1$ for all
$k\in\mathcal{E}\setminus\mathcal{C}$.

Figure 1: Configuration of a pipeline segment $k:i\mapsto j$ with
$i\in\mathcal{V}_{s}$ and $j\in\mathcal{V}_{w}$ with state and control
variables indicated.

The density of gas injected into the network at the supply node
$i\in\mathcal{V}_{s}$ is specified by the boundary condition profile
$\bm{s}_{i}(t)$ (kg/m3). The amount of gas withdrawn from the network at each
withdrawal node $j\in\mathcal{V}_{w}$ is specified by the mass outflow profile
$\bm{w}_{j}(t)$ (kg/s). For each $j\in\mathcal{V}_{w}$, define the nodal
density variable $\bm{\rho}_{j}(t)$. All of the nodal quantities in this work
are identified with bold symbols. Inlet and outlet edge variables are defined
by attaching superscripts “0” and “$\ell$”, respectively, to the associated
edge variables. For example, $\varphi_{k}^{0}(t)=\varphi_{k}(t,0)$ and
$\varphi^{\ell}_{k}(t)=\varphi_{k}(t,\ell_{k})$. Define $\chi_{k}=\pi
D_{k}^{2}/4$ to be the cross-sectional area of the pipeline $k\in\mathcal{E}$.
The boundary conditions of the network dynamics are given by

$\displaystyle\rho_{k}(t,0)$ $\displaystyle=$
$\displaystyle\mu_{k}(t)\bm{s}_{i}(t),\quad\rho_{k}(t,\ell_{k})=\bm{\rho}_{j}(t),$
(3) $\displaystyle\rho_{k}(t,0)$ $\displaystyle=$
$\displaystyle\mu_{k}(t)\bm{\rho}_{i}(t),\quad\rho_{k}(t,\ell_{k})=\bm{\rho}_{j}(t),$
(4) $\displaystyle\bm{w}_{j}(t)$ $\displaystyle=$
$\displaystyle\sum_{k\in_{\mapsto}j}\chi_{k}\varphi_{k}^{\ell}(t)-\sum_{k\in
j_{\mapsto}}\chi_{k}\varphi_{k}^{0}(t),$ (5)

where (3) is defined for $k:i\mapsto j$ with $i\in\mathcal{V}_{s}$, (4) is
defined for $k:i\mapsto j$ with $i\in\mathcal{V}_{w}$, and (5) is defined for
$j\in\mathcal{V}_{w}$. The conditions in (3)-(4) represent the effects of
compression and the conditions in (5) represent the conservation of mass flow
through withdrawal nodes. The initial condition of density and mass flux in
the network is assumed to be a steady-state solution given by [6]

$\displaystyle\rho_{k}(0,x)=\overline{\rho}_{k}(x),\qquad\varphi(0,x)=\overline{\varphi}_{k},$
(6)

where $\overline{\varphi}_{k}$ is constant for each $k\in\mathcal{E}$. We
provide details on the initial condition for the discretized system in the
following section. We assume standard conditions for well-posedness [28], and
specifically that the boundary conditions are smooth, slowly-varying, bounded
in their respective domains, and compatible with the initial condition to
ensure the existence of a smooth, slowly-varying, bounded solution. The flow
of natural gas in the network is defined by (1)-(6).

Gas network operators require pressure and compression to be within
engineering limits to ensure the safety of transportation and the quality of
gas delivered to customers. These limitations are modeled for all
$k\in\mathcal{E}$ with inequality constraints of the form

$p^{\min}_{k}\leq\sigma^{2}\rho_{k}\leq p^{\max}_{k},\quad 1\leq\mu_{k}\leq
2,$ (7)

where $p^{\min}_{k}$ and $p^{\max}_{k}$ are specified bounds on pressure. The
compressor control variables are designed to minimize the energy that they
expend. The accumulated energy is given by [2, 9]

$J=\sum_{k\in\mathcal{E}}\int_{0}^{T}c_{k}\varphi_{k}^{0}(t)\left((\mu_{k}(t))^{(\nu-1)/\nu}-1\right)dt,$
(8)

where $c_{k}$ is related to the efficiency of the compressor $\mu_{k}$ and
$\nu$ is the isentropic exponent of natural gas [29]. Therefore, the control
design is defined by the PDE-constrained OCP

$\begin{array}[]{ll}\text{min}&J\triangleq\text{compressor energy in
\eqref{J}},\\\ \text{s.t.}&\text{dynamic constraints: \eqref{1}-\eqref{2}},\\\
&\text{boundary conditions: \eqref{bc1}-\eqref{bc}},\\\ &\text{initial
condition: \eqref{ic}},\\\ &\text{inequality constraints:
\eqref{ineq}}.\end{array}$ (9)

The decision variables are densities, mass fluxes, and compressor ratios
throughout the network. Other OCPs of interest, such as carbon dioxide
reduction [30, 31], may be defined similarly by adjusting the objective
function in (8).

##  III Network Flow Control Discretization

The intial-boundary value system of PDEs from the previous section will be
discretized in space to obtain an initial-value system of ODEs using a popular
finite volume method for natural gas networks [20]. Discretization will be
formalized by refining the graph of the gas network. A graph refinement
$(\hat{\mathcal{E}},\hat{\mathcal{V}})$ of the graph
$(\mathcal{E},\mathcal{V})$ is made by adding auxiliary withdrawal nodes to
$\mathcal{V}$ that subdivide the edges of $\mathcal{E}$ so that
$\ell_{k}\leq\ell$ for all $k\in\hat{\mathcal{E}}$, where $\ell\leq 10$ (km)
is sufficiently small [15]. The refined graph inherits the prescribed
direction of the parent graph. For sufficiently fine network refinement, the
relative difference of the flow variables between adjacent nodes is small in
magnitude by continuity of the flow variables (assuming well-posedness). We
assume that the graph has been sufficiently refined and that the hats may be
omitted moving forward.

The system of ODEs is obtained by integrating the dynamic equations in (1)-(2)
along the length of each refined pipeline segment so that

$\displaystyle\int_{0}^{\ell}\partial_{t}\rho+\partial_{x}\varphi dx$
$\displaystyle=$ $\displaystyle 0,$
$\displaystyle\int_{0}^{\ell}\partial_{t}\varphi+\sigma^{2}\partial_{x}\rho
dx$ $\displaystyle=$
$\displaystyle-\frac{\lambda}{2D}\int_{0}^{\ell}\frac{\varphi|\varphi|}{\rho}dx,$

where edge subscripts have been removed for readability. The above integrals
of space derivatives are evaluated using the fundamental theorem of calculus.
The remaining integrals are evaluated by approximating pipeline density with
outlet density and pipeline flux with inlet flux. These approximations are
independent of $x$ and may be factored out of the integrals. The above
equations become

$\displaystyle\ell\dot{\rho}^{\ell}$ $\displaystyle=$
$\displaystyle\varphi^{0}-\varphi^{\ell},$ (10)
$\displaystyle\ell\dot{\varphi}^{0}+\sigma^{2}\left(\rho^{\ell}-\rho^{0}\right)$
$\displaystyle=$
$\displaystyle-\frac{\lambda\ell}{2D}\frac{\varphi^{0}\left|\varphi^{0}\right|}{\rho^{\ell}},$
(11)

where a dot above a variable represents the time-derivative of the variable.
We now write the discretized system in matrix form. Define the $E\times E$
diagonal matrices $L$, $K$, and $X$ with diagonal entries $L_{kk}=\ell_{k}$,
$K_{kk}=\lambda_{k}/(2D_{k})$, and $X_{kk}=\chi_{k}$. Define the time-varying
(transposed) incidence matrix $\Xi$ of size $E\times V$ componentwise by

$\displaystyle\Xi_{ki}=\begin{cases}-\mu_{k}(t),&\text{edge $k\in i_{\mapsto}$
leaves node $i$,}\\\ 1,&\text{edge $k\in_{\mapsto}i$ enters node $i$},\\\
0,&\text{else.}\end{cases}$ (12)

Define the $E\times r$ submatrix $N$ of $\Xi$ by the removal of columns
$i\in\mathcal{V}_{w}$ and the $E\times(V-r)$ submatrix $M$ of $\Xi$ by the
removal of columns $i\in\mathcal{V}_{s}$, where $r$ denotes the number of
supply nodes. Define the signed matrix $Q=\text{sign}(M)$ (which is well-
defined by the inequalities in (7)). Define the positive and negative parts of
$Q$ by $Q_{\ell}$ and $Q_{0}$, respectively, so that $Q=(Q_{\ell}+Q_{0})/2$
and $|Q|=(Q_{\ell}-Q_{0})/2$, where $|A|$ denotes the componentwise absolute
value of a matrix $A$.

Define inlet and outlet edge mass flux vectors by
$\varphi^{0}=(\varphi_{1}^{0},\dots,\varphi_{E}^{0})^{T}$ and
$\varphi^{\ell}=(\varphi_{1}^{\ell},\dots,\varphi_{E}^{\ell})^{T}$. Define the
vector of densities at supply nodes
$\bm{s}=(\bm{s}_{1},\dots,\bm{s}_{r})^{T}$, the vector of mass withdrawals
$\bm{w}=(\bm{w}_{r+1},\dots,\bm{w}_{V})^{T}$, and the vector of densities at
withdrawal nodes $\bm{\rho}=(\bm{\rho}_{r+1},\dots,\bm{\rho}_{V})^{T}$, where
the subscripts of the entries are indexed according to the node labels in
$\mathcal{V}$. Applying the above matrix definitions, the discretized
equations in (10)-(11) together with the boundary conditions in (3)-(5) become

$\displaystyle LQ_{\ell}\dot{\bm{\rho}}$ $\displaystyle=$
$\displaystyle\varphi^{0}-\varphi^{\ell},$ (13) $\displaystyle
L\dot{\varphi}^{0}+\sigma^{2}\left(M\bm{\rho}+N\bm{s}\right)$ $\displaystyle=$
$\displaystyle-LK\frac{\varphi^{0}\odot|\varphi^{0}|}{Q_{\ell}\bm{\rho}},$
(14) $\displaystyle\bm{w}$ $\displaystyle=$ $\displaystyle
Q_{\ell}^{T}X\varphi^{\ell}+Q_{0}^{T}X\varphi^{0},$ (15)

where $\odot$ is the Hadamard (componentwise) product and the ratio of vectors
on the right side of (14) is componentwise as well. Multiplying both sides of
(13) on the left by $Q_{\ell}^{T}X$ and using (15), we may combine (13) and
(15) to form the equation
$Q_{\ell}^{T}XLQ_{\ell}\dot{\bm{\rho}}=(Q^{T}X\varphi^{0}-\bm{w})$, where we
have used $Q=(Q_{0}+Q_{\ell})$. Therefore, outlet flux is a dependent variable
and we define the state vector of inlet flux by $\varphi=\varphi^{0}$.
Moreover, we define $I=Q_{\ell}$ for simplicity. Note that, for a connected
tree network with one supply node, the nodes and edges of the network may be
ordered so that $I$ is the identity matrix. The above equations become

$\displaystyle\dot{\bm{\rho}}$ $\displaystyle=$
$\displaystyle(I^{T}XLI)^{-1}\left(Q^{T}X\varphi-\bm{w}\right),$ (16)
$\displaystyle\dot{\varphi}$ $\displaystyle=$
$\displaystyle-\sigma^{2}L^{-1}\left(M\bm{\rho}+N\bm{s}\right)-K\frac{\varphi\odot|\varphi|}{I\bm{\rho}}.$
(17)

The matrices $L$ and $I^{T}XLI$ are invertible because $L$ and $X$ are
invertible and $I$ has full rank. In fact, each row $k$ of $I$ contains
exactly one nonzero component with $I_{kj}=1$ if and only if
$k\in$${}_{\mapsto}j$. Using the additional fact that $X$ and $L$ are
diagonal, it can be shown that the mass matrix $I^{T}XLI$ is diagonal with
positive diagonal components given by
$(I^{T}XLI)_{jj}=\sum_{k\in_{\mapsto}j}\chi_{k}\ell_{k}$ for
$j\in\mathcal{V}_{w}$. Therefore, the matrices $L$ and $I^{T}XLI$ may be
easily inverted to obtain the traditional control system presented above. The
state variables are densities $\bm{\rho}$ at withdrawal nodes and fluxes
$\varphi$ at the inlets of the edges. The compressor actuators are contained
in the matrices $M$ and $N$. The other matrices are known and constant.

The steady-state initial condition in (6), sampled at the refined nodes of the
network, is the solution of the time-invariant flow equations

$\displaystyle
Q^{T}X\overline{\varphi}=\overline{\bm{w}},\qquad\sigma^{2}\left(\overline{M}\overline{\bm{\rho}}+\overline{N}\overline{\bm{s}}\right)=-LK\frac{\overline{\varphi}\odot|\overline{\varphi}|}{I\overline{\bm{\rho}}},$
(18)

where $\overline{\bm{s}}=\bm{s}(0)$ and $\overline{\bm{w}}=\bm{w}(0)$.
Overlines attached to state variables, actuation matrices, and parameters are
used throughout to denote a time-invariant steady-state solution. We assume
that $\overline{M}$, $\overline{N}$, $\overline{\varphi}$, and
$\overline{\bm{\rho}}$ are optimally determined to minimize compressor energy
as in (8) while satisfying flow requirements in (18) and inequality
constraints as in (20). The initial condition of the system in (16)-(17) is
defined to be

$\displaystyle\bm{\rho}(0)=\overline{\bm{\rho}},\qquad\varphi(0)=\overline{\varphi}.$
(19)

The discretized system in (16)-(17) is numerically consistent in the sense
that its dynamics approach the continuous dynamics in (1)-(5) as the maximum
length of the edges of the refined network approaches zero. Pressure and
compression inequality constraints in (7) reduce to

$\bm{p}^{\min}_{j}\leq\sigma^{2}\bm{\rho}_{j}\leq\bm{p}^{\max}_{j},\quad
1\leq\mu_{k}\leq 2,$ (20)

where $\bm{p}^{\min}_{j}$ and $\bm{p}^{\max}_{j}$ are specified for each node
$j\in\mathcal{V}_{w}$. The ODE-constrained OCP is formulated as

$\begin{array}[]{ll}\text{min}&J\triangleq\text{compressor energy in
\eqref{J}},\\\ \text{s.t.}&\text{dynamic constraints:
\eqref{partial1}-\eqref{partial2}},\\\ &\text{initial condition:
\eqref{Disic}},\\\ &\text{inequality constraints:
\eqref{Disineq}}.\end{array}$ (21)

##  IV Linearization

In this section, we use classical linearization techniques to derive a linear
system that models the flow of gas in the network as the dynamics vary around
the steady-state solution. Flow variables, compressor actuators, and time-
varying supply and withdrawal profiles in (16)-(17) are first written in terms
of variations around their steady-state counterparts. The system of variations
is given by

$\displaystyle\dot{\bm{\rho}}$ $\displaystyle=$
$\displaystyle(I^{T}XLI)^{-1}\left(Q^{T}X\varphi-\bm{w}\right),$ (22)
$\displaystyle\dot{\varphi}$ $\displaystyle=$
$\displaystyle-\sigma^{2}L^{-1}\left(M\bm{\rho}+\overline{M}\bm{\rho}+M\overline{\bm{\rho}}+N\bm{s}+\overline{N}\bm{s}+N\overline{\bm{s}}\right)$
(23)
$\displaystyle+K\left(\frac{\overline{\varphi}\odot|\overline{\varphi}|}{I\overline{\bm{\rho}}}-\frac{\left(\overline{\varphi}+\varphi\right)\odot\left|\overline{\varphi}+\varphi\right|}{I\left(\overline{\bm{\rho}}+\bm{\rho}\right)}\right).$

The variation variables are expressed with the same symbols as their original
counterparts. All of the flow variables moving forward represent variations
around the steady-state. However, we adhere to convention and depict the non-
variation variables in the figures that follow. The system of variations in
(22)-(23) is equivalent to the original system in (16)-(17), where the
equivalence relation is represented as a translation of the equilibrium from
the origin to the steady-state solution.

We first rewrite the control input products $M\overline{\bm{\rho}}$ and
$N\overline{\bm{s}}$ in (23) into the standard form of a linear system. Each
row $k$ of $M$ and $N$ contain exactly one nonzero entry given by
$M_{ki}=-\mu_{k}$ and $N_{ki}=-\mu_{k}$ if and only if $k\in
i_{\mapsto}\cap\mathcal{C}$ with $i\in\mathcal{V}_{w}$ and
$i\in\mathcal{V}_{s}$, respectively. Define the $C$-dimensional vector
$\mu=(\\{\mu_{k}\\}_{k\in\mathcal{C}})$, where $C$ is the number of active
compressor stations (cardinality of $\mathcal{C}$). Define the matrix
$\overline{P}$ componentwise by

$\overline{P}_{ki}=\begin{cases}\overline{\bm{s}}_{i},&\text{if edge $k\in
i_{\mapsto}$ leaves node $i\in\mathcal{V}_{s}$,}\\\
\overline{\bm{\rho}}_{i},&\text{if edge $k\in i_{\mapsto}$ leaves node
$i\in\mathcal{V}_{w}$,}\\\ 0,&\text{otherwise,}\end{cases}$

and the $E\times C$ submatrix $\overline{B}$ of $\overline{P}$ by the removal
of every column $i$ of $\overline{P}$ for which $k\in i_{\mapsto}$ with
$k\not\in\mathcal{C}$. Then we have
$-(M\overline{\bm{\rho}}+N\overline{\bm{s}})=\overline{B}\mu$. Define the
$E\times E$ diagonal matrices $\overline{\alpha}$ and $\overline{\beta}$ that
result in linearizing the nonlinear frictional term in (23) by

$\displaystyle\overline{\alpha}=K\text{diag}\left(\frac{\overline{\varphi}\odot\left|\overline{\varphi}\right|}{I\overline{\bm{\rho}}^{2}}\right),\qquad\overline{\beta}=-2K\text{diag}\left(\frac{\left|\overline{\varphi}\right|}{I\overline{\bm{\rho}}}\right).$
(24)

Linearizing (22)-(23) in $\bm{\rho}$, $\varphi$, $\mu$, $\bm{s}$, and $\bm{w}$
around the origin gives the linear time-invariant system

$\displaystyle\dot{\begin{bmatrix}\bm{\rho}\\\
\varphi\end{bmatrix}}=A\begin{bmatrix}\bm{\rho}\\\
\varphi\end{bmatrix}+B\mu+W\begin{bmatrix}\bm{s}\\\ \bm{w}\end{bmatrix},$ (25)

where the state, input, and disturbance matrices are given in block matrix
form by

$\displaystyle A=\begin{bmatrix}0&(I^{T}XLI)^{-1}Q^{T}X\\\
\overline{\alpha}I-\sigma^{2}L^{-1}\overline{M}&\overline{\beta}\end{bmatrix}$

and

$B=\sigma^{2}L^{-1}\overline{B},\quad W=\begin{bmatrix}0&-(I^{T}XLI)^{-1}\\\
-\sigma^{2}L^{-1}\overline{N}&0\end{bmatrix}.$

The entries with zeros denote zero matrices of approapriate dimensions. The
initial condition of (25) is given by

$\bm{\rho}(0)=0,\qquad\varphi(0)=0.$ (26)

The initial-value system in (25)-(26) describes the gas dynamics for all
variations around the steady-state that are confined to a region in which
linearization is applicable. In subsequent sections, we will analyze this
applicable linearization region. Linearizing (8) around the steady-state gives
the linear objective

$\displaystyle
J_{\text{lin}}=\sum_{k\in\mathcal{E}}c_{k}\int_{0}^{T}\varphi_{k}(t)\left(\overline{\mu}_{k}^{(\gamma-1)/\gamma}-1\right)$
$\displaystyle+\left(\frac{\gamma-1}{\gamma}\overline{\varphi}_{k}\overline{\mu}_{k}^{-1/\gamma}\right)\mu_{k}(t)dt.$
(27)

The inequality constraints in (20) are linear in the state and control
variables. Translating these inequalities from the steady-state to the origin
results in

$\bm{p}^{\min}_{j}\leq\sigma^{2}\left(\bm{\rho}_{j}+\overline{\bm{\rho}}_{j}\right)\leq\bm{p}^{\max}_{j},\quad
1\leq\mu_{k}+\overline{\mu}_{k}\leq 2.$ (28)

The linear ODE-constrained OCP is given by

$\begin{array}[]{ll}\text{min}&J_{\text{lin}}\triangleq\text{compressor energy
in \eqref{Jlin}},\\\ \text{s.t.}&\text{dynamic constraints:
\eqref{linsys}},\\\ &\text{initial condition: \eqref{linic}},\\\
&\text{inequality constraints: \eqref{linineq}}.\end{array}$ (29)

##  V Transient Response

Although the gas network is physically stable because frictional effects
strongly dissipate kinetic energy, a quantitative analysis of the eigenvalues
of the state matrix $A$ is beneficial for the design of controllers and the
analysis of transient and frequency responses. However, a closed-form
representation of the eigenvalues may not exist for even the simplest of
networks, such as a refined pipeline. We first provide a result that
characterizes the “center of gravity” of the eigenvalues of $A$.

Proposition 1. Let $\\{\zeta_{m}\\}_{m=1}^{E+V_{w}}$ be the set of eigenvalues
of $A$ in (25), where $E$ and $V_{w}$ are the cardinalities of $\mathcal{E}$
and $\mathcal{V}_{w}$. Then

$\sum_{m=1}^{E+V_{w}}\zeta_{m}=-\sum_{k=1}^{E}\frac{\lambda_{k}|\overline{\varphi}_{k}|}{D_{k}(I\overline{\bm{\rho}})_{k}}.$
(30)

Proof. The trace formula for the sum of the eigenvalues of the matrix $A$
gives

$\sum_{m=1}^{E+V_{w}}\zeta_{m}=\sum_{k=1}^{E}\overline{\beta}_{kk}.$

The result follows from (24). $\square$

The center of gravity of the eigenvalues of $A$ is defined to be the average
value $c=\sum_{m=1}^{E+V_{w}}\zeta_{m}/(E+V_{w})$. Proposition 1 shows that
the center of gravity depends on $\overline{\beta}$, but not on
$\overline{\alpha}$. Moreover, the center of gravity does not explicitly
depend on the connectivity of the network graph described by the incidence
submatrices $\overline{M}$ and $\overline{N}$. The left and right sides of
Figure 2 depict the eigenvalues of $A$ for a single (discretized) pipeline
with two different mass outflow boundary conditions. The figure also depicts a
subset of poles of an irrational transfer matrix $G(s)$ obtained by
linearizing the PDEs in (1)-(2) and transforming the linear PDEs into the
frequency domain [32, 33, 34]. Similarly as in [32], the Laplace transforms of
the linearized approximations of $\rho(t,x)$ and $\varphi(t,x)$ in (1)-(2),
denoted by $P(s,x)$ and $\Phi(s,x)$, have the transfer matrix representation

$\displaystyle\left[\begin{matrix}P^{\ell}\\\
\Phi^{0}\end{matrix}\right]=G\left[\begin{matrix}P^{0}\\\
\Phi^{\ell}\end{matrix}\right],\quad
G=\left[\begin{matrix}\text{sech}\left(\frac{L}{\sigma}z\right)&-\frac{z_{c}}{\sigma}\text{tanh}\left(\frac{L}{\sigma}z\right)\\\
\frac{\sigma}{z_{c}}\text{tanh}\left(\frac{L}{\sigma}z\right)&\text{sech}\left(\frac{L}{\sigma}z\right)\end{matrix}\right].$

Here, $z_{c}(s)=\sqrt{(s-\beta)/s}$ and $z(s)=\sqrt{s(s-\beta)}$, where
$\beta$ is a scalar approximation of the diagonal elements of
$\overline{\beta}$. It was illustrated in [32] that just one of the diagonal
elements of $\overline{\beta}$ may be used to approximate the total dynamics
of the matrix $\overline{\beta}$.

Although the eigenvalues of $A$ and poles of $G(s)$ are, in general, not
equal, the poles have a simple closed formula that can be used with caution to
gain insight into the structure of the eigenvalues of $A$. The poles of $G(s)$
are given by

$\zeta^{\pm}_{m}=\frac{\beta}{2}\pm{\bf
j}\sqrt{\left(\frac{\pi\sigma}{2\ell}\right)^{2}(2m+1)^{2}-\left(\frac{\beta}{2}\right)^{2}},$
(31)

for $m=0,1,\dots$. Here, $\bf{j}$ denotes the imaginary unit. The poles of
$G(s)$ in Figure 2 are constructed using the average value
$\beta=\sum_{k=1}^{E}\overline{\beta}_{kk}/E$ in (31). From Proposition 1, it
can be shown that this average value of the diagonal elements of
$\overline{\beta}$ results in $\beta=[(E+V_{w})/E]c$. In the case of a single
pipeline with one supply node, this result simplifies to $\beta=2c$ (since
$E=V_{w}$). Thus, the imaginary asymptote of the poles of $G(s)$ intersects
the real axis at the center of gravity of the eigenvalues of $A$ for this
average choice of $\beta$. This is demonstrated in Figure 2, where the poles
and eigenvalues gravitate toward the imaginary asymptote $c=\beta/2$.

Figure 2: Eigenvalues ($\ast$) of $A$ and poles ($\circ$) of $G(s)$ for a
single pipeline with 30 refined edges. The pipeline parameters are $\ell=100$
(km), $D=0.75$ (m), $\lambda=0.01$, and $\sigma=377$ (m/s).

All of the poles of $G(s)$ have negative real parts, irrespective of the
pipeline parameters, because the diagonal elements of $\overline{\beta}$ are
always negative. This observation agrees with the physical stability of the
system. However, the time response may vary significantly for different
parameter values. The settling time $t_{s}$ is defined to be the time elapsed
from the application of an ideal instantaneous step input to the time at which
the system has entered and will remain within a specified percent error of the
final state. The settling time of a transient input depends largely on the
real parts of one or more of the eigenvalues of $A$. For our purpose, we
approximate the settling time by $t_{s}=1/\min|\text{Re}(\lambda_{m})|$. We
investigate how the settling time changes as the parameters of the pipeline
vary for two cases.

  1. 1.

Varying $\ell$. First, assume that $\ell<\pi\sigma/|\beta|$. In this case, all
of the poles of $G(s)$ lie on the asymptote $c=\beta/2$ and the settling time
is approximated with $t_{s}=1/|c|$. Second, assume that
$\ell>\pi\sigma/|\beta|$. In this case, there are at least two purely real
poles, the largest of all of which is given by
$\zeta_{1}^{+}=[\beta/2+\sqrt{(\beta/2)^{2}-(\pi\sigma/2\ell)^{2}}]$. As
$\ell$ increases without bound, the number of purely real poles of $G(s)$
increases without bound and $\zeta_{1}^{+}$ increases toward the origin of the
complex plane. This indicates that the settling time $t_{s}=1/|\zeta_{1}^{+}|$
increases as $\ell$ increases, and that the system will never settle
everywhere in the theoretical case of an infinitely long pipeline.

  2. 2.

Varying $\overline{\varphi}$. For a single pipeline, the values
$\overline{\varphi}_{k}$ for $k\in\mathcal{E}$ in (30) are all the same value,
which we denote by $\overline{\varphi}$. The poles of $G(s)$ are, therefore,
given by

$\displaystyle\zeta_{m}^{\pm}=b|\overline{\varphi}|\pm{\bf
j}\sqrt{\left(\frac{\pi\sigma}{2\ell}\right)^{2}(2m+1)^{2}-(b|\overline{\varphi}|)^{2}},$

where $b=-\lambda/(2DE)\sum_{k=1}^{E}1/(I\overline{\bm{\rho}})_{k}$. As
$|\overline{\varphi}|\rightarrow 0$, it is clear that all of the poles, hence
the center of gravity, approach the imaginary axis of the complex plane. Thus,
as $|\overline{\varphi}|\rightarrow 0$, the friction term in the state
equation goes to zero and the variations are undamped waves, leading to an
infinite settling time. As $|\overline{\varphi}|$ theoretically increases
without bound, the center of gravity $c$ decreases without bound along the
real axis of the complex plane. The decrease of the center of gravity as
$|\bm{w}|$, hence $|\overline{\varphi}|$, increases is illustrated in Figure
2.

Let us return to the approximation of the eigenvalues. Recall that the center
of gravity of the eigenvalues of $A$ depends on
$\overline{\beta}_{kk}=-\lambda_{k}|\overline{\varphi}_{k}|/[D_{k}(I\overline{\bm{\rho}})_{k}]$.
These values usually do not significantly change for changing indices $k$ over
a single refined pipe within a typical natural gas network. Moreover, the
values $\lambda_{k}/[D_{k}(I\bm{\rho})_{k}]$ do not change significantly for
changing indices $k$ over the network. Therefore, we expect to see two or more
imaginary asymptotes in the eigenvalues of $A$ if the network contains two or
more pipelines that deliver significantly different amounts of gas through
them. We extend the eigenvalue approximation from a single pipe to a connected
network of pipelines by approximating the eigenvalues for each pipeline in the
network separately using the average value of $\overline{\beta}$ corresponding
to each pipe and then collecting each of these subsets of approximate
eigenvalues to form the total set of approximate eigenvalues. We demonstrate
this approach for two network structures.

Figure 3: Eigenvalues ($\ast$) of $A$ and poles ($\circ$) of $G(s)$ for (left)
a single pipeline with three withdrawal nodes and (right) the cyclic network
from Section VII.

First, we consider the same pipeline in Figure 2 with supply at the inlet and
withdrawal at the outlet but with two additional withdrawal nodes located with
equidistant spacing. The network graph of this pipeline is actually a cascade
connection of three pipes of equal length. The withdrawal node closest to the
supply node withdraws 500 (kg/m2s), the subsequent node withdraws 200
(kg/m2s), and the outlet node withdraws 70 (kg/m2s). From (5), the initial
steady-state mass flux of gas in the first part of the pipeline is 770, that
in the second part is 270, and that in the last part is 70. Since these mass
flux values are significantly different from one another, we expect to see
three distinct imaginary asymptotes in the eigenvalues of $A$. The eigenvalues
of $A$ and combined poles of the three individual pipe segments are depicted
on the left side of Figure 3. The right side of the figure depicts the
eigenvalues and poles of the cyclic network in Figure 4 that we study later on
in Section VII. This cyclic network contains five pipelines. Three of the
pipelines have nearly equal center of gravity values and this is evident by
the heavily-weighted imaginary asymptote that intersects the real axis of the
complex plane at apprroximately $s=-0.025$. One imaginary asymptote is more
isolated at approximately $s=-0.040$. Lastly, there is another imaginary
asymptote at approximately $s=-0.034$, but the pipeline corresponding to this
asymptote contains only one internal refined node.

##  VI Digital Implementation

The nonlinear and linearized ODE-constrained OCPs in (21) and (29) may be
generally expressed as

$\displaystyle\min_{x(t),u(t)}\qquad$
$\displaystyle\int_{0}^{T}\mathcal{F}(x(t),u(t))dt$ (32a) s.t.
$\displaystyle\dot{x}(t)=f(x(t),u(t),d(t)),$ (32b) $\displaystyle x(0)=0,$
(32c) $\displaystyle\text{l.b.}\leq(x(t),u(t))\leq\text{u.b.},$ (32d)

where $x=(\bm{\rho},\varphi)^{T}$ is the state, $u=\mu$ is the input, and
$d=(\bm{s},\bm{w})^{T}$ is a known vector of disturbances. The symbols l.b.
and u.b. represent vector bounds that constrain the components of the state
and input vectors, where the inequalities are applied componentwise. Nonlinear
and linear programs are obtained by discretizing the time interval $[0,T]$
into $N$ subintervals $(t_{n},t_{n+1})$, for $n=0,\dots,(N-1)$, where the
sampling times are defined by $t_{n}=nT/N$.

The vector-valued functions $x(t)$, $u(t)$, and $d(t)$ are sampled to form
finite sequences of discrete vector-values $x[n]=x(t_{n})$, $u[n]=u(t_{n})$,
and $d[n]=d(t_{n})$, for $n=0,\dots,N$. The integral in the objective function
is approximated using the left-endpoint integration method, resulting in

$\int_{0}^{T}\mathcal{F}(x(t),u(t))dt\approx\sum_{n=0}^{N-1}\frac{T}{N}\mathcal{F}(x[n],u[n]).$

The time derivative of $x(t)$, for $t\in(t_{n},t_{n+1})$, is approximated with
Euler’s method

$\dot{x}(t)\approx\frac{x[n+1]-x[n]}{T/N}.$

The nonlinear and linear optimization programs are then defined by

$\displaystyle\min_{x[n],u[n]}\quad$
$\displaystyle\sum_{n=0}^{N-1}\frac{T}{N}\mathcal{F}(x[n],u[n])$ (33) s.t.
$\displaystyle\frac{x[n+1]-x[n]}{T/N}=f(x[n],u[n],d[n]),$ (34) $\displaystyle
x[0]=0,$ (35) $\displaystyle\text{l.b.}\leq(x[n],u[n])\leq\text{u.b.},$ (36)

where (34) is defined for $n=0,\dots,(N-1)$ and (36) is defined for
$n=1,\dots,N$. These linear and nonlinear optimization problems are performed
on a HP Spectre x360 4-core CPU with 16GB of unified memory, and are
implemented in Matlab with the sequential quadratic program algorithm using
the function fmincon. The gradient of the objective function is supplied for
an improvement in the performance.

##  VII Network Illustration

Figure 4: Top: Network configuration (not to scale). The triangles represent
compressor stations. Edge lengths in kilometer units: $\ell_{1}=20$,
$\ell_{2}=70$, $\ell_{3}=10$, $\ell_{4}=80$, $\ell_{5}=60$. The pipelines have
uniform diameter (0.9144 m) and friction factor (0.01), except for edge
$5\in\mathcal{E}$ that has diameter $D_{5}=0.635$ (m) and friction factor
$\lambda_{5}=0.015$. The speed of sound is $\sigma=377$ (m/s). Bottom: Mass
outflow boundary condition profiles at the color-coordinated nodes. Figure 5:
Optimal solution in pressure, mass flow at the inlets of the edges, and
compressor actions. Pressure and compression are color-coordinated with the
nodes and compressor stations of the network in Figure 4. Solid lines
represent the nonlinear solution and marker symbols represent the linearized
solution.

The solutions of the linear and nonlinear optimization problems are examined
for the test network and boundary conditions depicted in Figure 4. This
network was used in a previous study to validate a staggered grid
discretization method for simulation [24]. The blue node is the only supply
node. The pressure of gas entering this node is specified to be 5 (MPa) by
writing $\bm{s}_{\text{blue}}=5\times 10^{6}/\sigma^{2}$. The five edges of
the network are discretized into 48 refined edges, so that $\ell_{k}=5$ (km)
for all $k\in\hat{\mathcal{E}}$. The size of the state matrix $A$ is $95\times
95$ and the size of the input matrix $B$ is $95\times 3$. The time interval
$[0,T]$ is discretized into 24 subintervals using 25 evenly-spaced time
samples, as described in the previous section. This results in 2,450
optimization variables for both the nonlinear and linearized programs.

The optimal solutions of the nonlinear and linearized programs are shown in
Figure 5, where the solutions are translated back to non-variation variables
by adding their associated steady-state components. Figure 5 shows that the
linear program performs reasonably well at deciding the optimal solution for
these boundary conditions. The percent relative error between the nonlinear
and linearized solutions is $1.54\%$ in pressure, $0.57\%$ in mass flow, and
$2.3\%$ in compression, where the percent relative error of pressure between
the nonlinear solution $\bm{p}_{\text{non}}$ and the linearized solution
$\bm{p}_{\text{lin}}$ is defined by

$\text{Error in
}\bm{p}=\max_{j\in\mathcal{V}_{w},t\in[0,T]}\left(\frac{|(\bm{p}_{\text{non}})_{j}(t)-(\bm{p}_{\text{lin}})_{j}(t)|}{\overline{\bm{p}}_{j}}\right)\times
100.$ (37)

Here, pressure is defined by $\bm{p}=\sigma^{2}\bm{\rho}$ with appropriate
symbols for steady-state, linear, and nonlinear solutions. Errors in mass flow
and compression are defined similarly.

##  VIII Error Analysis

If $\bm{w}$, $\bm{s}$, and $\mu$ are time-invariant, then the initial-value
systems in (22)-(23) and (25) are undisturbed and unforced, in the sense that
the solutions of these systems remain at the origin for all time. Assuming
that the eigenvalues of $A$ have negative real parts, the origin of the linear
and nonlinear systems is an asymptotically and a locally asymptotically stable
equilibrium point, respectively. In this section, we assume that $\bm{s}=0$
and $\mu=0$, which is common for gas pipelines in practice.

Figure 6: Percent relative error (as in (37)) as a function of percent flow
variation $(100\kappa)$ as defined in (39) for the single pipeline studied in
Figure 2.

Proposition 2. Suppose that the real parts of the eigenvalues of $A$ are
negative, that $\mu=0$, $\bm{s}=0$, and $\overline{\varphi}_{k}>0$ for all
$k\in\mathcal{E}$. If the solution of (22)-(23) satisfies
$|\bm{\rho}_{j}|\leq\kappa\overline{\bm{\rho}}_{j}$ for all
$j\in\mathcal{V}_{w}$ and $|\varphi_{k}|\leq\kappa\overline{\varphi}_{k}$ for
all $k\in\mathcal{E}$ with $\kappa\in(0,\kappa_{\max})$ and $\kappa_{\max}<1$
sufficiently small, then there exists positive constants $a$, $b$, and $r$
such that

$\|\bm{e}\|\leq r\left(\frac{\kappa^{3}}{(1-\kappa)^{2}}a+\kappa^{2}b\right),$
(38)

where
$\bm{e}=[(\bm{\rho},\varphi)-(\bm{\rho}_{\text{lin}},\varphi_{\text{lin}})]$
and $\|\bm{e}\|$ is the Euclidean norm of the vector $\bm{e}$. Here
$(\bm{\rho},\varphi)$ is the solution of (22)-(23) and
$(\bm{\rho}_{\text{lin}},\varphi_{\text{lin}})$ is the solution of (25).

Proof. See the Appendix. $\square$

The bound in (38) provides a rate by which the norm of the error changes in
terms of the size of the variation around the steady-state. If $\kappa=0$,
then $\|\bm{e}\|=0$, as expected, since the solutions of the linear and
nonlinear systems are undisturbed and remain in equilibrium at the origin. As
$\kappa$ increases from $\kappa=0$ to $\kappa=\kappa_{\max}$, the norm of the
error increases. Intuitively, if the variation in density, $\bm{\rho}_{j}$,
approaches the negative of $\overline{\bm{\rho}}_{j}$ for some
$j\in\mathcal{V}_{w}$ (i.e., $\kappa\rightarrow 1$), then the magnitude of the
nonlinear term in (23) increases without bound. However, the solution of the
linear system is bounded due to the global asymptotic stability of the
equilibrium of the linear system. The bound of this extreme difference as
$\kappa\rightarrow 1$ is captured in the ratio on the right side of (38).
Although the bound in (38) provides some intuition about how the error changes
as a function of the size of the flow variation, it is derived using
conservative Lyapunov bounds, and, therefore, its quantitative values may be
impractical for real gas systems.

To improve the error estimate, we analyze the error in each flow variable of
the linear and nonlinear optimization problems as functions of the size of the
withdrawal variation. In particular, a series of optimization routines are
performed on the single pipeline that was introduced in Section V with
pipeline parameters given in Figure 2. The goal is to design the compressors
to minimize the expended energy for both the linear and nonlinear OCPs and
then compare each of the resulting optimal flow variables for a series of
different boundary conditions. The boundary conditions are constant 5 (MPa)
pressure at the inlet of the pipeline and varying mass outflow profiles at the
outlet given by

$\bm{w}_{\text{out}}(t)=200\left(1+\kappa\,\text{step}(t-T/4)-\kappa\,\text{step}(t-3T/4)\right),$
(39)

where $0\leq\kappa\leq 1$ and $\text{step}(t)=0.5\text{tanh}(7.2t)$ for
$t\in[0,T]$ with $T=24$ (hr). The percent relative error in pressure, mass
flux, and compression between the solutions of the nonlinear and linearized
optimization problems, as in (37), are depicted in Figure 6 as functions of
$\kappa\times 100$ (percent flow variation).

##  IX Conclusion

We have developed a linear control system of ODEs for natural gas flows in
pipeline networks. The eigenvalues of the linear state matrix were computed
numerically for a single pipe and compared to a subset of poles of an
irrational transfer matrix. The analysis was extended to networks of
pipelines, from which several conjectures can be made. First, the eigenvalues
of the state matrix are qualitatively similar to the collection of poles of
the transfer functions for some natural gas networks with typical pipeline
parameters. Second, the eigenvalues of the state matrix may have greater
dependence on the parameter values of each individual pipeline of the network
than the connectivity structure of the network. Third, the number of imaginary
asymptotes toward which the eigenvalues of the state matrix gravitate is at
most the number of edges of the original network graph. These conjectures are
based on the observations presented in Section V. The simple closed-form
representation of the poles of the transfer function was used to qualitatively
analyze the transient behavior of the natural gas pipeline system.

Control actions of compressor units were designed to minimize the energy that
the compressors expend. Although this OCP has received significant attention
in the literature, there has not been much work on comparing the linear and
nonlinear OCPs. The error between the solutions of the LP and NLP is analyzed
computationally as a function of the size of the variation in withdrawal flow.
It was demonstrated that the solution of the LP is within 5% relative error
from the solution of the NLP for sufficiently slow withdraw transients that
vary by less than 100% of the steady-state solution. This result generally
depends on the network structure, but the methodology may be applied by
control system operators using real data on sections of the network to
determine flow variation conditions under which the LP is sufficiently
accurate. We also derived an analytical bound on the error between the
solutions of the linear and nonlinear flow equations using Lyapunov functions.
One may use the LP for other applications as well. For example, if the error
is too large for certain flow variations near highly transient electric
generators, the control engineer may design adaptive switching between the LP
and NLP to reduce complexity and improve convergence of optimizing gas
dynamics for large network systems [35]. The results presented in this work
could open new capabilities to analyze very large pipeline network systems
that transport natural gas, or more generally, infrastructure networks
involving natural gas, water distribution, carbon capture and storage,
hydrogen blending, or traffic flow.

## References

  * [1] U.S. EIA, “US Energy Information Administration Annual Energy Outlook 2020,” _US Department of Energy: Washington, DC, USA_ , 2020.
  * [2] P. L. Wong and R. Larson, “Optimization of natural-gas pipeline systems via dynamic programming,” _IEEE Transactions on Automatic Control_ , vol. 13, no. 5, pp. 475–481, 1968.
  * [3] P. B. Percell and M. J. Ryan, “Steady state optimization of gas pipeline network operation,” in _PSIG annual meeting_. OnePetro, 1987.
  * [4] S. Wu, R. Z. Rios-Mercado, E. A. Boyd, and L. R. Scott, “Model relaxations for the fuel cost minimization of steady-state gas pipeline networks,” _Mathematical and Computer Modelling_ , vol. 31, no. 2-3, pp. 197–220, 2000\.
  * [5] D. De Wolf and Y. Smeers, “The gas transmission problem solved by an extension of the simplex algorithm,” _Management Science_ , vol. 46, no. 11, pp. 1454–1465, 2000.
  * [6] S. Misra, M. W. Fisher, S. Backhaus, R. Bent, M. Chertkov, and F. Pan, “Optimal compression in natural gas networks: A geometric programming approach,” _Transactions on Control of Network Systems_ , vol. 2, no. 1, pp. 47–56, 2014.
  * [7] M. Abbaspour and K. S. Chapman, “Nonisothermal Transient Flow in Natural Gas Pipeline,” _Journal of Applied Mechanics_ , vol. 75, no. 3, 2008.
  * [8] A. Osiadacz, “Simulation of transient gas flows in networks,” _International Journal for Numerical Methods in Fluids_ , vol. 4, no. 1, pp. 13–24, 1984.
  * [9] A. Zlotnik, M. Chertkov, and S. Backhaus, “Optimal control of transient flow in natural gas networks,” in _Conference on Decision and Control_. IEEE, 2015, pp. 4563–4570.
  * [10] H. H. Rachford and R. G. Carter, “Optimizing pipeline control in transient gas flow,” in _PSIG annual meeting_. OnePetro, 2000.
  * [11] A. Zlotnik, K. Sundar, A. M. Rudkevich, A. Beylin, and X. Li, “Optimal control for scheduling and pricing intra-day natural gas transport on pipeline networks,” in _Conference on Decision and Control_. IEEE, 2019, pp. 4887–4884.
  * [12] A. J. Osiadacz and M. Chaczykowski, “Comparison of isothermal and non-isothermal pipeline gas flow models,” _Chemical Engineering Journal_ , vol. 81, no. 1-3, pp. 41–51, 2001.
  * [13] A. R. D. Thorley and C. H. Tiley, “Unsteady and transient flow of compressible fluids in pipelines—a review of theoretical and some experimental studies,” _International journal of heat and fluid flow_ , vol. 8, no. 1, pp. 3–15, 1987.
  * [14] R. Z. Ríos-Mercado _et al._ , “A reduction technique for natural gas transmission network optimization problems,” _Annals of Operations Research_ , vol. 117, no. 1, pp. 217–234, 2002.
  * [15] S. Grundel, N. Hornung, B. Klaassen, P. Benner, and T. Clees, “Computing surrogates for gas network simulation using model order reduction,” in _Surrogate-Based Modeling and Optimization_. Springer, 2013, pp. 189–212.
  * [16] M. Behbahani-Nejad and Y. Shekari, “The accuracy and efficiency of a reduced-order model for transient flow analysis in gas pipelines,” _Journal of Petroleum Science and Engineering_ , vol. 73, no. 1-2, pp. 13–19, 2010.
  * [17] A. Beylin, A. M. Rudkevich, and A. Zlotnik, “Fast transient optimization of gas pipelines by analytic transformation to linear programs,” in _PSIG annual meeting_. OnePetro, 2020.
  * [18] P. Domschke, B. Geißler, O. Kolb, J. Lang, A. Martin, and A. Morsi, “Combination of nonlinear and linear optimization of transient gas networks,” _INFORMS Journal on Computing_ , vol. 23, no. 4, pp. 605–617, 2011.
  * [19] M. Sirvent, N. Kanelakis, B. Geißler, and P. Biskas, “Linearized model for optimization of coupled electricity and natural gas systems,” _Journal of Modern Power Systems and Clean Energy_ , vol. 5, no. 3, pp. 364–374, 2017.
  * [20] C. Himpe, S. Grundel, and P. Benner, “Model order reduction for gas and energy networks,” _Journal of Mathematics in Industry_ , vol. 11, no. 1, pp. 1–46, 2021.
  * [21] A. Zlotnik, S. Dyachenko, S. Backhaus, and M. Chertkov, “Model reduction and optimization of natural gas pipeline dynamics,” in _Dynamic Systems and Control Conference_ , vol. 57267. American Society of Mechanical Engineers, 2015.
  * [22] S. K. K. Hari, K. Sundar, S. Srinivasan, A. Zlotnik, and R. Bent, “Operation of natural gas pipeline networks with storage under transient flow conditions,” _IEEE Transactions on Control Systems Technology_ , vol. 30, no. 2, pp. 667–679, 2021.
  * [23] A. Zlotnik, L. Roald, S. Backhaus, M. Chertkov, and G. Andersson, “Coordinated scheduling for interdependent electric power and natural gas infrastructures,” _IEEE Transactions on Power Systems_ , vol. 32, no. 1, pp. 600–610, 2016.
  * [24] V. Gyrya and A. Zlotnik, “An explicit staggered-grid method for numerical simulation of large-scale natural gas pipeline networks,” _Applied Mathematical Modelling_ , vol. 65, pp. 34–51, 2019.
  * [25] A. J. Osiadacz and M. Gburzyńska, “Selected mathematical models describing flow in gas pipelines,” _Energies_ , vol. 15, no. 2, 2022.
  * [26] K. Sundar and A. Zlotnik, “State and parameter estimation for natural gas pipeline networks using transient state data,” _IEEE Transactions on Control Systems Technology_ , vol. 27, no. 5, pp. 2110–2124, 2019.
  * [27] S. Misra, M. Vuffray, and A. Zlotnik, “Monotonicity properties of physical network flows and application to robust optimal allocation,” _Proceedings of the IEEE_ , vol. 108, no. 9, pp. 1558–1579, 2020.
  * [28] M. Gugat, M. Herty, A. Klar, G. Leugering, and V. Schleper, “Well-posedness of networked hyperbolic systems of balance laws,” _Constrained optimization and optimal control for partial differential equations_ , vol. 160, pp. 123–146, 2012.
  * [29] I. Marić, A. Galovic, and T. Smuc, “Calculation of natural gas isentropic exponent,” _Flow Measurement and Instrumentation_ , vol. 16, pp. 13–20, 03 2005.
  * [30] E. Feldhoff, T. Duphorn, S. Schlegel, and D. Westermann, “Cross-sectoral provision of ancillary services in coupled electricity and gas sectors on a real-time simulator,” in _2022 IEEE PES Innovative Smart Grid Technologies Conference Europe (ISGT-Europe)_. IEEE, 2022, pp. 1–5.
  * [31] E. Feldhoff, S. Schlegel, D. Westermann, S. Naumann, and P. Bretschneider, “Co 2-and cost-optimal multi-energy system modelling considering the provision of operating reserve,” in _2022 18th International Conference on the European Energy Market (EEM)_. IEEE, 2022, pp. 1–6.
  * [32] L. Baker, D. Armbruster, A. Scaglione, and R. B. Platte, “Analysis of a model of a natural gas pipeline—a transfer function approach,” _Transactions of Mathematics and Its Applications_ , vol. 5, no. 1, 11 2021\.
  * [33] A. C. Zecchin, A. R. Simpson, M. F. Lambert, L. B. White, and J. P. Vítkovskỳ, “Transient modeling of arbitrary pipe networks by a laplace-domain admittance matrix,” _Journal of Engineering Mechanics_ , vol. 135, no. 6, pp. 538–547, 2009.
  * [34] J. Králik, P. Stiegler, Z. Vostrú, and J. Závorka, “Modeling the dynamics of flow in gas pipelines,” _IEEE transactions on systems, man, and cybernetics_ , vol. 14, no. 4, pp. 586–596, 1984.
  * [35] P. Bales, O. Kolb, and J. Lang, “Hierarchical modelling and model adaptivity for gas flow on networks.” in _ICCS (1)_ , 2009, pp. 337–346.
  * [36] H. K. Khalil, _Nonlinear control_. Pearson New York, 2015, vol. 406.

##  X Appendix

Proof of Proposition 2. The proof is based on standard Lyapunov arguments and
we provide a brief outline here (see, e.g. [36] for more examples). Assume
that $(\bm{\rho},\varphi)$ satisfies (22)-(23) and that
$(\bm{\rho}_{\text{lin}},\varphi_{\text{lin}})$ satisfies (25) under the
conditions stated in Proposition 2. Define the error vector
$\bm{e}=[(\bm{\rho},\varphi)-(\bm{\rho}_{\text{lin}},\varphi_{\text{lin}})]$.

Since $\mu=0$, it follows that $M=0$ and $N=0$. The function
$f(\bm{\rho},\varphi,\bm{w})$ on the right-hand side of (22)-(23) may be
uniquely written as
$f(\bm{\rho},\varphi,\bm{w})=F(\bm{\rho},\varphi)+\tilde{W}\bm{w}$, where
$\tilde{W}$ is defined to be the second block column of $W$ as

$\tilde{W}=\begin{bmatrix}-(I^{T}XLI)^{-1}\\\ 0\end{bmatrix}.$

The error satisfies

$\displaystyle\dot{\bm{e}}$
$\displaystyle=\left(\begin{bmatrix}\dot{\bm{\rho}}\\\
\dot{\varphi}\end{bmatrix}-\begin{bmatrix}\dot{\bm{\rho}}_{\text{lin}}\\\
\dot{\varphi}_{\text{lin}}\end{bmatrix}\right)$
$\displaystyle=F(\bm{\rho},\varphi)+\tilde{W}\bm{w}-A\begin{bmatrix}\bm{\rho}_{\text{lin}}\\\
\varphi_{\text{lin}}\end{bmatrix}-\tilde{W}\bm{w}$
$\displaystyle=A\bm{e}+F(\bm{\rho},\varphi)-A\begin{bmatrix}\bm{\rho}\\\
\varphi\end{bmatrix}$ $\displaystyle=A\bm{e}+\left(\partial
F(\bm{\rho}^{*},\varphi^{*})-A\right)\begin{bmatrix}\bm{\rho}\\\
\varphi\end{bmatrix},$ (40)

where we have used the mean-value theorem to replace $F$ by

$F(\bm{\rho},\varphi)=F(0,0)+\partial
F(\bm{\rho}^{*},\varphi^{*})\begin{bmatrix}\bm{\rho}\\\ \varphi\end{bmatrix},$

for some point $(\bm{\rho}^{*},\varphi^{*})$ on the line segment connecting
$(\bm{\rho},\varphi)$ to the origin. In the above equation, $F(0,0)=0$ and the
symbol $\partial F(\bm{\rho}^{*},\varphi^{*})$ represents the Jacobian matrix
of the vector-valued function $F$ evaluated at the point
$(\bm{\rho}^{*},\varphi^{*})$. Therefore, the dynamics of $\bm{e}$ are
represented as a perturbed linear system with a nonvanishing perturbation
term.

The difference between the Jacobians in (40) is given by

$\displaystyle\partial F(\bm{\rho}^{*},\varphi^{*})-A=\begin{bmatrix}0&0\\\
(\alpha^{*}-\overline{\alpha})I&\beta^{*}-\overline{\beta}\end{bmatrix},$

where $\overline{\alpha}$ and $\overline{\beta}$ are defined in (24) and

$\displaystyle\alpha^{*}$ $\displaystyle=$ $\displaystyle
K\text{diag}\left(\frac{\left(\varphi^{*}+\overline{\varphi}\right)\odot\left|\varphi^{*}+\overline{\varphi}\right|}{I\left(\bm{\rho}^{*}+\overline{\bm{\rho}}\right)^{2}}\right),$
(41) $\displaystyle\beta^{*}$ $\displaystyle=$
$\displaystyle-2K\text{diag}\left(\frac{\left|\varphi^{*}+\overline{\varphi}\right|}{I\left(\bm{\rho}^{*}+\overline{\bm{\rho}}\right)}\right).$
(42)

The matrices $K$, $\overline{\alpha}$, $\alpha^{*}$, $\overline{\beta}$, and
$\beta^{*}$ are all diagonal. The $k$-th diagonal entries of the matrices
$(\alpha^{*}-\overline{\alpha})$ and $(\beta^{*}-\overline{\beta})$ with
$k:i\mapsto j$ are given by

$\displaystyle(\alpha_{kk}^{*}-\overline{\alpha}_{kk})$ $\displaystyle=$
$\displaystyle
K_{kk}\left(\frac{(\varphi_{k}^{*}+\overline{\varphi}_{k})|\varphi_{k}^{*}+\overline{\varphi}_{k}|}{\left(\bm{\rho}_{j}^{*}+\overline{\bm{\rho}}_{j}\right)^{2}}-\frac{\overline{\varphi}_{k}|\overline{\varphi}_{k}|}{\overline{\bm{\rho}}_{j}^{2}}\right),$
$\displaystyle(\beta_{kk}^{*}-\overline{\beta}_{kk})$ $\displaystyle=$
$\displaystyle
K_{kk}\left(-\frac{2|\varphi_{k}^{*}+\overline{\varphi}_{k}|}{\bm{\rho}_{j}^{*}+\overline{\bm{\rho}}_{j}}+\frac{2|\overline{\varphi}_{k}|}{\overline{\bm{\rho}}_{j}}\right).$

Since $|\bm{\rho}_{j}|\leq\kappa\overline{\bm{\rho}}_{j}$ and
$|\varphi_{k}|\leq\kappa\overline{\varphi}_{k}$, it can be shown with
elementary inequalities that

$\displaystyle\|\alpha^{*}-\overline{\alpha}\|\leq\frac{2\kappa^{2}}{(1-\kappa)^{2}}\|\overline{\alpha}\|,\qquad\|\beta^{*}-\overline{\beta}\|\leq
2\kappa\|\overline{\beta}\|,$

where the norm of a matrix is defined to be the induced Euclidean norm. It
then follows from the inequality
$\|(\bm{\rho},\varphi)\|\leq\kappa(\|\overline{\bm{\rho}}\|+\|\overline{\varphi}\|)$
that

$\displaystyle\|(\partial
F(\bm{\rho}^{*},\varphi^{*})-A)(\bm{\rho},\varphi)^{T}\|$
$\displaystyle\qquad\leq\kappa\left(\|\alpha^{*}-\overline{\alpha}\|\|I\|+\|\beta^{*}-\overline{\beta}\|\right)(\|\overline{\bm{\rho}}\|+\|\overline{\varphi}\|)$
$\displaystyle\qquad\leq\frac{\kappa^{3}}{(1-\kappa)^{2}}a+\kappa^{2}b,$

where
$a=2\|\overline{\alpha}\|\|I\|(\|\overline{\bm{\rho}}\|+\|\overline{\varphi}\|)$
and
$b=2\|\overline{\beta}\|(\|\overline{\bm{\rho}}\|+\|\overline{\varphi}\|)$.

Since the eigenvalues of $A$ have negative real parts, we may define the
Lyapunov function $V(\bm{e})=\bm{e}^{T}R\bm{e}$, where $R$ is a positive
definite and symmetric matrix that satisfies the Lyapunov equation
$(A^{T}R+RA)=-I$. The time derivative of $V(\bm{e})$ along the trajectories of
(40) is bounded by $\dot{V}(\bm{e})\leq(-\|\bm{e}\|^{2}+L\|\bm{e}\|)$, where
$L=2\|R\|(a\kappa^{3}/(1-\kappa)^{2}+b\kappa^{2})$. Define
$Z(\bm{e})=\sqrt{V(\bm{e})}$. Since $R$ is symmetric, it has real eigenvalues.
It follows from Rayleigh’s quotient that $\lambda_{1}\|\bm{e}\|^{2}\leq
V(\bm{e})\leq\lambda_{2}\|\bm{e}\|^{2}$, where $\lambda_{1}$ and $\lambda_{2}$
are the minimum and maximum eigenvalues of $R$. For $V(\bm{e})>0$,
$\dot{Z}(\bm{e})\leq(-\frac{1}{2\lambda_{2}}Z(\bm{e})+\frac{L}{2\sqrt{\lambda_{1}}})$.
By the comparison lemma,

$\displaystyle
Z(t)\leq\frac{L}{2\sqrt{\lambda_{1}}}\int_{0}^{t}e^{-(t-\tau)/2\lambda_{2}}d\tau\leq\frac{\lambda_{2}L}{\sqrt{\lambda_{1}}}.$

It can be shown that the above inequality also holds for $V=0$ [36]. The
proposition then follows from
$\|\bm{e}\|\leq\sqrt{V/\lambda_{1}}\leq\lambda_{2}L/\lambda_{1}$ with
$r=2\|R\|\lambda_{2}/\lambda_{1}$. Throughout this proof, we assume that
$\kappa\leq\kappa_{\text{max}}$ to ensure that the nonlinear solution is well-
defined.