Iterative LP-based Methods for the Multiperiod Optimal Electricity and Gas Flow Problem

Sleiman Mhanna, Member, IEEE, Isam Saedi, Student Member, IEEE,
Pierluigi Mancarella, Senior Member, IEEE

In light of the increasing coupling between electricity and gas networks, this paper introduces two novel iterative methods for efficiently solving the multiperiod optimal electricity and gas flow (MOEGF) problem. The first is an iterative MILP-based method and the second is an iterative LP-based method with an elaborate procedure for ensuring an integral solution. The convergence of the two approaches is founded on two key features. The first is a penalty term with a single, automatically tuned, parameter for controlling the step size of the gas network iterates. The second is a sequence of supporting hyperplanes together with an increasing number of carefully constructed halfspaces for controlling the convergence of the electricity network iterates. Moreover, the two proposed algorithms use as a warm start the solution from a novel polyhedral relaxation of the MOEGF problem, for a noticeable improvement in computation time as compared to a cold start. Unlike the first method, which invokes a branch-and-bound algorithm to find an integral solution, the second method implements an elaborate steering procedure that guides the continuous variables to take integral values at the solution. Numerical evaluation demonstrates that the two proposed methods can converge to high-quality feasible solutions in computation times at least two orders of magnitude faster than both a state-of-the-art nonlinear branch-and-bound (NLBB) MINLP solver and a mixed-integer convex programming (MICP) relaxation of the MOEGF problem. The experimental setup consists of five test cases, three of which involve the real electricity and gas transmission networks of the state of Victoria with actual linepack and demand profiles.

Sequential linear programming, polyhedral envelopes, integrated electricity and gas systems, piecewise linear approximations, mixed-integer nonlinear programming, mixed-integer second-order cone programming, branch-and-bound.

§ NOTATION

§.§ Sets

[$\phi_{mn}^{{\rm in/o},t}$]

$\mathcal{B}$ Set of buses in the electricity network.

$\mathcal{B}_{i}$ Set of buses adjacent to bus $i$.

$\mathcal{C}$ Set of compressors in the gas network.

$\mathcal{D}^{\rm g}$ Set of gas-powered generators (GPGs) on the gas network.

$\mathcal{D}_{m}^{\rm g}$ Set of GPGs connected to node $m$.

$\mathcal{E}$ Set of non-pipe elements.

$\mathcal{G}^{\rm g/ng}$ Set of GPGs/non-GPGs in the electricity network.

$\mathcal{G}_{i}^{\rm g/ng}$ Set of GPGs/non-GPGs connected to bus $i$.

$\mathcal{L}/\mathcal{L}^{\rm t}$ Set of all branches $ij$/$ji$ where $i$/$j$ is the “from”/“to” bus.

$\mathcal{N}$ Set of nodes in the gas network.

$\mathcal{P}$ Set of pipelines in the gas network.

$\mathcal{R}$ Set of pressure regulators in the gas network.

$\mathcal{S}$ Set of gas supplies in the gas network.

$\mathcal{S}_{m}$ Set of gas supplies connected to node $m$.

$\mathcal{T}$ Decision time horizon.

§.§ Parameters and input data

[$\phi_{mn}^{{\rm in/o},t}$]

$c_{0,gi}$ Constant coefficient ($\SI{}{\$\per\hour}$) term of GPG/non-GPG $g$'s cost function.

$c_{1,gi}$ Coefficient ($\SI{}{\$\per\mega\watt\hour}$) of the linear term of GPG/non-GPG $g$'s cost function.

$c_{2,gi}$ Coefficient ($\SI{}{\$\per\mega\watt\hour\squared}$) of the quadratic term of GPG/non-GPG $g$'s cost function.

$c_{sm}$ Cost ($\SI{}{\$\per\meter\cubed}$) of gas production of gas supply $s$ at node $m$.

$\Delta \tau,\Delta t$ Time resolution ($\SI{}{\second},\SI{}{\hour}$).

$D_{mn}$ Diameter ($\SI{}{\meter}$) of pipeline $mn$.

$g^{\rm sh}_{i}$ Shunt conductance (pu) at bus $i$.

$\mathrm{i}$ Imaginary unit.

$k$ Iteration number.

$\eta_{gi}$ Efficiency of GPG $g$ at bus $i$.

$\phi_{m}^{{\rm d},t}$ Gas demand ($\SI{}{\meter\cubed\per\second}$) at node $m$.

$L_{mn}$ Length ($\SI{}{\meter}$) of pipeline $mn$.

$HHV$ Higher heating value ($\SI{38.07}{\mega\joule\per\meter\cubed}$) of natural gas.

$\nu_{mn}$ Gas consumption coefficient of compressor $mn$.

$p_{i}^{{\rm d},t}$ Active power demand ($\SI{}{\mega\watt}$) at bus $i$.

$p_{gi}^{\rm RD/RU}$ Ramp down/up rate ($\SI{}{\mega\watt}/\SI{}{\hour}$) of GPG/non-GPG $g$ at bus $i$.

$R$ Specific gas constant ($\SI{478.42}{\joule\per\kilogram\per\kelvin}$) at standard conditions.

$\rho$ Gas density ($\SI{0.735}{\kilogram\per\meter\cubed}$) at standard conditions.

$\overline{\gamma}_{mn}/\underline{\gamma}_{mn}$ Upper/Lower limit on the compression (pressure drop) ratio of a compressor (pressure regulator).

$S$ Specific gravity (relative density) of natural gas, $S = 0.6$ (dimensionless).

$T$ Gas temperature ($\SI{288.15}{\kelvin}$) at standard conditions.

$T_{ij}$ Complex tap ratio of a phase shifting transformer ($T_{ij}=\tau_{ij}\mathrm{e}^{\mathrm{i} \theta_{ij}^{\rm shift}}$).

$Y_{ij}$ Series admittance (pu) in the $\pi$-model of branch $ij$.

$Z_{mn}$ Compressibility factor (dimensionless) of gas in pipeline $mn$.

§.§ Operators

[$\phi_{mn}^{{\rm in/o},t}$]

$\bullet^*$ Conjugate operator.

$\Im/\Re\left\{\bullet\right\}$ Imaginary/Real value operator.

$\underline{\bullet}/\overline{\bullet}$ Minimum/Maximum magnitude operator.

$\left|\bullet\right|$ Magnitude operator/Cardinality of a set.

$\nint{\bullet}$ Round-to-nearest-integer operator.

$\times$ Cross product of two vectors.

§.§ Variables

[$\phi_{mn}^{{\rm in/o},t}$]

$p_{gi}^{t}$ Active power ($\SI{}{\mega\watt}$) generation of GPG/non-GPG $g$ at bus $i$.

$p_{ij}^{t}$ Active power ($\SI{}{\mega\watt}$) flow along branch $ij$.

$\theta_{i}^{t}$ Voltage angle ($\SI{}{\radian}$) at bus $i$.

$\wp_{m}^{t}$ Gas pressure ($\SI{}{\pascal}$) at node $m$.

$\tilde{\wp}_{mn}^{t}$ Average gas pressure ($\SI{}{\pascal}$) across pipeline $mn$.

$\phi_{mn}^{t}$ Gas flow rate ($\SI{}{\meter\cubed\per\second}$) across edge $mn$.

$\phi_{mn}^{{\rm in/out},t}$ Inflow/Outflow rate of gas ($\SI{}{\meter\cubed\per\second}$) of pipeline $mn$.

$\phi_{mn}^{{\rm c},t}$ Gas trapped in compressor $mn$ when it is boosting pressure.

$\phi_{mn}^{+,t}$ Gas flowing across compressor $mn$ when it is operating in the direction of pressure boost.

$\phi_{mn}^{-,t}$ Gas flowing across compressor $mn$ when it is operating in a reverse direction.

$\phi_{gm}^{t}$ Gas consumption ($\SI{}{\meter\cubed\per\second}$) of GPG $g$ at node $m$.

$\phi_{sm}^{t}$ Gas flow rate ($\SI{}{\meter\cubed\per\second}$) from gas supply $s$ at node $m$.

$\ell_{mn}^{t}$ Linepack ($\SI{}{\meter\cubed}$) of pipeline $mn$.

$z_{mn}^{t}$ Binary variables for non-pipe elements (also for pipelines in the MICP relaxation).

§ INTRODUCTION

Traditionally, energy systems like electricity and gas were modeled, operated, and managed separately. However, recent developments in power-to-gas (PtG) and co-generation technologies, and the heavy reliance on gas-powered generators (GPGs) to balance intermittent generation from renewable energy sources, have all prompted a paradigm shift towards jointly modeling and operating electricity and gas systems [1]. The main sources of complexity in the modeling of integrated electricity and gas systems (IEGS) are the alternating current (AC) power flow equations in the electricity transmission network, the equations describing the dynamic behavior of gas flow in pipelines, and the disjoint sets describing the operation of non-pipe elements such as compressors and pressure regulators in the gas transmission network. This modeling of IEGS results in an NP-hard MINLP problem that has a nonconvex continuous relaxation, which is extremely challenging to solve, even to local optimality, using current state-of-the-art MINLP technology.

Therefore, tractable alternatives such as linear programming (LP) approximations [2, 3, 4], second-order cone programming (SOCP) relaxations [5, 6, 7], semidefinite programming (SDP) relaxations [8], mixed-integer linear programming (MILP) approximations [9, 10, 11, 12, 13, 14, 15, 16], and mixed-integer second-order cone programming (MISOCP) relaxations [17, 18, 19, 20, 21, 22, 23, 24], are garnering considerable attention in the research community.

However, unless realistic approximations and assumptions are involved, achieving tractability may come at the price of infeasibility.

On the gas network side, one common oversimplification is to adopt a steady-state gas flow model instead of a dynamic gas flow model. In reality, the behavior of gas flow is characterized by much slower dynamics as compared to electric power flow, which means that gas demand and supply are not balanced instantaneously. These slow dynamics give rise to the linepack, which is instrumental in assessing gas network flexibility [25].[The linepack is the volume of gas that can be stored in a pipeline.] In contrast, steady-state gas flow assumes that supply and demand are balanced instantaneously, which does not account for the effect of the linepack, and therefore leads to an overestimation of the gas injections from gas suppliers [26].

Another common unrealistic assumption is related to the direction of flow in pipes and non-pipe elements. The assumption of known direction of gas flow in pipe and non-pipe elements is conducive for two main reasons. First, it obviates the disjoint sets describing the operation of non-pipe elements. Second, it renders the equation of gas flow in a pipe easier to convexify or approximate [5, 6, 7]. This assumption is valid under certain conditions where the demand predictably fluctuates within the day. However, in a multiperiod setting, this assumption would no longer be valid in the planning and operation of future IEGS under different clean-fuel (e.g., hydrogen, synthetic methane, etc.) injection and storage scenarios [27], and different electrification versus decarbonization scenarios of the heating and transportation sectors [28, 29]. Without any assumptions on the direction of flow in a pipe, there are two approaches to approximate or convexify the underlying nonconvex equality constraint. The first is through piecewise linear (PWL) techniques, which can be divided into two types. In the first type, each univariate term is approximated by PWL segments [9, 10, 12, 14, 15]. The second type consists of higher dimensional PWL techniques, such as the 2-D grid triangulation domains (in the 3-D Euclidean space) in [11], and the 3-D PWL techniques which are based on the Taylor series expansion of the 3-variable equation of gas flow in pipelines [13, 16]. The aim of these methods is to approximate the original MINLP problem by the more appealing MILP, as it is computationally more efficient and provides a measure of optimality. However, even for a modest approximation accuracy, these MILP approximations would require a computationally prohibitive number of PWL segments [30].[This claim is also true for PWL methods employing special ordered sets of type 2 (SOS2). SOS2 constraints, which are ordered sets of variables where at most two variables in the set may take non-zero values, instruct the branch and bound algorithm to branch on sets of variables, rather than individual variables. Examples of the latter include PWL methods that use binary variables for each segment or a logarithmic number of binary variables. A numerical comparison of state-of-the-art PWL methods to the proposed algorithms in this paper is given in [31] as part of this work.] The second approach consists of relaxing the nonconvex constraint to an MISOCP one [17, 18, 20, 21, 22]. However, these methods can still result in an infeasible solution as the original nonconvex equality constraint is relaxed into two disjoint SOC inequality constraints. In other words, the MISOCP approach is only feasible if one of the two SOC constraints is active at the optimum. One way to recover a feasible solution is to implement the convex-concave procedure [32], which in this setting entails solving a series of MISOCP problems [19, 23, 24], but now at the expense of more computation time.

On the electricity network side, common practice is to approximate the AC power flow model by a DC power flow model (see [33, 10, 34, 16, 35, 36, 4, 19, 18, 11, 12, 13, 14, 15, 37, 21, 26]), as it is linear and therefore amenable to convex programming, MILP, and MISOCP. The main drawback of the DC power flow model is its inability to capture losses in transmission lines (and transformers), which may result in an oversimplification of the problem. One way to include these losses is to approximate the cosine terms in the AC power flow constraints by their second-order Maclaurin series [2, 3, 7, 22].

Against this background, this paper proposes two iterative LP-based methods to directly solve the multiperiod optimal electricity and gas flow (MOEGF) problem with quasi-dynamic gas constraints,[In a quasi-dynamic gas flow model, the continuity equation and the motion equation are discretized over the full length of the pipe, whereby the input flow and the output flow are different. In contrast, a steady-state model assumes constant flow across a pipe, i.e., the input and output flows are equal. As a consequence, the steady-state model fails to capture the linepack, which results in unrealistic gas supply profiles and pressures across the network. Interested readers are referred to [38] and [39] for more details on the derivation of the quasi-dynamic gas flow model.] strengthened DC power flow, bidirectional pipes, compressors, and pressure regulators. Without the procedures that address the integrality constraints, the two methods are similar in concept to sequential linear programming (SLP), which has been applied to the optimal power flow (OPF) problem in [40], and the gas transmission problem in [41], but has not been applied to the MOEGF problem with quasi-dynamic gas constraints and disjoint sets. The step size in [40] and [41] is controlled by tightening the bounds on selected variables at each iteration, which requires multiple parameters and problem-dependent tuning, and can result in infeasible LP subproblems. In contrast, in this work, the step size for the gas network iterates is controlled in the objective function by a penalty term with a single, automatically tuned, parameter. More interestingly, electricity network iterates are not controlled by a step size but by dynamically introducing affine and linear cuts as a sequence of supporting hyperplanes and an increasing number of carefully constructed supporting halfspaces. Relevant previous works involving LP consisted mainly of first-order Taylor series approximations around a fixed pre-determined operating point [2, 3, 4]. However, although computationally efficient, those works do not yield feasible solutions.

In a nutshell, this paper advances the state of the art in the following ways:

* It presents a novel iterative MILP-based method that is numerically demonstrated to converge to near-optimal and feasible solutions in at least one order of magnitude faster than both a state-of-the-art NLBB MINLP solver and a mixed-integer convex programming (MICP) relaxation of the problem.

* It presents a novel fast iterative LP-based method that is numerically demonstrated to converge to high-quality feasible solutions in at least two orders of magnitude faster than a state-of-the-art NLBB MINLP solver. Unlike the iterative MILP-based method which invokes a branch-and-bound algorithm to find an integral solution, this approach implements a fast iterative LP-based steering procedure that guides the continuous variables to take integral values at the solution.

* It presents novel polyhedral envelopes for the direction of flow and average pressure in a pipeline and proves that the nonconvex equality constraint describing the average pressure in a pipeline becomes convex when relaxed into an inequality constraint.

* It numerically shows that the two proposed methods are robust against the choice of starting point but warm-starting them from the solution of the proposed polyhedral relaxation substantially improves computational speed. The proposed polyhedral relaxation also provides a valid starting point when no prior information is available.

While convex relaxations such as the SOCP, SDP, and the MISOCP are necessary to better understand the structure and complexity of the problem, they generally do not yield feasible solutions. Moreover, existing SOCP relaxations in [5, 6, 7] can only be applied when the direction of flow in pipeline is known beforehand, and the existing SDP relaxation in [8] can only be applied in steady-state modeling where the square of the pressure terms can be substituted by linear ones. On the other hand, the MISOCP relaxation in the transient domain in [17, 18, 19, 20, 21, 22], although able to model bi-directional gas flow in pipes, also yields infeasible solutions in general. This, along with the heavy computational burden of the MISOCP relaxation, are numerically confirmed in Section <ref> of this paper on real-world systems. Nonetheless, the MICP-relaxed problem provides a good-quality lower bound than can be used to assess the optimality of the two proposed methods. In contrast, the two proposed method are not convex relaxations, but are rather founded on a computationally efficient SLP method that exploits the inherent engineering limits dictated by the design parameters of components, as well as the structure and other properties of the MOEGF problem with quasi-dynamic gas flow modeling, to obtain high-quality feasible solutions. The validity of both the modeling and the solutions are demonstrated on the real Victorian gas transmission system with actual linepack and pressure profiles from the Australian Energy Market Operator (AEMO) [42].

The purpose and innovation of the proposed two iterative algorithms are two fold. First, they leverage the superior computational efficiency of state-of-the-art MILP and LP solvers, while retaining the accuracy of interior-point methods (IPM). Second, since the market dispatch engines of most Independent System Operators (ISO) around the world use MILP or LP solvers to obtain the locational marginal prices (LMP), by solving a problem that incorporates an LP approximation of the OPF problem, this approach is designed to allow ISO to retain their MILP or LP solvers but now with the added benefits of capturing the couplings between electricity and gas networks, gas network flexibility, and LMP that reflect transmission line losses, all in an integrated optimization framework.

It is worth noting that current electricity and gas markets are still operated separately. However, extensive recent literature has shown that an integrated operation of both systems leads to notable operational cost savings and a better quantification of the flexibility of both systems [1, 38, 39, 43, 25]. In addition, the advent of power-to-gas technologies will further increase the interactions between the two networks, thereby reinforcing the need for an integrated modeling [27]. Nonetheless, the concepts underlying the proposed algorithms in this paper can still be straightforwardly applied to each market/system independently.

The paper is organized as follows. Section <ref> describes the MOEGF problem and Section <ref> presents an MICP relaxation and an MISOCP relaxation of the MOEGF problem. Section <ref> introduces the polyhedral envelopes of the nonconvex sets, and the resulting polyhedral relaxation of the MOEGF problem. The proposed iterative MILP-based method and the iterative LP-based method are introduced in Sections <ref> and <ref>, respectively. The optimality, feasbility, and computational effort of the two proposed methods are compared to those of a state-of-the-art NLBB solver and the MICP relaxation of the MOEGF problem in Section <ref>, which also assesses the output linepack of the two proposed methods against historical ones from AEMO. The paper concludes in Section <ref>.

§ MULTIPERIOD OPTIMAL ELECTRICITY AND GAS FLOW

The MOEGF problem consists of finding the least-cost dispatch of power from electric generators, and gas from gas supplies, to satisfy both electrical and gas demands at all buses in the electrical network and nodes in the gas network. Electric power flow is governed by physical laws such as Ohm's law and Kirchhoff's current law (KCL), and other technical restrictions. On the other hand, the flow of gas is governed by physical laws, such as the conservation of flow at each node, the quasi-dynamic behavior of gas flow in pipelines, and other operational requirements such as valve switching in compressor and pressure regulation stations. The MOEGF can be mathematically formulated as

\begin{align}
\hspace{-1.25cm} \underset{\substack{p_{gi}^{t},p_{ij}^{t},\theta_{i}^{t},\wp_{m}^{t}\\\tilde{\wp}_{mn}^{t},\phi_{mn}^{t},\ell_{mn}^{t},\phi_{mn}^{{\rm in},t}\\\phi_{nm}^{{\rm out},t},z_{mn}^{t},\phi_{gm}^{t},\phi_{sm}^{t}\\\phi_{nm}^{{\rm c},t},\phi_{nm}^{+,t},\phi_{nm}^{-,t}}}
{\mbox{ minimize}} & \sum_{t \in \mathcal{T}} \Bigg( \sum_{gi \in \mathcal{G}} \hspace{-0cm} f_{gi}\left(p_{gi}^{t}\right)\Bigg. + \hspace{-0.2cm} \hspace{-0cm} & \hspace{-0cm} & \nonumber
\end{align}

\begin{align}
& \hspace{0cm} \Bigg. \sum_{sm\in \mathcal{S}} f_{sm}\left(\phi_{sm}^{t}\right)\Bigg) \hspace{-3cm} & \hspace{-0cm} & \label{OEGF_objective} \\
\text{subject to } \ \ \underline{p}_{gi} \leq p_{gi}^{t} & \leq \overline{p}_{gi}, \hspace{-1.5cm} & \hspace{-0cm} gi \in \mathcal{G} \label{OEGF_Pminmax} \\
- p_{gi}^{\rm RD} \Delta t \le p_{gi}^{t} - p_{gi}^{t-1} & \leq p_{gi}^{\rm RU} \Delta t, \hspace{-1.5cm} & \ \ \quad \quad gi \in \mathcal{G} \label{OEGF_Pramp} \\
-\overline{p}_{ij} \leq p_{ij}^{t} & \leq \overline{p}_{ij}, \hspace{-1.5cm} & ij \in \mathcal{L} \cup \mathcal{L}&^{\rm t} \label{OEGF_MWlimits} \\
\theta_{ij}^{t} = \theta_{i}^{t}-\theta_{j}^{t}, \ \underline{\theta}_{ij} & \leq \theta_{ij}^{t} \leq \overline{\theta}_{ij}, \hspace{-1.5cm} & ij \in \mathcal{L} & \label{OEGF_anglediff} \\
\sum_{g \in \mathcal{G}_{i}^{\rm g} \cup \mathcal{G}_{i}^{\rm ng}} p_{gi}^{t} = p_{i}^{{\rm d},t} & + \sum_{j \in \mathcal{B}_{i}} p_{ij}^{t} + g^{\rm sh}_{i}, \hspace{-1.5cm} & i \in \mathcal{B} & \label{OEGF_KCLp} \\
p_{ij}^{t}= g_{ij} \hspace{0.5cm} & \hspace{-0.5cm} 0.5(\theta_{ij}^{t})^2 + b_{ij}\theta_{ij}^{t}, \hspace{-1.5cm} & ij \in \mathcal{L} & \label{OEGF_Pij} \\
p_{ji}^{t}= g_{ji} \hspace{0.5cm} & \hspace{-0.5cm} 0.5(\theta_{ij}^{t})^2 - b_{ji}\theta_{ij}^{t}, \hspace{-1.5cm} & ji \in \mathcal{L}&^{\rm t} \hspace{-0.1cm} \label{OEGF_Pji} \\
p_{gi}^{t}=\phi_{gm}^{t} HHV &\eta_{gi}, & \hspace{-3cm} gi \in \mathcal{G}^{\rm g}, \ gm \in \mathcal{D}&^{\rm g} \label{OEGF_gasflow_GFG}\\
\underline{\phi}_{sm}^{\rm s} \leq \phi_{sm}^{t} & \leq \overline{\phi}_{sm}^{\rm s}, \hspace{-0.5cm} & sm \in \mathcal{S}& \label{OEGF_Fminmax} \\
%			- \phi_{sm}^{\rm RD}\Delta t \le \phi_{sm}^{t} & - \phi_{sm}^{t-1} \leq \phi_{sm}^{\rm RU}\Delta t, \hspace{-1.5cm} & sm \in \mathcal{S}& \label{OEGF_Framp} \\
\underline{\wp}_{m} \leq \wp_{m}^{t} & \leq \overline{\wp}_{m} & m \in \mathcal{N} & \label{OEGF_PressureMinMax}\\
&\hspace{-3.45cm} \sum_{s \in \mathcal{S}_{m}} \phi_{sm}^{t} = \sum_{mn \in \mathcal{P}} \phi_{mn}^{{\rm in},t} - \sum_{nm \in \mathcal{P}} \phi_{nm}^{{\rm out},t} + \hspace{-1.5cm} & & \nonumber \\
&\hspace{-3.5cm} \sum_{mn \in \mathcal{E}} \phi_{mn}^{t} - \hspace{-0.1cm} \sum_{nm \in \mathcal{E}} \phi_{nm}^{t} + \hspace{-0.1cm} \sum_{g \in \mathcal{D}_{m}^{\rm g}} \phi_{gm}^{t} + \phi_{m}^{{\rm d},t} , \hspace{-1.5cm} & m \in \mathcal{N} & \label{OEGF_Nodal} \\
&\hspace{-3cm} \phi_{mn}^{t}\left|\phi_{mn}^{t} \right| = \Phi_{mn}\left((\wp_{m}^{t})^2 - (\wp_{n}^{t})^2 \right), \hspace{-1.5cm} & mn \in \mathcal{P} & \label{OEGF_gasflow_pipes}\\
\phi_{mn}^{t} = 0.5 & \left(\phi_{mn}^{{\rm in},t} + \phi_{mn}^{{\rm out},t} \right), \hspace{-1.5cm} & mn \in \mathcal{P} & \label{OEGF_averageflow}\\
&\hspace{-3cm} \tilde{\wp}_{mn}^{t} = \frac{2}{3} \left(\wp_{m}^{t} + \wp_{n}^{t} - \frac{\wp_{m}^{t}\wp_{n}^{t}}{\wp_{m}^{t} + \wp_{n}^{t}} \right), \hspace{-1.5cm} & mn \in \mathcal{P} & \label{OEGF_averagepressure} \\
\ell_{mn}^{t}=&\Psi_{mn}\tilde{\wp}_{mn}^{t}, & mn \in \mathcal{P} & \label{OEGF_linepack1} \\
\ell_{mn}^{t}=\ell_{mn}^{t-1} + \Delta \tau & \left( \phi_{mn}^{{\rm in},t} - \phi_{mn}^{{\rm out},t} \right), \hspace{-1.5cm} & mn \in \mathcal{P} & \label{OEGF_continuity} \\
-\overline{\phi}_{mn} \leq \phi_{mn}^{t},\phi_{mn}^{{\rm in},t}&,\phi_{mn}^{{\rm out},t} \leq \overline{\phi}_{mn}, \hspace{-1.5cm} & mn \in \mathcal{P} & \label{OEGF_FlowPipeMinMax}\\
&\hspace{-1.5cm} \phi_{mn}^{t} = \phi_{mn}^{+,t} + \phi_{mn}^{-,t}, \hspace{-1.5cm} & mn \in \mathcal{C} & \label{OEGF_compflowplusminus} \\
&\hspace{-1.5cm} \phi_{mn}^{{\rm c},t} = \nu_{mn} \phi_{mn}^{+,t} , \hspace{-1.5cm} & mn \in \mathcal{C} & \label{OEGF_compgastrap} \\
&\hspace{-2cm} 0 \le \phi_{mn}^{+,t} \le z_{mn}^{t}\overline{\phi}_{mn} , \hspace{-1.5cm} & mn \in \mathcal{C} & \label{OEGF_compplus} \\
&\hspace{-2.5cm} -\left(1-z_{mn}^{t}\right)\overline{\phi}_{mn} \le \phi_{mn}^{-,t} \le 0 , \hspace{-1.5cm} & mn \in \mathcal{C} & \label{OEGF_compminus} \\
&\hspace{-3.5cm} -\left(1-z_{mn}^{t}\right)\overline{\phi}_{mn} \leq \phi_{mn}^{t} \leq z_{mn}^{t}\overline{\phi}_{mn}, \hspace{-1.5cm} & mn \in \mathcal{E} & \label{OEGF_compflow} \\
&\hspace{-3.5cm} \left(\overline{\wp}_{n} - \underline{\wp}_{m} \overline{\gamma}_{mn}\right)\left(z_{mn}^{t} - 1\right) + \wp_{n}^{t} \le \wp_{m}^{t} \overline{\gamma}_{mn}, \hspace{-1cm} & mn \in \mathcal{E} & \label{OEGF_compboost1} \\
&\hspace{-3.5cm} \left(\overline{\wp}_{m}\underline{\gamma}_{mn} - \underline{\wp}_{n} \right)\left(z_{mn}^{t} - 1\right) + \wp_{m}^{t}\underline{\gamma}_{mn} \le \wp_{n}^{t}, \hspace{-1cm} & mn \in \mathcal{E} & \label{OEGF_compboost2} \\
&\hspace{-2.5cm} \wp_{n}^{t} - \wp_{m}^{t} \leq z_{mn}^{t}\left( \overline{\wp}_{n} - \underline{\wp}_{m} \right), \hspace{-1.5cm} & mn \in \mathcal{E} & \label{OEGF_compnoboost1} \\
&\hspace{-2.5cm} \wp_{m}^{t} - \wp_{n}^{t} \leq z_{mn}^{t}\left(\overline{\wp}_{m} - \underline{\wp}_{n} \right), \hspace{-1.5cm} & mn \in \mathcal{E} & \label{OEGF_compnoboost2} \\
&\hspace{-1.5cm} z_{mn}^{t} \in \left\lbrace 0 , 1 \right\rbrace, & \hspace{-3cm} mn \in \mathcal{E} & \label{OEGF_binary}
\end{align}

for $t \in \mathcal{T}=\{t_{0},t_{0}+\Delta t,\ldots,t_{0} + \Delta t \left( H - 1 \right) \}$, where $\Delta \tau = 3600 \Delta t$, and $\mathcal{G} = \mathcal{G}^{\rm g} \cup \mathcal{G}^{\rm ng}$. The cost functions of electrical generators in the objective function in (<ref>) are assumed to be quadratic, and of the form $f_{gi}\left(p_{gi}^{t}\right)=c_{2,gi}\left(p_{gi}^{t}\Delta t \right)^2 + c_{1,gi}\left(p_{gi}^{t}\Delta t \right) + c_{0,gi} \Delta t$. The cost functions of the gas supplies are assumed to be linear, and of the form $f_{sm}\left(\phi_{sm}^{t}\right) = c_{sm}\phi_{sm}^{t} \Delta \tau$. Electrical network constraints are delineated by <ref>, where <ref> and <ref> capture the limits on active power generation and the generator ramp rates, respectively, whereas (<ref>) and (<ref>) capture the limits on branch active power and angle difference, respectively. KCL is represented by <ref> and the strengthened DC power flow constraints are delineated by <ref> and <ref>, where $g_{ij}=\Re\left\{Y_{ij}^*/T_{ij}\right\}$, $b_{ij}=\Im\left\{Y_{ij}^*/T_{ij}\right\}$, $g_{ji}=\Re\left\{Y_{ji}^* / T_{ji}^*\right\}$, and $b_{ji}=\Im\left\{Y_{ji}^* / T_{ji}^*\right\}$. The strengthened DC OPF formulation is intended to approximate transmission line losses by replacing the $\cos\left( \theta_{ij} \right)$ terms in the original AC constraints by their second-order Maclaurin series; i.e., $\cos\left( \theta_{ij} \right) \approx 1 - 0.5\theta_{ij}^2 $, and is therefore a better approximation compared to its vanilla DC OPF counterpart. Gas system constraints are delineated by <ref>, where <ref> delineates capacity limits of gas supplies, (<ref>) captures the nodal pressure limits, (<ref>) is the gas flow nodal balance, <ref> describe the quasi-dynamic behavior of gas flow in a pipe, and <ref> describe the operation of non-pipe elements, i.e., compressors and pressure regulators. In more detail, (<ref>) is the discretized equation of motion along the full length of the pipe [44], where

\begin{align*}
\Phi_{mn} = \frac{\pi^2 D_{mn}^5}{16 \rho^2 Z_{mn} R T L_{mn} f_{mn}},
\end{align*}

and $f_{mn} = 4\left(20.621D_{mn}^{1/6} \right)^{-2}$ defines the Weymouth friction factor [45]. In this work, the compressibility factor is defined as

\begin{multline*}
\hspace{-0.3cm} Z_{mn} \hspace{-0.1cm} = \hspace{-0.1cm} \left( 1 + \frac{49.9511 \left( 10^{1.785 S}\right) \left(\tilde{\wp}_{m} + \tilde{\wp}_{n} - \frac{\tilde{\wp}_{m}\tilde{\wp}_{n}}{\tilde{\wp}_{m} + \tilde{\wp}_{n}} \right) }{\frac{3\left(1.8 T \right)^{3.825}}{2} }\right)^{-1} \hspace{-0.5cm},
\end{multline*}

where $\tilde{\wp}_{m} = 0.5 \left(\underline{\wp}_{m} + \overline{\wp}_{m} \right) $, as opposed to the gross simplification of assuming a constant compressibility factor across the whole network. The average flow and pressure across a pipe are captured by (<ref>) and (<ref>), respectively. Constraint (<ref>) is the linepack equation with

\begin{align*}
\Psi_{mn} = \frac{\pi D_{mn}^2 L_{mn}}{4\rho Z_{mn} R T},
\end{align*}

and (<ref>) is the discretized continuity equation over the full length of the pipe. The limits on the gas flowing in a pipeline are captured by (<ref>). The operation of bidirectional compressors and pressure regulators is captured by <ref>, where $z_{mn}^{t} = 1$ if the flow is in the direction of pressure boost (drop) for a compressor (pressure regulator), and $z_{mn}^{t} = 0$ otherwise. Note that $\underline{\gamma}_{mn} \geq 1$ for a compressor and $\overline{\gamma}_{mn} \leq 1$ for a pressure regulator. Similarly, $\overline{\gamma}_{mn} > 1$ for a compressor and $0 < \underline{\gamma}_{mn} < 1$ for a pressure regulator. The gas trapped in a compressor ($\phi_{mn}^{{\rm c},t}$) when it is boosting pressure is captured by (<ref>), where parameter $\nu_{mn}$ encapsulates the percentage of $\phi_{mn}^{+,t}$ that is trapped in the compressor. It is identified in [46] that the gas trapped by a compressor typically ranges between 3% to 5% of the gas flowing across it. The two systems are coupled by (<ref>), which assumes a linear relationship between the power output of a GPG and its input gas consumption.

The nonconvexity of the problem stems from constraints <ref>, <ref>, and the disjoint sets in <ref>.

Problem <ref> can be written in the general form

\begin{align}
\underset {\substack{x}}
{\mbox{ minimize}} \quad f_{0}\left(x\right) & & \label{OEGF_Generalproblem_objective}\\
\text{ subject to} \quad h_{i}\left(x\right) & = 0, & & i=1,\ldots,p_{\rm e} \label{OEGF_Generalproblem_eq_elec}\\
h_{i}\left(x\right) & = 0, & & i=p_{\rm e}+1,\ldots,p \label{OEGF_Generalproblem_eq_gas}\\
a_{i}^{T}x & \leq b_{i}, & & i=1,\ldots,o \label{OEGF_Generalproblem_linear} \\
x_{i} & \in \left\lbrace 0,1 \right\rbrace, & & i=1,\ldots,r \label{OEGF_Generalproblem_integrality}
\end{align}

where $x \in \reals^{n}$, $f_{0}$,$h_{1},\ldots,h_{p}: \reals^{n} \rightarrow \reals$, $a_{1},\ldots,a_{o} \in \reals^{n}$, $b_{1},\ldots,b_{o} \in \reals$, and $r < n$.[Recall that an affine equality constraint of the form $a_{i}^{T}x = b_{i}$ can be rewritten as $\left\lbrace a_{i}^{T}x \leq b_{i}\right\rbrace \cap \left\lbrace - a_{i}^{T}x \leq -b_{i} \right\rbrace $.] Function $f_{0}\left(x\right)$ represents the cost functions in (<ref>), whereas (<ref>) represent <ref> and <ref>, and (<ref>) represent <ref> and <ref>. Equality constraints <ref>, <ref>, <ref>, and <ref> are nonlinear and therefore nonconvex. The linear constraints in (<ref>) represent <ref>, <ref>, <ref>, and <ref>, whereas (<ref>) represents integrality constraints (<ref>). Problem <ref> is an MINLP problem with a nonconvex continuous relaxation,[A continuous relaxation is obtained by relaxing the integrality constraints $z \in \left\lbrace 0,1 \right\rbrace$ into box constraints $z \in \left[ 0,1 \right] $.] which, compared to convex MINLP problems, is more challenging to solve using current state-of-the-art MINLP technology.
This work therefore introduces (i) an iterative MILP-based method that leverages the computational efficiency of MILP solvers and the accuracy of IPM, and (ii) a fast iterative LP-based method that can be used in situations where only LP solvers are available. The optimality of the solution is measured by comparing it to the one obtained from the MICP relaxation described in the next section.

§ MIXED-INTEGER CONVEX RELAXATION

An MICP relaxation of Problem <ref> can be obtained by convexifying constraints <ref>, <ref>, and <ref>. Quadratic equality constraints <ref> and <ref> can be straightforwardly convexified by relaxing them into inequality constraints of the form

\begin{align}
p_{ij}^{t} \ge g_{ij} \hspace{0.5cm} & \hspace{-0.5cm} 0.5(\theta_{ij}^{t})^2 + b_{ij}\theta_{ij}^{t}, \hspace{-1.5cm} & ij \in \mathcal{L} & \label{OEGF_MISOCP_Pij} \\
p_{ji}^{t} \ge g_{ji} \hspace{0.5cm} & \hspace{-0.5cm} 0.5(\theta_{ij}^{t})^2 - b_{ji}\theta_{ij}^{t}, \hspace{-1.5cm} & ji \in \mathcal{L}&^{\rm t}. \hspace{-0.1cm} \label{OEGF_MISOCP_Pji}
\end{align}

Constraints <ref> can be equivalently rewritten as

\begin{align}
\left( \phi_{mn}^{t}\right)^2 = & \Phi_{mn} \left( 2z_{mn}^{t}-1 \right) \left((\wp_{m}^{t})^2 - (\wp_{n}^{t})^2 \right), \hspace{-0.2cm} & mn \in \mathcal{P} & \label{OEGF_gasflow_pipes_alt} \\
-\overline{\phi}_{mn} & \left( 1 - z_{mn}^{t}\right) \le \phi_{mn} \le z_{mn}^{t} \overline{\phi}_{mn}, & mn \in \mathcal{P} & \label{OEGF_MISOCP_FlowPipeMinMax} \\
\wp_{m}^{t} & - \wp_{n}^{t} \le z_{mn}^{t}\left( \overline{\wp}_{m} - \underline{\wp}_{n} \right) , & mn \in \mathcal{P} & \label{OEGF_MISOCP_PressureMax} \\
\wp_{m}^{t} & - \wp_{n}^{t} \ge \left( 1 - z_{mn}^{t}\right) \left( \underline{\wp}_{m} - \overline{\wp}_{n}\right) , & mn \in \mathcal{P} & \label{OEGF_MISOCP_PressureMin} \\
& z_{mn}^{t} \in \left\lbrace 0 , 1 \right\rbrace, & \hspace{-3cm} mn \in \mathcal{P} & \label{OEGF_MISOCP_binary}
\end{align}

where $z_{mn}^{t}$ is a binary variable that takes the value of 1 when the gas is flowing from node $m$ to node $n$ and the value of 0 when the gas is flowing from node $n$ to node $m$. Next, (<ref>) can be transformed into an SOC constraint by introducing two new variables, $\xi_{mn}^{t}$ and $\zeta_{mn}^{t}$, and constraints of the form

\begin{align}
\left( \phi_{mn}^{t}\right)^2 = & \Phi_{mn} \left((\xi_{mn}^{t})^2 - (\zeta_{mn}^{t})^2 \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_gasflow_pipes}\\
\xi_{mn}^{t} - \wp_{m}^{t} & \le \left( 1 - z_{mn}^{t}\right) \left( \overline{\wp}_{n} - \underline{\wp}_{m} \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_xim_max} \\
\xi_{mn}^{t} - \wp_{m}^{t} & \ge \left( 1 - z_{mn}^{t}\right) \left( \underline{\wp}_{n} - \overline{\wp}_{m} \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_xim_min} \\
\xi_{mn}^{t} - \wp_{n}^{t} & \le z_{mn}^{t} \left( \overline{\wp}_{m} - \underline{\wp}_{n} \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_xin_max} \\
\xi_{mn}^{t} - \wp_{n}^{t} & \ge z_{mn}^{t} \left( \underline{\wp}_{m} - \overline{\wp}_{n} \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_xin_min} \\
\zeta_{mn}^{t} - \wp_{n}^{t} & \le \left( 1 - z_{mn}^{t}\right) \left( \overline{\wp}_{m} - \underline{\wp}_{n} \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_zetan_max} \\
\zeta_{mn}^{t} - \wp_{n}^{t} & \ge \left( 1 - z_{mn}^{t}\right) \left( \underline{\wp}_{m} - \overline{\wp}_{n} \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_zetan_min} \\
\zeta_{mn}^{t} - \wp_{m}^{t} & \le z_{mn}^{t} \left( \overline{\wp}_{n} - \underline{\wp}_{m} \right), & mn \in \mathcal{P} & \label{OEGF_MISOCP_zetam_max} \\
\zeta_{mn}^{t} - \wp_{m}^{t} & \ge z_{mn}^{t} \left( \underline{\wp}_{n} - \overline{\wp}_{m} \right). & mn \in \mathcal{P} & \label{OEGF_MISOCP_zetam_min}
\end{align}

An MISOCP relaxation of (<ref>) can now be obtained by relaxing (<ref>) into the convex SOC constraint

\begin{align}
\left( \xi_{mn}^{t} \right)^2 \ge & \left( \zeta_{mn}^{t} \right)^2 + \left( \phi_{mn}^{t}\right)^2/\Phi_{mn} . & mn \in \mathcal{P} & \label{OEGF_MISOCP_gasflow_pipes_r}
\end{align}

More interestingly, and perhaps less obvious at first glance, because the nodal pressures are strictly positive, constraint <ref> becomes convex when it is relaxed into an inequality constraint of the form

\begin{align}
& \tilde{\wp}_{mn}^{t} \ge \frac{2}{3} \left(\wp_{m}^{t} + \wp_{n}^{t} - \frac{\wp_{m}^{t}\wp_{n}^{t}}{\wp_{m}^{t} + \wp_{n}^{t}} \right). & mn \in \mathcal{P} & \label{OEGF_averagepressure_r}
\end{align}

The function $f(x,y) = - xy/(x+y)$ with $\dom f = \reals_{++}^{2} = \left\lbrace (x,y) \in \reals^{2}| x > 0, y > 0 \right\rbrace $ is a convex function, which entails that its epigraph $\epi f = \left\{ (x,y,z) | (x,y) \in \dom f, f(x,y) \le z \right\}$ is a convex set.

The proof proceeds by showing that the Hessian of $f(x,y)$ is positive semidefinite, i.e.,

\begin{align}
\hspace{-0.3cm} \nabla^{2}f(x,y) =
\begin{bmatrix}
2y^{2}/(x+y)^{3} & -2xy/(x+y)^{3} \\
-2xy/(x+y)^{3} & 2x^{2}/(x+y)^{3}
\end{bmatrix} \succeq 0.
\end{align}

Since $x > 0$ and $y > 0$, one can show that $z^{T}\nabla^{2}f(x,y)z \ge 0$ for all $z \in \reals^{2}$. In more detail, the above reduces to proving that $z_{1}^2y^2 + z_{2}^2x^2 - 2z_{1}z_{2}xy \ge 0$, which can be straightforwardly verified since $x^2 + y^2 -2xy \ge 0$ and $z_{1}^2 + z_{2}^2 - 2z_{1}z_{2} \ge 0$. This completes the proof.

An MICP relaxation of Problem <ref> can now be formally written as

\begin{align}
\hspace{-1.25cm} \underset{\substack{p_{gi}^{t},p_{ij}^{t},\theta_{i}^{t},\wp_{m}^{t}\\\tilde{\wp}_{mn}^{t},\phi_{mn}^{t},\ell_{mn}^{t},\phi_{mn}^{{\rm in},t}\\\phi_{nm}^{{\rm out},t},z_{mn}^{t},\phi_{gm}^{t},\phi_{sm}^{t}\\\phi_{nm}^{{\rm c},t},\phi_{nm}^{+,t},\phi_{nm}^{-,t},\zeta_{mn}^{t},\xi_{mn}^{t}}}
{\mbox{ minimize}} & \sum_{t \in \mathcal{T}} \Bigg( \sum_{gi \in \mathcal{G}} \hspace{-0cm} f_{gi}\left(p_{gi}^{t}\right)\Bigg. + \hspace{-0.2cm} \hspace{-0cm} & \hspace{-0cm} & \nonumber
\end{align}

\begin{align}
& \hspace{3cm} \Bigg. \sum_{sm\in \mathcal{S}} f_{sm}\left(\phi_{sm}^{t}\right)\Bigg) \hspace{-3cm} & & \label{OEGF_MISOCP_objective} \\
& \hspace{-0.75cm} \text{subject to \cref{OEGF_Pminmax,OEGF_Pramp,OEGF_MWlimits,OEGF_anglediff,OEGF_KCLp}, \cref{OEGF_gasflow_GFG,OEGF_Fminmax,OEGF_PressureMinMax,OEGF_Nodal},  \cref{OEGF_linepack1,OEGF_continuity,OEGF_FlowPipeMinMax,OEGF_compflowplusminus,OEGF_compgastrap,OEGF_compplus,OEGF_compminus,OEGF_compflow,OEGF_compboost1,OEGF_compboost2,OEGF_compnoboost1,OEGF_compnoboost2}}, \hspace{-2cm} & & \nonumber \\
& \hspace{0.75cm} \text{\cref{OEGF_MISOCP_Pij,OEGF_MISOCP_Pji}, \cref{OEGF_MISOCP_FlowPipeMinMax,OEGF_MISOCP_PressureMax,OEGF_MISOCP_PressureMin,OEGF_MISOCP_binary}, \cref{OEGF_MISOCP_xim_max,OEGF_MISOCP_xim_min,OEGF_MISOCP_xin_max,OEGF_MISOCP_xin_min,OEGF_MISOCP_zetan_max,OEGF_MISOCP_zetan_min,OEGF_MISOCP_zetam_max,OEGF_MISOCP_zetam_min,OEGF_MISOCP_gasflow_pipes_r,OEGF_averagepressure_r}}. \hspace{-2cm} & & \label{OEGF_MISOCP_common}
\end{align}

Problem <ref> is a tractable convex MINLP problem[An MINLP problem is called “convex MINLP” if it becomes a convex NLP when the integrality constraints are relaxed.] that can be solved using the (iterative) outer approximation (OA) [47]. Alternatively, if constraint (<ref>) is replaced by a tight polyhedral envelope using a sufficiently large number of halfspaces, Problem <ref> can be transformed into an MISOCP that can be directly handled by powerful MISOCP solvers such as Gurobi [48]. The construction of such a tight polyhedral envelope is discussed in the next section. Numerical evaluation in Section <ref> will show that, although tractable compared to solving Problem <ref> directly using NLBB solvers, the MICP problem in (<ref>) converges to infeasible solutions to the original problem in (<ref>). Nonetheless, this MICP relaxation provides a good-quality lower bound for assessing the optimality of solutions to Problem <ref> from methods such as NLBB solvers and the ones proposed in this work.

§ POLYHEDRAL ENVELOPES

As a first step towards formulating a computationally efficient LP approximation of Problem <ref>, the below discourse introduces polyhedral envelopes for the nonlinear terms in <ref>, <ref>, and a PWL approximation of the quadratic cost functions in (<ref>).

§.§ Polyhedral envelope of $x^{2}$

A polyhedral envelope of a nonconvex set of the form $\mathcal{V}= \left\{x^{2} | x \in \left[\underline{x},\overline{x}\right]\right\}$, can be obtained by, (i), finding the supporting hyperplanes passing through points $\left( \underline{x},\underline{x}^2 \right) $ and $\left( \overline{x},\overline{x}^2 \right) $, and (ii), finding the supporting hyperplanes obtained by outer PWL approximations of the set. In more detail, $l$ points $x_1,\ldots,x_l$ are selected in the interval $\left[\underline{x},\overline{x}\right]$, which allows adding $(l-1) + 1$ halfspaces of the form

\[$$
\begin{subnumcases}{\label{convV} \hspace{-0mm} {\rm conv}\left(\mathcal{V}\right)= \!}
\! v \geq \left(2 x_{h}\right)x-x_{h}^{2}, \quad h=\left\{1,\ldots,l\right\}, \label{convV_v1}\\
\! v \leq \left(\overline{x}+\underline{x}\right)x-\overline{x} \underline{x} . \label{convV_v2}
\end{subnumcases}
$$ \nonumber\]

§.§ Polyhedral envelope of $x \left| x \right|$

Polyhedral envelope of $\mathcal{\mathcal{U}}= \left\{x \left| x \right| | x \in \left[\underline{x},\overline{x}\right]\right\}$.

A polyhedral envelope of a nonconvex set of the form $\mathcal{\mathcal{U}}= \left\{x \left| x \right| | x \in \left[\underline{x},\overline{x}\right]\right\}$ is obtained by determining the parameters of the four supporting hyperplanes passing through points $\left( \underline{x},\underline{x} \left|\underline{x} \right| \right) $ and $\left( \overline{x},\overline{x}\left|\overline{x} \right| \right) $. This can be done by separately solving two nonlinear equations in one dimension for each gas pipe. The first supporting hyperplane passing through $\underline{x}$ has the form $\left.g(x)\right|_{\underline{x}} = \underline{x} \left|\underline{x} \right| + 2\left| \underline{x}\right|\left(x - \underline{x} \right)$. The second supporting hyperplane passing through $\underline{x}$ can be obtained by solving the one-dimensional equation $\left.g(z)\right|_{\underline{x}} = z\left| z\right| + 2\left| z\right| \left(\underline{x}- z \right) = \underline{x} \left|\underline{x} \right|$, and	whose solution is denoted by $z^{\star} = \alpha = \underline{x}(1-\sqrt{2})$. Analogously, the first supporting hyperplane passing through $\overline{x}$ has the form $\left.g(x)\right|_{\overline{x}} = \overline{x} \left|\overline{x} \right| + 2\left| \overline{x}\right| \left(x - \overline{x} \right)$. The second supporting hyperplane passing through $\overline{x}$ can be obtained by solving $\left.g(z)\right|_{\overline{x}} = z\left| z\right| + 2\left| z\right| \left(\overline{x}- z \right) = \overline{x} \left|\overline{x} \right|$, and whose solution is denoted by $z^{\star} = \beta = \overline{x}(1-\sqrt{2})$. More supporting hyperplanes can be obtained in intervals $[\underline{x}, \beta]$ and $[\alpha, \overline{x}]$ to obtained a tighter polyhedral envelope but in this work only $\left.g(x)\right|_{\gamma}$ and $\left.g(x)\right|_{\delta}$ are added, where $\gamma$ and $\delta$ are obtained from the intersections of $\left.g(x)\right|_{\alpha}$ and $\left.g(x)\right|_{\overline{x}}$, and $\left.g(x)\right|_{\beta}$ and $\left.g(x)\right|_{\underline{x}}$, respectively. As a result, the polyhedral envelope of the nonconvex set $\mathcal{U}$ can be written as
conv(𝒰) =
u ≥x( 2 - √(8) )x - x^2( 3 - √(8) ),

u ≥2xx - x^2,

u ≥γ| γ| + 2| γ| (x- γ),

u ≤x(√(8) - 2)x + x^2(3 - √(8)),

u ≤-2xx + x^2,

u ≤δ| δ| + 2| δ| (x - δ).

An illustration of (<ref>) is shown in Fig. <ref>.[This formulation is valid for $\underline{x} < 0$ and $\overline{x} > 0$.]

§.§ Polyhedral envelope of the average pressure

Polyhedral envelope of $\mathcal{W}= \left\{ f(x,y) = -xy/(x+y) \right.$ $ \left. | \left(x,y\right) \in \left[\underline{x},\overline{x}\right] \times\left[\underline{y},\overline{y}\right]\right\}$.

Following Theorem <ref>, a tight polyhedral envelope of a nonconvex set of the form $\mathcal{W}= \left\{ f(x,y) \right. $ $\left. = -xy/(x+y) | \left(x,y\right) \in \left[\underline{x},\overline{x}\right] \times\left[\underline{y},\overline{y}\right]\right\}$ can be obtained from the halfspaces associated with a sufficiently large number of uniformly spaced points in the 2-D grid defined by $\dom f$. However, since the aim of this construction is computational efficiency and not tightness, only the four extreme points in the 2-D grid defined by $\dom f$ are chosen for the construction of the polyhedral envelope. These four extreme points are $p1 = \left(\underline{x},\underline{y}, f\left(\underline{x},\underline{y}\right)\right) $, $p2 = \left(\overline{x},\underline{y}, f\left(\overline{x},\underline{y}\right)\right) $, $p3 = \left(\underline{x},\overline{y}, f\left(\underline{x},\overline{y}\right)\right) $, and $p4 = \left(\overline{x},\overline{y}, f\left(\overline{x},\overline{y}\right)\right) $. In more detail, the first supporting hyperplane in the 3-D Euclidean space, which can be obtained from points $p1$, $p2$, and $p3$, is written as $g\left( x,y\right) = f\left(\underline{x},\underline{y}\right) - \frac{a_{x}}{a_{z}}\left( x - \underline{x} \right) - \frac{a_{y}}{a_{z}}\left( y - \underline{y} \right)$, where $a = \left(p2-p1 \right) \times \left(p3-p1 \right) $ is the normal vector. Similarly, the second supporting hyperplane can be obtained from points $p2$, $p3$, and $p4$, and can be written as $g\left( x,y\right) = f\left(\overline{x},\overline{y}\right) - \frac{b_{x}}{b_{z}}\left( x - \overline{x} \right) - \frac{b_{y}}{b_{z}}\left( y - \overline{y} \right)$, where $b = \left(p2-p4 \right) \times \left(p3-p4 \right) $ is the normal vector. The remaining four supporting hyperplanes can be obtained from the first-order Taylor series approximation of $f(x,y)$ at the four extreme points $p1$, $p2$, $p3$, and $p4$. The polyhedral envelope of $\mathcal{W}$ can now be written as

conv(𝒲) =

w ≥f(x,y) + ∂f(x,y)/∂x ( x - x ) + ∂f(x,y)/∂y ( y - y ),

w ≥f(x,y) + ∂f(x,y)/∂x ( x - x ) + ∂f(x,y)/∂y ( y - y ),

w ≥f(x,y) + ∂f(x,y)/∂x ( x - x ) + ∂f(x,y)/∂y ( y - y ),

w ≥f(x,y) + ∂f(x,y)/∂x ( x - x ) + ∂f(x,y)/∂y ( y - y ),

w ≤f(x,y) - a_x/a_z( x - x ) - a_y/a_z( y - y ),

w ≤f(x,y) - b_x/b_z( x - x ) - b_y/b_z( y - y ).

An illustration of (<ref>) is shown in Fig. <ref>.

§.§ PWL approximation of a quadratic cost function

The last step is to substitute the quadratic terms in $f_{gi}\left( p_{gi}^{t} \right)$ by corresponding variables and rotated SOC constraints which can then be tightly approximated by a lifted polyhedron as in [49]. This construction requires far less inequality constraints than standard piecewise linear (PWL) approximations for the same accuracy. Specifically, $\left|\mathcal{G}\right| \times \left|\mathcal{T}\right|$ variables $p_{g}^{t}$ and associated linear constraints of the form $p_{g}^{t}=\Delta t \sqrt{c_{2,gi}}p_{gi}^{t}$, for all $ gi \in \mathcal{G}$, $t \in \mathcal{T}$, are introduced along with $N \times \left|\mathcal{T}\right|$ variables $\delta_{n}^{t}$ and constraints of the form $\delta_{n}^{t} \geq ( p_{2n-1}^{t} )^2 + ( p_{2n}^{t} )^2 $, for which the lifted polyhedral construction is denoted as

\begin{align}\label{RSOC_Obj}
\hspace{-0.2cm} \mathcal{P}_{k}^{\rm R} \left( \left( p_{2n-1}^{t}\right)^2 + \left( p_{2n}^{t} \right)^2 \leq \delta_{n}^{t} \right), n \in \left\{1,\ldots,N\right\}, t \in \mathcal{T}
\end{align}

where $N = \left\lfloor \left|\mathcal{G}\right| /2\right\rfloor+\left\lceil \left| \mathcal{G} \right| /2 - \left\lfloor \left| \mathcal{G} \right| /2\right\rfloor \right\rceil$.

§.§ Polyhedral relaxation of the MOEGF

In light of the above, the LP-relaxed MOEGF with the polyhedral envelopes and PWL objective function can be written as

\begin{align}
\hspace{-0.5cm} \underset{\substack{p_{gi}^{t},\delta_{n}^{t},p_{ij}^{t},\theta_{i}^{t},v_{ij}^{t},\wp_{m}^{t}\\v_{m}^{t},\tilde{\wp}_{mn}^{t},w_{mn}^{t},\phi_{mn}^{t},u_{mn}^{t}\\\ell_{mn}^{t},\phi_{mn}^{{\rm in},t},\phi_{nm}^{{\rm out},t},z_{mn}^{t}\\\phi_{gm}^{t},\phi_{sm}^{t},\phi_{nm}^{{\rm c},t},\phi_{nm}^{+,t},\phi_{nm}^{-,t} }}
{\mbox{ minimize}} & \sum_{t \in \mathcal{T}} \Bigg( \sum_{n=1}^{N} \delta_{n}^{t} + \hspace{-0.0cm} \sum_{gi \in \mathcal{G}} \hspace{-0cm} c_{1,gi}\left(p_{gi}^{t}\Delta t \right) \Bigg. \hspace{-0.2cm} & \hspace{-0cm} & \nonumber
\end{align}

\begin{align}
& \hspace{2cm} + c_{0,gi} \Delta t + \Bigg. \sum_{sm\in \mathcal{S}} f_{sm}\left(\phi_{sm}^{t}\right)\Bigg) \hspace{-3cm} & & \label{OEGF_Poly_objective} \\
& \hspace{-0.75cm} \text{subject to \cref{OEGF_Pminmax,OEGF_Pramp,OEGF_MWlimits,OEGF_anglediff,OEGF_KCLp}, \cref{OEGF_gasflow_GFG,OEGF_Fminmax,OEGF_PressureMinMax,OEGF_Nodal}, \cref{OEGF_averageflow}, \cref{OEGF_linepack1,OEGF_continuity,OEGF_FlowPipeMinMax,OEGF_compflowplusminus,OEGF_compgastrap,OEGF_compplus,OEGF_compminus,OEGF_compflow,OEGF_compboost1,OEGF_compboost2,OEGF_compnoboost1,OEGF_compnoboost2}, \cref{RSOC_Obj}} \hspace{-2cm} & & \label{OEGF_Poly_common} \\
p_{ij}^{t}&=g_{ij}0.5{\rm conv}\left( \left(\theta_{ij}^{t}\right)^2 \right) + b_{ij}\theta_{ij}^{t}, & \hspace{-0.0cm} ij \in \mathcal{L} & \label{OEGF_Poly_Pij} \\
p_{ji}^{t}&= g_{ji}0.5{\rm conv}\left( \left(\theta_{ij}^{t}\right)^2 \right) - b_{ji}\theta_{ij}^{t}, & \hspace{-0.0cm} ji \in \mathcal{L}&^{\rm t} \label{OEGF_Poly_Pji} \\
{\rm conv}&\left( \phi_{mn}^{t}\left|\phi_{mn}^{t} \right|\right) = \Phi_{mn} {\rm conv}\left( \left(\wp_{m}^{t}\right)^2 \right) \hspace{-0.5cm} & & \nonumber \\
& \hspace{2.5cm} - {\rm conv}\left( \left(\wp_{n}^{t}\right)^2 \right), & mn \in \mathcal{P} & \label{OEGF_Poly_gasflow_pipes}\\
& \hspace{-1cm} \tilde{\wp}_{mn}^{t} = \frac{2}{3} \left(\wp_{m}^{t} + \wp_{n}^{t} + {\rm conv}\left( \frac{-\wp_{m}^{t}\wp_{n}^{t}}{\wp_{m}^{t} + \wp_{n}^{t}} \right) \right), \hspace{-0.3cm} & mn \in \mathcal{P} & \label{OEGF_Poly_averagepressure} \\
&\hspace{2cm} z_{mn}^{t} \in \left[ 0 , 1 \right] , & \hspace{-0cm} mn \in \mathcal{E} & \label{OEGF_Poly_nobinary}
\end{align}

for $t \in \mathcal{T}$. These polyhedral envelopes depend only on the properties and topology of the network and can thus be formed once, in a single computationally-cheap preprocessing stage. The solution of Problem <ref> can be used as a warm start for the two iterative algorithms discussed in the next two sections.

§ ITERATIVE MILP-BASED METHOD

Iterative MILP-based method
[1]
Set $k=1$, $\overline{K} = 40$, $\epsilon \in \left[ 10^{-10}, 10^{-4} \right] $, $\sigma >> 0$ , $\overline{\alpha} > \alpha_{mn}^{t,(1)} > 0$ for all $ mn \in \mathcal{P}$ and $t \in \mathcal{T}$, $\gamma > 0$, and $\mathcal{I}_{i} = \left\lbrace \emptyset \right\rbrace$ for all $ i \in \left\lbrace 1,\ldots,p_{\rm e} \right\rbrace $. Obtain the initial point $x^{(1)}$ by either solving Problem <ref> (warm start) or by choosing $x^{(1)}$ randomly such that $x^{(1)}$ is in the domain of the all constraint functions of Problem <ref>. The point $x^{(1)}$ need not be feasible.
Phase I (LP)
${\rm max}\left( \left( \left| h_{i}\left( x^{(k)}\right) \right| \right)_{i=1,\ldots,p}\right) > \epsilon$

\begin{align}
& \hspace{-0.6cm} x^{(k+1)} = & & \nonumber \\
& \hspace{-0.6cm} \underset {\substack{x,r_{ij}^{t}}}
{\argmin} \ \eqref{OEGF_Poly_objective} + \sum_{t \in \mathcal{T}} \left( \sum_{ij \in \mathcal{L}} \sigma r_{ij}^{t} + \hspace{-0.1cm} \sum_{mn \in \mathcal{P}} \hspace{-0.0cm} \alpha_{mn}^{t,(k)} \left|\phi_{mn}^{t} - \phi_{mn}^{t,(k)} \right| \right) \hspace{-5cm} & & \label{OEGF_Taylor_objectiveI}\\
& \text{subject to } a_{i}^{T}x \leq b_{i}, \hspace{-0.0cm} & & i=1,\ldots,o \label{OEGF_Taylor_linear} \\
& \hspace{-0.5cm} h_{i}(x^{(k)}) + \nabla h_{i}(x^{(k)})^{T}(x - x^{(k)}) = r_{i}, \hspace{-0.0cm} & & i = 1,\ldots,p_{\rm e} \label{OEGF_Taylor_eq_elec_eq}\\
& \hspace{-0.5cm} h_{i}(x^{(\kappa)}) + \nabla h_{i}(x^{(\kappa)})^{T}(x - x^{(\kappa)}) \leq 0, \hspace{-0.0cm} & & i = 1,\ldots,p_{\rm e}, \kappa \in \mathcal{I}_{i} \label{OEGF_Taylor_eq_elec_ineq}\\
& \hspace{-0.5cm} r_{i} \ge 0, \hspace{-0.0cm} & & i = 1,\ldots,p_{\rm e} \label{OEGF_Taylor_eq_elec_ri}\\
& \hspace{-0.5cm} h_{i}(x^{(k)}) + \nabla h_{i}(x^{(k)})^{T}(x - x^{(k)}) = 0, \hspace{-0.31cm} & & i=p_{\rm e}+1,\ldots,p \label{OEGF_Taylor_eq_gas}\\
& \hspace{1cm} x_{i} \in \left[ 0,1 \right] , \hspace{-0.31cm} & & i=1,\ldots,r. \label{OEGF_Taylor_01_box}
\end{align}

$i = 1,\ldots,p_{\rm e}$
$h_{i}\left( x^{(k+1)}\right) < 0$
$\mathcal{I}_{i} := \mathcal{I}_{i} \cup \{k\} $.
$i = p_{\rm e}+1,\ldots,p$
$ \left| h_{i}\left( x^{(k+1)}\right) \right| \geq \left| h_{i}\left( x^{(k)}\right) \right| \textbf{and} \left| h_{i}\left( x^{(k+1)}\right) \right| > \epsilon $
$\alpha_{mn}^{t,(k+1)} := \min\left( \overline{\alpha}, \gamma \alpha_{mn}^{t,(k)} \right) $.
$k \leftarrow k + 1 $.
Phase II (MILP)
${\rm max}\left( \left( \left| h_{i}\left( x^{(k)}\right) \right| \right)_{i=1,\ldots,p}\right) > \epsilon$

\begin{align}
& \hspace{-1cm} x^{(k+1)} = \underset {\substack{x,r_{ij}^{t}}}
{\argmin} \ \eqref{OEGF_Taylor_objectiveI} \hspace{-5cm} & & \label{OEGF_Taylor_objectiveI_MILP}\\
& \hspace{-1cm} \text{subject to \cref{OEGF_Taylor_linear,OEGF_Taylor_eq_elec_eq,OEGF_Taylor_eq_elec_ineq,OEGF_Taylor_eq_elec_ri,OEGF_Taylor_eq_gas}} & & \label{OEGF_Taylor_Common} \\
& \hspace{0cm} x_{i} \in \left\lbrace 0,1 \right\rbrace, \hspace{-0.31cm} & & i=1,\ldots,r. \label{OEGF_Taylor_01}
\end{align}

Same as lines 4 to 13.
$k \leftarrow k + 1 $.

The proposed MILP-based method for directly solving Problem <ref> is described in Algorithm <ref>, which is divided into two phases. Phase I consists of a sequence of carefully  coordinated LP problems for solving the continuous relaxation of Problem <ref>, i.e., Problem <ref> with (<ref>) instead of (<ref>). The solution from Phase I is then used as a warm start for Phase II which consists of a sequence of carefully coordinated MILP problems that converges to a feasible solution to Problem <ref>. Phase I is akin to SLP which was originally introduced in [50, 51]. SLP, which is also known at the cutting plane method, consists of solving the original NLP problem by solving a series of LP problems generated by approximating all the nonlinear constraints by their first-order Taylor series expansion around the current point $x^{(k)}$, where $k$ is the iteration number. The iterative procedure in Phase I can either be initialized from the solution of Problem <ref> or from a random point that lies in the domain of the all constraint functions of Problem <ref>. The purpose of Problem <ref> is two fold. First, it provides a valid starting point for Algorithm <ref> when no prior information is available. Second, the superior quality of this starting point, compared to a random point, manifests in an improvement in overall convergence. The small computational overhead of Problem <ref> is far outweighed by the improvement in overall convergence, as will be demonstrated in Section <ref>.

The convergence of Phase I is enabled by two key features. The first consists of a finite set of closed halfspaces, delineated by (<ref>), implemented in conjunction with a finite set of hyperplanes of the form (<ref>). However, instead of directly using supporting hyperplanes of the form $\left\lbrace x | h_{i}(x^{(k)}) + \nabla h_{i}(x^{(k)})^{T}(x - x^{(k)}) = 0, i=1,\ldots,p_{\rm e} \right\rbrace $ (to the nonconvex set of the form $\mathcal{V}$), which might lead to infeasible LP problems, Algorithm <ref> introduces non-negative slack variables $r_{ij}^{t}$ whose purpose is to prevent infeasible LP problems, especially during the first few iterations of the algorithm if it is initialized from a poor-quality starting point (random $x^{(1)}$). These variables are then minimized in the objective by assigning a relatively large value for parameter $\sigma$. This construction can be viewed as an iteratively refined polyhedral outer approximation of the nonconvex sets of the form $\mathcal{V}= \left\{x^{2} | x \in \left[\underline{x},\overline{x}\right]\right\}$, but with the use of both supporting hyperplanes and halfspaces. Set $\mathcal{I}_{i}$ registers all previous iteration numbers that qualify as supporting halfspaces for $\mathcal{V}= \left\{x^{2} | x \in \left[\underline{x},\overline{x}\right]\right\}$. The condition on line 5 is satisfied when constraints (<ref>)-(<ref>) are violated at iteration $k$, in which case the supporting halfspace of the form (<ref>) is added to Problem <ref> at $k+1$. The second key feature is the term $\alpha_{mn}^{t,(k)} \left|\phi_{mn}^{t} - \phi_{mn}^{t,(k)} \right|$ in (<ref>), which controls the step size of the LP approximation of (<ref>). The step size parameter $\alpha_{mn}^{t,(k)}$ is automatically tuned on lines 9 to 13, predicated on the condition on line 10 which is true when the violations of constraints (<ref>) and (<ref>) do not decrease as the algorithm iterates.

Phase I terminates when the constraint with the largest violation, $C_{\rm max}^{k} = {\rm max}\left( \left( \left| h_{i}\left( x^{(k)}\right) \right| \right)_{i=1,\ldots,p}\right) $, does not exceed a certain tolerance $\epsilon$.[Feasibility here is measured with respect to Problem <ref> and not to the original problem with AC power flow constraints (which is not shown here).] The solution obtained at the termination of Phase I is likely to be non-integral, i.e., some of the $z_{mn}^{t}$ variables may not be strictly 0 or 1. This solution is therefore used as a warm start for Phase II which is iterative MILP-based algorithm that closely resembles Phase I but with the integrality constraints (<ref>) instead of the relaxed ones in (<ref>). Phase II, and therefore Algorithm <ref>, terminates under the same conditions in Phase I but now with a guarantee that the integrality constraints are satisfied. The next section introduces a fast iterative LP-based method that can be used as an alternative to Algorithm <ref> in situations where only LP solvers are available.

§ ITERATIVE LP-BASED METHOD

Iterative LP-based mehod
[1]
Set $k=1$, $\overline{K} = 40$, $k_{2} = 0$, $K_{f} = 2$, $\epsilon \in \left[ 10^{-10}, 10^{-4} \right] $, $\sigma >> 0$, $\beta > 0$, $\overline{\alpha} > \alpha_{mn}^{t,(1)} > 0$ for all $ mn \in \mathcal{P}$ and $t \in \mathcal{T}$, $\gamma > 0$, and $\mathcal{I}_{i} = \left\lbrace \emptyset \right\rbrace$ for all $ i \in \left\lbrace 1,\ldots,p_{\rm e} \right\rbrace $. Obtain the initial point $x^{(1)}$ by either solving Problem <ref> (warm start) or by choosing $x^{(1)}$ randomly such that $x^{(1)}$ is in the domain of the all constraint functions of Problem <ref>. The point $x^{(1)}$ need not be feasible.
Phase I (LP)
Same as lines 2 to 15 of Algorithm <ref>.
Phase II (LP)
$mn \in \mathcal{E}, t \in \mathcal{T}$
$\phi_{mn}^{t,(k)} \leq 0$
$I_{mn}^{t} \leftarrow 0 $.
$\phi_{mn}^{t,(k)} > 0$
$I_{mn}^{t} \leftarrow 1 $.
$ k \le \overline{K}$

\begin{align}
& \hspace{-1cm} x^{k+1} = \underset {\substack{x},r_{ij}^{t}}
{\argmin} \ \eqref{OEGF_Taylor_objectiveI} + \sum_{mn \in \mathcal{E}, t \in \mathcal{T}} \beta \left|z_{mn}^{t} - I_{mn}^{t} \right| \hspace{0cm} & & \label{OEGF_Taylor_objectiveII}\\
& \hspace{-0cm} \text{subject to \cref{OEGF_Taylor_linear,OEGF_Taylor_eq_elec_eq,OEGF_Taylor_eq_elec_ineq,OEGF_Taylor_eq_elec_ri,OEGF_Taylor_eq_gas,OEGF_Taylor_01_box}}. & &
\end{align}

${\rm max}\left( \left( \left| h_{i}\left( x^{(k)}\right) \right| \right)_{i=1,\ldots,p}\right) \le \epsilon$ and all $z_{mn}^{t,(k+1)} \in \left\lbrace 0,1 \right\rbrace $
$ \text{mod}\left(k_{2},K_{f} \right) = 0 $
$mn \in \mathcal{E}, t \in \mathcal{T}$
$I_{mn}^{t} == 0$ and $z_{mn}^{t,(k+1)} > \epsilon$
$I_{mn}^{t} \leftarrow 1 $.
$I_{mn}^{t} == 1$ and $z_{mn}^{t,(k+1)} < 1 - \epsilon$
$I_{mn}^{t} \leftarrow 0 $.
$k \leftarrow k + 1 $.
$k_{2} \leftarrow k_{2} + 1 $.

The LP-based heuristic, described in Algorithm <ref>, starts exactly like Algorithm <ref> by solving the continuous relaxation of Problem <ref> in Phase I, whose solution is then passed on to the LP-based heuristic in Phase II. The solution obtained at the termination of Phase I is likely to be non-integral, i.e., some of the $z_{mn}^{t}$ variables may not be strictly 0 or 1. Phase II therefore consists of a heuristic that ensures that the solution is integral and feasible within a tolerance $\epsilon$. In a nutshell, instead of a branch-and-bound algorithm, the heuristic in Phase II implements an elaborate steering procedure controlled by the term $\beta \left|z_{mn}^{t} - I_{mn}^{t} \right|$ in (<ref>). More specifically, lines 3 to 8 set parameter $I_{mn}^{t}$ for both compressors and pressure regulators to either 0 or 1 based on the direction of flow $\phi_{mn}^{t,(k)}$ from the solution of Phase I. These values of $I_{mn}^{t}$ are then passed on to Problem <ref>. The purpose of lines 15 to 23 is to handle any remaining non-integral values of $z_{mn}^{t}$ after solving (<ref>). There are two reasons a variable $z_{mn}^{t}$ can remain non-integral after solving (<ref>) at iteration $k$. The first pertains to the potential for decreasing the objective function value, and the second is related to feasibility. For instance, if $I_{mn}^{t}$ was set to 1 but (<ref>) resulted in $z_{mn}^{t,(k+1)} = 0.85$, the heuristic uses this as a clue to steer $z_{mn}^{t,(k+2)}$ towards 0. In other words, what this means is that had $z_{mn}^{t,(k+1)}$ been forced to 1, the problem would have either converged to a solution with a larger objective function value, or to an infeasible point. Analogously, if $I_{mn}^{t}$ was set to 0 but (<ref>) resulted in $z_{mn}^{t,(k+1)} = 0.25$ for instance, the heuristic uses this as a clue to steer $z_{mn}^{t,(k+2)}$ towards 1. In other words, what this means is that had $z_{mn}^{t,(k+1)}$ been forced to 0, the problem would have either converged to a solution with a larger objective function value, or to an infeasible point. This procedure is only invoked every $K_{f}$ iterations, as opposed to at every iteration, in an effort to prevent potential oscillatory behavior. In summary, the procedure in Phase II does not force the binary variables to either 0 or 1, but instead uses the term $\beta \left|z_{mn}^{t} -I_{mn}^{t} \right|$ to steer these variables to take binary values predicated on a careful tuning of parameter $\beta$. Therefore, the value of $\beta$ is chosen to strike a good tradeoff between objective value minimization and infeasibility handling. Phase II is terminated when $C_{\rm max}^{k} = {\rm max}\left( \left( \left| h_{i}\left( x^{(k)}\right) \right| \right)_{i=1,\ldots,p}\right) \leq \epsilon$ and when all the $z_{mn}^{t}$ variables are integral.

The procedure in Phase II, although technically a heuristic, is successful for this type of problem mainly because the direction of gas flow in compressors and regulators is predominantly determined by pipeline constraints (<ref>) and the location of the gas demands. The heuristic in Algorithm 1 therefore exploits this property of the problem to recover an integral solution. Just like other local MINLP solvers such as Juniper and KNITRO [52], Algorithm <ref> cannot guarantee global optimality.

§ NUMERICAL EVALUATION

Elements of the test cases.

Case	   	$\left| \mathcal{B} \right| $	   	$\left| \mathcal{L} \right| $	   	$\left| \mathcal{G}^{\rm ng} \right| $	   	$\left| \mathcal{G}^{\rm g} \right| $	   	$\left| \mathcal{N} \right| $	   	$\left| \mathcal{P} \right| $	   	$\left| \mathcal{C} \right| $	   	$\left| \mathcal{R} \right| $	   	$\left| \mathcal{S} \right| $
A	   	5	   	6	   	3	   	2	   	7	   	4	   	2	   	0	   	2
B	   	14	   	20	   	3	   	2	   	25	   	24	   	6	   	0	   	6
C	   	72	   	122	   	14	   	17	   	47	   	39	   	10	   	7	   	5

In this experimental setup, Julia v1.4.0 [53] is used as a programming language along with JuMP v0.21.1 [54] as a frontend modeling language for the optimization problems. Gurobi v9.0.2 [48] with the Barrier method (without a crossover strategy) is used to solve all the LP problems (including Problem <ref>). The MICP optimization problems in (<ref>) are also solved using Gurobi v9.0.2 but now with the linearized outer approximation algorithm. The continuous relaxation of Problem <ref>, which is a nonconvex NLP problem, is solved using IPOPT v3.12.11 [55] and linear solver MA57 [56]. The original MINLP problem in (<ref>) is solved using the NLBB solver Juniper [57] with IPOPT v3.12.11 [55] and linear solver MA57 [56] for the NLP subproblems, and Gurobi v9.0.2 in the feasibility pump heuristic at the root node. All simulations are conducted on a computing platform with an Intel Core i7-6820HK CPU at 2.7GHz, 64-bit operating system, and 32GB RAM. Five test systems are considered for the numerical evaluation of Algorithm <ref> with hourly ($\Delta t = 1$) demand and linepack data over $H = 24$ hours. Test case A consists of a 7-node gas system connected to a 5-bus electrical system, test case B consists of the Belgian gas system [15] connected to the IEEE 14-bus electrical system [40], and test cases C1, C2, and C3 consist of the real electricity and gas system of the state of Victoria, Australia for low (23/11/2019), medium (21/08/2019), and high (09/08/2019) demand days, respectively. Data for the actual gas network of the state of Victoria, the Victorian Declared Transmission System (DTS), was developed from scratch with the help of industry support within the Future Fuels CRC project [58]. The electrical network data for the state of Victoria is obtained from [59] and updated to reflect the generation mix and network augmentations of 2019. A summary description of the three IEGS test cases is shown in Table <ref>, and detailed data can be found in [31].

In Problem <ref>, functions $f_{gi}\left( p_{gi}^{t}\right) $ in (<ref>) are approximated by a lifted polyhedron with accuracy of $1.15 \times 10^{-9}$ (see [49]), whereas $l=10$ is chosen for the outer PWL approximation in the polyhedral envelopes of the quadratic terms in <ref>, and <ref>. The angle limits are set to $\underline{\theta}_{ij} = -45^{\circ}$ and $\overline{\theta}_{ij} = 45^{\circ}$, and are considered generous as in practice $\theta_{ij} = \overline{\theta}_{ij} - \underline{\theta}_{ij}$ typically does not exceed $\pm 10^{\circ}$ [60]. Furthermore, in the aim of improving the numerical conditioning of the problem, all the units are nondimensionalized with a base apparent power of $\SI{100}{\mega\VA}$, base flow rate of $\SI{100}{\meter\cubed\per\second}$, and base pressure of $10^6$$\SI{}{\pascal}$. Algorithm <ref> is initialized with $\epsilon = 10^{-4}$, $\alpha_{mn}^{t,(1)} = 0.1$, $\gamma = 10$, $\overline{\alpha} = 1000$, and $\sigma = 0.1{\max} \left( \left( {\max}\left( c_{2,gi},c_{1,gi}\right)  \right) _{gi \in \mathcal{G}} \right)$ for all test systems. In Algorithm <ref>, $\beta = 0.1 {\max} \left( \left(c_{sm} \right) _{sm \in \mathcal{S}} \right) $.

The below starts with a numerical assessment of Phase I of Algorithm <ref> in Section <ref>, followed by an assessment of Algorithm <ref> and Algorithm <ref> in Sections <ref> and <ref>. Two metrics are used for measuring optimality. The first is the optimality gap (%), defined as

\begin{align*}
Ogap= \left( f_{0}\left( x^{\star} \right) - f_{0}\left( x^{ \rm cvx} \right) \right) /f_{0}\left( x^{\star} \right) \times 100,
\end{align*}

where $x^{\star}$ is the solution obtained by a local solver such as IPOPT, Juniper or Algorithms <ref> and <ref>, and $x^{\rm cvx}$ is the solution obtained by an MICP relaxation such as the one in Problem <ref>.[Note that IPOPT, Juniper, and Algorithms <ref> and <ref> cannot guarantee global optimality in this case. IPOPT is a local primal-dual IPM that consists of a sequence of second-order approximations coordinated by a filter line search method whereas Phase I of Algorithms <ref> and <ref> consists of a sequence of first-order approximations coordinated by (<ref>), (<ref>), and the second and third terms in (<ref>).] The second is the relative optimality gap (%), defined as

\begin{align*}
ROgap= \left( f_{0}\left( x^{\dagger} \right) - f_{0}\left( x^{\star} \right) \right) /f_{0}\left( x^{\dagger} \right) \times 100,
\end{align*}

where $x^{\star}$ is the solution obtained by a local IPM solver (IPOPT in this case), and $x^{\rm \dagger}$ is the solution obtained by Phase I of Algorithm <ref>.

§.§ Evaluation of Phase I of Algorithm <ref>

Optimality and feasibility of Phase I of Algorithm <ref> compared to MOEGFc (solved using IPOPT [55]), and MICPc (solved using Gurobi [48]).

1c|     2c0pt7pt MICPc (Gurobi)     2|c|0pt7pt MOEGFc (IPOPT)     6c|Phase I (Gurobi)

1-1 4-11
1|c|Test     2c|     2c|Cold start     3c|Cold start     3c|Warm start (Problem <ref>)

2-3 4-11
0pt10pt case     Cost ($)     $C_{\rm max}^{\rm cvx}$     Cost ($)     $Ogap$ (%)     Cost ($)     $ROgap$ (%)     $C_{\rm mean}^{\dagger}$     Cost ($)     $ROgap$ (%)     $C_{\rm mean}^{\dagger}$
A	   	99946.5	   	394.67	   	99946.5	   	-1.26E-06	   	99946.5	   	-1.02E-06	   	1.86E-09	   	99946.5	   	-1.15E-06	   	2.16E-10
B	   	93999.6	   	195.00	   	94000.4	   	8.87E-04	   	94003.9	   	-3.70E-03	   	9.73E-08	   	94000.4	   	-3.94E-05	   	1.67E-07
C1	   	2975270.4	   	10.25	   	2975990.2	   	2.42E-02	   	2975988.5	   	5.73E-05	   	2.51E-09	   	2975988.2	   	6.79E-05	   	3.98E-10
C2	   	6315991.7	   	8.96	   	6318017.0	   	3.21E-02	   	6318445.6	   	-6.78E-03	   	2.90E-08	   	6318097.4	   	-1.27E-03	   	8.97E-09
C3	   	11358689.3	   	9.96	   	11359946.4	   	1.11E-02	   	11360188.9	   	-2.14E-03	   	3.78E-08	   	11359700.5	   	2.16E-03	   	1.42E-08

Computational effort of Phase I of Algorithm <ref> compared to MOEGFc (solved using IPOPT [55]), and MICPc (solved using Gurobi [48]). The numbers in parenthesis denote the number of iterations of Phase I.

1|c|Test     4c|CPU time (s)

0pt10pt Case     MICPc     IPOPT     Phase I (Cold)     Phase I (Warm)
A	   	8.1	   	1.5	   	0.2 (3)	   	0.3 (5)
B	   	2856.3	   	32.0	   	4.1 (13)	   	2.2 (7)
C1	   	8873.6	   	937.0	   	12.7 (10)	   	14.8 (12)
C2	   	10751.7	   	993.0	   	19.1 (11)	   	17.2 (11)
C3	   	12069.1	   	551.3	   	34.4 (11)	   	18.3 (10)

This section assesses the optimality, feasibility, and computational efficiency of Phase I of Algorithm <ref> (which is the same as Phase I of Algorithm <ref>), IPOPT on the continuous relaxation of Problem <ref>, and the MICP relaxation in (<ref>) with (<ref>) instead of (<ref>). For ease of exposition, these three approaches are called “Phase I”, “MOEGFc”, and “MICPc”, respectively. Two starting point strategies are considered in the initialization of Phase I and IPOPT (MOEGFc). The first is a cold start strategy that assumes $\theta_{ij}^{t,(1)} = 0.01$ for all $ij \in \mathcal{L}$, $\phi_{mn}^{t,(1)} = 0.1\overline{\phi}_{mn}$ for all $mn \in \mathcal{P}$, and $\wp_{m}^{t,(1)} = 0.5(\overline{\wp}_{m}+\underline{\wp}_{m})$ for all $m \in \mathcal{N}$, $t \in \mathcal{T}$. The warm start strategy consists of using the solution of Problem <ref> in the initialization of Phase I. Table <ref> compares the optimality and feasibility of Phase I to those of IPOPT on the continuous relaxation of Problem <ref> (MOEGFc), and Table <ref> compares the computational effort of the three.[Warm starting IPOPT did not noticeably affect the quality or the computational effort of the solution. Those results are therefore not shown in Table <ref>.] Table <ref> shows the objective function values of MICPc, MOEGFc, and Phase I under cold start, and Phase I under warm start in columns 2, 4, 6, and 9, respectively. Table <ref> also shows the maximum constraint violation of the solution of Problem <ref> with (<ref>) instead of (<ref>) (i.e., MICPc) in column 3 under $C_{\rm max}^{\rm cvx}$. These large constraint violations confirm the infeasibility of the MICPc solutions. Nonetheless, the MICPc can still provide a good-quality lower bound on the solution of MOEGFc as shown in column 5. More specifically, it can be inferred that the MOEGFc solution to test Case A is optimal, and the solution to test case C3 is within $1.65 \times 10^{-2}$% of the optimum. Additionally, since the largest $ROgap$ in columns 7 and 10 does not exceed $5 \times 10^{-3}$, and the largest constraint violation at the termination of Phase I is around $10^{-5}$ (i.e., $C_{\rm max}^{\dagger} \approx 10^{-5}$), it can be concluded that Phase I reaches near-optimal and feasible solutions to MOEGFc. The mean violations, $C_{\rm mean}^{\dagger} = {\rm mean}\left( \left( \left| h_{i}\left( x^{\dagger}\right) \right| \right)_{i=1,\ldots,p}\right)$, at the termination of Phase I are in the order of $10^{-8}$ on average for all test cases. Decreasing the constraint feasibility tolerance beyond $10^{-4}$ results in negligible change in the objective function and in the total linepack in all test cases.

The evolution of the maximum constraint violation $C_{\rm max}^{k}$ for test case C3 is shown in Fig. <ref>.

C$C_{\rm max}^{k}$
Victorian IEGSTest case C3

Phase I's $C_{\rm max}^{k}$ vs $k$ for test case C3 under a warm start.

Furthermore, Table <ref> shows that Phase I solves MOEGFc in computation times that are at least one order of magnitude smaller than those of IPOPT. On the other hand, the computational effort of solving MICPc is disproportionately larger than those of IPOPT and Phase I due to the presence of binary variables (one for each pipeline as in (<ref>)), which increases the complexity of the problem. Finally, although only two starting point strategies are shown in Tables <ref> and <ref>, the SLP algorithm underpinning Phase I is in fact robust to the choice of starting point $x^{(1)}$, owing to (i) the introduction slack variables to prevent infeasible problems, (ii) an automatic tuning of the convergence parameter $\alpha_{mn}^{t,(k)}$, and (iii) an iteratively refined polyhedral outer approximation of the nonconvex sets of the form $\mathcal{V}= \left\{x^{2} | x \in \left[\underline{x},\overline{x}\right]\right\}$ with the use of both supporting hyperplanes and halfspaces. The computational advantages on the warm start (using the solution of Problem <ref>) become more prominent in the solution of Problem <ref>, as will be shown in the next two sections.

§.§ Evaluation of Algorithm <ref>

Optimality and feasibility of Algorithm <ref> compared to MOEGF (solved using Juniper [57]), and MICP (solved using Gurobi [48]). “NA” designates instances that did not converge after one day of run time.

1c     2|c|0pt7pt MICP (Gurobi)     2|c|0pt7pt MOEGF (Juniper)     6c|Algorithm <ref> (Gurobi)

1-1 4-11
1|c|Test     2c|     2c|Cold start     3c|Cold start     3c|Warm start (Problem <ref>)

2-3 4-11
0pt10pt case     Cost ($)     $C_{\rm max}^{\rm cvx}$     Cost ($)     $Ogap$ (%)     Cost ($)     $Ogap$ (%)     $C_{\rm mean}^{\dagger}$     Cost ($)     $Ogap$ (%)     $C_{\rm mean}^{\dagger}$
A	   	99946.5	   	400.66	   	99946.5	   	-1.71E-06	   	99946.7	   	2.06E-04	   	3.41E-08	   	99946.3	   	-2.47E-04	   	2.93E-08
B	   	93999.8	   	189.17	   	NA	   	NA	   	94003.8	   	4.25E-03	   	2.60E-08	   	94001.6	   	1.88E-03	   	2.13E-08
C1	   	2975274.4	   	10.72	   	NA	   	NA	   	2975884.3	   	2.05E-02	   	7.70E-08	   	2975975.6	   	2.36E-02	   	9.91E-09
C2	   	6316166.5	   	8.16	   	NA	   	NA	   	6318099.4	   	3.06E-02	   	1.76E-08	   	6318009.4	   	2.92E-02	   	1.76E-08
C3	   	11358722.1	   	10.36	   	NA	   	NA	   	11360828.9	   	1.85E-02	   	5.72E-08	   	11360640.7	   	1.69E-02	   	3.08E-08

[Performance Alg_SMILP]Computational effort of Algorithm <ref> compared to Juniper [57], and MICP (solved using Gurobi [48]). The numbers in parenthesis denote the total number of iterations of  Algorithm <ref>. “DNC” designates instances that did not converge after one day of run time.

1|c|Test     4c|CPU time (s)

0pt10pt Case     MICP     Juniper     Alg. <ref> (Cold)     Alg. <ref> (Warm)
A	   	9.6	   	2.9	   	3.3 (6)	   	1.6 (7)
B	   	3986.2	   	DNC	   	14.4 (16)	   	11 (10)
C1	   	12750.4	   	DNC	   	59.8 (24)	   	40.5 (17)
C2	   	14867.6	   	DNC	   	77.8 (19)	   	65.3 (16)
C3	   	17205.7	   	DNC	   	536.7 (24)	   	305.6 (21)

The computation times of Algorithm <ref> with the warm start include the CPU times of Problem <ref>.
This section assesses the optimality, feasibility, and computational efficiency of Algorithm <ref>, Juniper on Problem <ref>, and the MICP relaxation in (<ref>). The starting point strategies adopted here are similar to the ones discussed in the previous section. The optimality and feasibility, and the computational effort of Algorithm <ref> compared to those of Juniper and the MICP relaxation in (<ref>) are shown in Tables <ref> and <ref>, respectively. It is evident from those two tables that the original MINLP problem in (<ref>) is extremely challenging to solve using an NLBB solver such as Juniper. In fact, Juniper did not converge (after one day of run time) on the practical-size systems in test cases B, C1, C2, C3. As a result, Table <ref> only shows the optimality gap, as opposed to the relative optimality gap, of the solution of Algorithm <ref> as it is the only available feasible solution to Problem <ref> for test cases B, C1, C2, C3. The MICP relaxation in (<ref>) is also more challenging to solve due the additional number of binary variables associated with the non-pipe elements (constraint (<ref>)). On test case C3, the MICP relaxation now takes more than 4 hours to converge. On the other hand, Algorithm <ref> with the warm start strategy takes 5 minutes. Not only Algorithm <ref> is substantially faster than the MICP relaxation in (<ref>), it also converges to near-optimal and feasible solutions ($C_{\rm max}^{\dagger} \approx 10^{-5}$) with an average constraint violation in the order of $10^{-8}$. The solution to test case C3 is within $2.67 \times 10^{-2}$% of the optimum. In comparison, the MICP relaxation converges to infeasible solutions (with large constraint violations $C_{\rm max}^{\rm cvx}$) in substantially longer computation times. The evolution of the maximum constraint violation $C_{\rm max}^{k}$ for test case C3 is shown in Fig. <ref>.

C$C_{\rm max}^{k}$
Victorian IEGSTest case C3
Phase IPhase I
Phase IIPhase II

Algorithm <ref>'s $C_{\rm max}^{k}$ vs $k$ for test case C3 under a warm start.

Finally, the computational advantage of warm starting Algorithm <ref> from the solution of Problem <ref> as compared to a cold start is more salient in this case, as the last column of Table <ref> shows a noticeable improvement across all the test cases. In fact, Algorithm <ref> with the warm start strategy is at least one order of magnitude faster than Juniper and the MICP relaxtion in (<ref>).

§.§ Evaluation of Algorithm <ref>

Optimality and feasibility of Algorithm <ref> compared to MOEGF (solved using Juniper [57]), and MICP (solved using Gurobi [48]). “NA” designates instances that did not converge after one day of run time.

1c     2|c|0pt7pt MICP (Gurobi)     2|c|0pt7pt MOEGF (Juniper)     6c|Algorithm <ref> (Gurobi)

1-1 4-11
1|c|Test     2c|     2c|Cold start     3c|Cold start     3c|Warm start (Problem <ref>)

2-3 4-11
0pt10pt case     Cost ($)     $C_{\rm max}^{\rm cvx}$     Cost ($)     $Ogap$ (%)     Cost ($)     $Ogap$ (%)     $C_{\rm mean}^{\dagger}$     Cost ($)     $Ogap$ (%)     $C_{\rm mean}^{\dagger}$
A	   	99946.5	   	400.66	   	99946.5	   	-1.71E-06	   	99948.8	   	2.30E-03	   	2.71E-07	   	99948.8	   	2.29E-03	   	2.40E-07
B	   	93999.8	   	189.17	   	NA	   	NA	   	94004.7	   	5.22E-03	   	7.41E-08	   	94001.9	   	2.24E-03	   	2.23E-08
C1	   	2975274.4	   	10.72	   	NA	   	NA	   	2979358.7	   	1.37E-01	   	1.63E-08	   	2979435.8	   	1.40E-01	   	4.20E-08
C2	   	6316166.5	   	8.16	   	NA	   	NA	   	6322072.4	   	9.34E-02	   	1.50E-07	   	6321209.9	   	7.98E-02	   	4.31E-08
C3	   	11358722.1	   	10.36	   	NA	   	NA	   	11364530.9	   	5.11E-02	   	3.56E-08	   	11364465.6	   	5.05E-02	   	2.27E-08

[Performance Alg_SLPH]Computational effort of Algorithm <ref> compared to Juniper [57], and MICP (solved using Gurobi [48]). The numbers in parenthesis denote the total number of iterations of Algorithm <ref>. “DNC” designates instances that did not converge after one day of run time.

1|c|Test     4c|CPU time (s)

0pt10pt Case     MICP     Juniper     Alg. <ref> (Cold)     Alg. <ref> (Warm)
A	   	9.6	   	2.9	   	0.4 (6)	   	0.5 (8)
B	   	3986.2	   	DNC	   	4.9 (16)	   	2.9 (10)
C1	   	12750.4	   	DNC	   	23.3 (19)	   	21.1 (21)
C2	   	14867.6	   	DNC	   	37.2 (17)	   	30.7 (17)
C3	   	17205.7	   	DNC	   	71.3 (20)	   	37.3 (18)

The computation times of Algorithm <ref> with the warm start include the CPU times of Problem <ref>.
The optimality and feasibility, and the computational effort of Algorithm <ref> compared to those of Juniper (on Problem <ref>) and the MICP relaxation in (<ref>) are shown in Tables <ref> and <ref>, respectively. In this case, the optimality gaps on the Victorian test cases C1, C2, and C3 are still less than 0.2%, which corroborates the high quality of the solutions. The solutions are also feasible to within at most $1.8 \times 10^{-5}$ (i.e., $C_{\rm max}^{\dagger} \approx 1.8 \times 10^{-5}$) with an average constraint violation in the order of $10^{-8}$. The evolution of the maximum constraint violation $C_{\rm max}^{k}$ for test case C3 is shown in Fig. <ref>.

C$C_{\rm max}^{k}$
Victorian IEGSTest case C3
Phase IPhase I
Phase IIPhase II

Algorithm <ref>'s $C_{\rm max}^{k}$ vs $k$ for test case C3 under a warm start.

Although it converges to slightly higher objective function values compared to Algorithm <ref>, the LP-based method in Algorithm <ref> is substantially faster. In particular, the speed-up is more than 8$\times$ on test case C3 under the same warm start strategy. And once again, the warm start strategy is generally noticeably faster than the cold start strategy, which underscores the value of using the solution of Problem <ref> as a starting point. Finally, it is worth noting that both Algorithms <ref> and <ref> are exact only with respect to Problem <ref>, which incorporates a quasi-dynamic gas flow model, and not with respect to the original partial differential equations (PDEs) describing the full dynamic flow of gas.

§.§ Linepack backtesting

To assess the validity of the solutions of Test cases C1, C2, and C3 from a physical gas network standpoint, they are backtested against actual 24-hour total linepack profiles for the state of Victoria, Australia. The 24-hour total linepack profiles of both Algorithm <ref> and the actual one from (AEMO) [42] are shown in Fig. <ref> for the high demand day (09/08/2019).

t$t$ ($\SI{}{\hour}$)
LinepackLinepack ($\SI{}{\tera\joule}$)
ActualActual (AEMO)
AlgorithmAlgorithm <ref>

Linepack backtesting for Test case C3.

The maximum and mean errors between the two are $1.21\%$ and $0.38\%$, respectively, which could be attributed to the fact that the actual linepack on that high-demand day was not a result of optimized operation. The linepack backtesting not only serves as a testament to the validity of the solution of Algorithm <ref>, but also to the gas network model of the Victorian DTS which was developed from scratch as part of this work.

§ CONCLUSION

This paper introduced two novel SLP-based algorithms for efficiently solving the MOEGF problem. The first is an iterative MILP-based algorithm and the second is an iterative LP-based algorithm with an elaborate procedure for ensuring an integral solution. Numerical evaluation demonstrates that both algorithms can solve the MOEGF problem to high-quality feasible solutions in computation times at least two orders of magnitude faster than both a state-of-the-art NLBB MINLP solver and an MICP relaxation, on the real electricity and gas transmission networks of the state of Victoria with actual linepack and demand profiles. Moreover, both approaches are warm-started from the solution of a novel polyhedral relaxation of the problem for a noticeable improvement in computation time compared to a cold start. Finally, while not claiming it is superior to existing MINLP solvers, the proposed iterative LP-based method represents a fast alternative for the specific type of problem addressed here, especially in settings where only LP solvers can be used (as is the case with existing dispatch engines of ISO like AEMO). To do so, the iterative LP-based method exploits the structure and engineering properties of the MOEGF problem and tailors a novel SLP approach, compounded by a fast heuristic that ensures an integral solution.

§ ACKNOWLEDGMENT

This work is supported by Future Fuels CRC as part of the R.P1.1-02: “Regional Case Studies on Multi-Energy System Integration” project. The cash and in-kind support from the industry participants is gratefully acknowledged.

[1]
M. Geidl and G. Andersson, “Optimal Power Flow of Multiple Energy
Carriers,” IEEE Transactions on Power Systems, vol. 22, no. 1, pp.
145–155, 2007.

[2]
J. Qiu, H. Yang, Z. Y. Dong, J. H. Zhao, K. Meng, F. J. Luo, and K. P. Wong,
“A Linear Programming Approach to Expansion Co-Planning in Gas and
Electricity Markets,” IEEE Transactions on Power Systems, vol. 31,
no. 5, pp. 3594–3606, 2016.

[3]
S. Luo, L. Yang, X. Zhang, W. Chen, K. Wang, and Y. Xu, “A Fully
Linear-Constrained Optimal Electricity-Gas Flow in an Integrated Energy
System,” in 2nd IEEE Conference on Energy Internet and Energy System
Integration, EI2 2018 - Proceedings.1em plus 0.5em minus 0.4emIEEE, oct 2018, pp. 1–6.

[4]
X. Fang, M. T. Craig, and B.-M. Hodge, “Linear Approximation Line Pack Model
for Integrated Electricity and Natural Gas Systems OPF,” IEEE PES
General Meeting Atlanta, no. August, 2019.

[5]
C. Wang, W. Wei, J. Wang, L. Bai, Y. Liang, and T. Bi, “Convex Optimization
Based Distributed Optimal Gas-Power Flow Calculation,” IEEE
Transactions on Sustainable Energy, vol. 9, no. 3, pp. 1145–1156, 2018.

[6]
S. Chen, A. J. Conejo, R. Sioshansi, and Z. Wei, “Unit Commitment with an
Enhanced Natural Gas-Flow Model,” IEEE Transactions on Power
Systems, vol. 34, no. 5, pp. 3729–3738, 2019.

[7]
L. Yang, Y. Xu, S. Hongbin, and X. Zhao, “Two-Stage Convexification-Based
Optimal Electricity-Gas Flow,” IEEE Transactions on Smart Grid,
vol. 11, no. 2, pp. 1465–1475, 2020.

[8]
S. D. Manshadi and M. E. Khodayar, “A Tight Convex Relaxation for the Natural
Gas Operation Problem,” IEEE Transactions on Smart Grid, vol. 9,
no. 5, pp. 5467–5469, 2018.

[9]
C. M. Correa-Posada and P. Sánchez-Martín, “Gas Network
Optimization: A comparison of Piecewise Linear Models,” Optimization
online, pp. 133–159, 2014.

[10]
——, “Security-constrained unit commitment with dynamic gas
constraints,” IEEE Power and Energy Society General Meeting, vol.
2015-Septe, pp. 1–5, 2015.

[11]
X. Zhang, M. Shahidehpour, A. Alabdulwahab, and A. Abusorrah, “Hourly
Electricity Demand Response in the Stochastic Day-Ahead Scheduling of
Coordinated Electricity and Natural Gas Networks,” IEEE Transactions
on Power Systems, vol. 31, no. 1, pp. 592–601, 2016.

[12]
C. He, L. Wu, T. Liu, and M. Shahidehpour, “Robust Co-Optimization Scheduling
of Electricity and Natural Gas Systems via ADMM,” IEEE Transactions
on Sustainable Energy, vol. 8, no. 2, pp. 658–670, 2017.

[13]
C. Shao, X. Wang, M. Shahidehpour, X. Wang, and B. Wang, “An MILP-Based
Optimal Power Flow in Multicarrier Energy Systems,” IEEE Transactions
on Sustainable Energy, vol. 8, no. 1, pp. 239–248, 2017.

[14]
S. Chen, Z. Wei, G. Sun, Y. Sun, H. Zang, and Y. Zhu, “Optimal Power and Gas
Flow with a Limited Number of Control Actions,” IEEE Transactions on
Smart Grid, vol. 9, no. 5, pp. 5371–5380, 2018.

[15]
Y. Zhang, Y. Hu, J. Ma, and Z. Bie, “A Mixed-Integer Linear Programming
Approach to Planning of Natural Gas and Electricity,” IEEE
Transactions on Power Systems, vol. 33, no. 6, pp. 6368–6378, 2018.

[16]
C. He, X. Zhang, T. Liu, and L. Wu, “Distributionally Robust Scheduling of
Integrated Gas-Electricity Systems with Demand Response,” IEEE
Transactions on Power Systems, vol. 34, no. 5, pp. 3791–3803, 2019.

[17]
C. Borraz-Sánchez, R. Bent, S. Backhaus, H. Hijazi, and P. Van
Hentenryck, “Convex Relaxations for Gas Expansion Planning,”
INFORMS Journal on Computing, vol. 28, no. 4, pp. 645–656, 2016.

[18]
Y. Wen, X. Qu, W. Li, X. Liu, and X. Ye, “Synergistic operation of
electricity and natural gas networks via ADMM,” IEEE Transactions on
Smart Grid, vol. 9, no. 5, pp. 4555–4565, 2018.

[19]
Y. He, M. Yan, M. Shahidehpour, Z. Li, C. Guo, L. Wu, and Y. Ding,
“Decentralized Optimization of Multi-Area Electricity-Natural Gas Flows
Based on Cone Reformulation,” IEEE Transactions on Power Systems,
vol. 33, no. 4, pp. 4531–4542, 2018.

[20]
M. K. Singh and V. Kekatos, “Natural gas flow equations: Uniqueness and an
MI-SOCP solver,” Proceedings of the American Control Conference,
vol. 2019-July, pp. 2114–2120, 2019.

[21]
A. Schwele, C. Ordoudis, J. Kazempour, and P. Pinson, “Coordination of Power
and Natural Gas Systems: Convexification Approaches for Linepack Modeling,”
2019 IEEE Milan PowerTech, pp. 1–6, 2019.

[22]
L. Yang, X. Zhao, and Y. Xu, “A convex optimization and iterative solution
based method for optimal power-gas flow considering power and gas losses,”
International Journal of Electrical Power and Energy Systems, vol.
121, no. December 2019, p. 106023, 2020.

[23]
A. R. Sayed, C. Wang, T. Bi, M. A. Abdelbaky, and A. Masood, “Optimal
Power-Gas Flow of Integrated Electricity and Natural Gas System: A Sequential
MISOCP Approach,” 2019 IEEE 3rd Conference on Energy Internet and
Energy System Integration (EI2), pp. 283–288, 2020.

[24]
L. Yang, Y. Xu, W. Gu, and H. Sun, “Distributionally Robust
Chance-constrained Optimal Power-Gas Flow under Bidirectional Interactions
Considering Uncertain Wind Power,” IEEE Transactions on Smart Grid,
vol. 3053, no. c, pp. 1–1, 2020.

[25]
S. Clegg and P. Mancarella, “Integrated Electrical and Gas Network
Flexibility Assessment in Low-Carbon Multi-Energy Systems,” IEEE
Transactions on Sustainable Energy, vol. 7, no. 2, pp. 718–731, 2016.

[26]
J. Yang, N. Zhang, C. Kang, and Q. Xia, “Effect of Natural Gas Flow Dynamics
in Robust Generation Scheduling under Wind Uncertainty,” IEEE
Transactions on Power Systems, vol. 33, no. 2, pp. 2087–2097, 2018.

[27]
S. Clegg and P. Mancarella, “Integrated Modeling and Assessment of the
Operational Impact of Power-to-Gas (P2G) on Electrical and Gas Transmission
Networks,” IEEE Transactions on Sustainable Energy, vol. 6, no. 4,
pp. 1234–1244, oct 2015.

[28]
——, “Integrated electricity-heat-gas modelling and assessment, with
applications to the Great Britain system. Part II: Transmission network
analysis and low carbon technology and resilience case studies,”
Energy, 2018.

[29]
I. Saedi, S. Mhanna, H. Wang, and P. Mancarella, “Integrated Electricity and
Gas Systems Modelling : Assessing the Impacts of Electrification of
Residential Heating in Victoria,” in Australasian Universities Power
Engineering Conference, no. December, 2020, pp. 1–6.

[30]
J. Huchette and J. P. Vielma, “Nonconvex piecewise linear functions: Advanced
formulations and simple modeling tools,” arXiv, pp. 1–40, 2017.

[31]
“IEGS data.” [Online]. Available:

[32]
T. Lipp and S. Boyd, “Variations and extension of the convex–concave
procedure,” Optimization and Engineering, vol. 17, no. 2, pp.
263–287, 2016.

[33]
C. Liu, M. Shahidehpour, Y. Fu, and Z. Li, “Security-constrained unit
commitment with natural gas transmission constraints,” IEEE
Transactions on Power Systems, vol. 24, no. 3, pp. 1523–1536, aug 2009.

[34]
A. Alabdulwahab, A. Abusorrah, X. Zhang, and M. Shahidehpour, “Stochastic
security-constrained scheduling of coordinated electricity and natural gas
infrastructures,” IEEE Systems Journal, vol. 11, no. 3, pp.
1674–1683, 2017.

[35]
B. Zhao, A. J. Conejo, and R. Sioshansi, “Coordinated Expansion Planning of
Natural Gas and Electric Power Systems,” IEEE Transactions on Power
Systems, vol. 33, no. 3, pp. 3064–3075, 2018.

[36]
C. He, L. Wu, T. Liu, and Z. Bie, “Robust Co-Optimization Planning of
Interdependent Electricity and Natural Gas Systems with a Joint N-1 and
Probabilistic Reliability Criterion,” IEEE Transactions on Power
Systems, vol. 33, no. 2, pp. 2140–2154, 2018.

[37]
C. Wang, R. Gao, W. Wei, M. Shafie-Khah, T. Bi, and J. P. Catalão,
“Risk-based distributionally robust optimal gas-power flow with wasserstein
distance,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp.
2190–2204, 2019.

[38]
C. M. Correa-Posada and P. Sánchez-Martín, “Integrated Power and
Natural Gas Model for Energy Adequacy in Short-Term Operation,” IEEE
Transactions on Power Systems, vol. 30, no. 6, pp. 3347–3355, 2015.

[39]
M. Chaudry, N. Jenkins, and G. Strbac, “Multi-time period combined gas and
electricity network optimisation,” Electric Power Systems Research,
vol. 78, no. 7, pp. 1265–1279, jul 2008.

[40]
A. Castillo, P. Lipka, J.-P. Watson, S. S. Oren, and R. P. O'Neill, “A
successive linear programming approach to solving the IV-ACOPF,”
IEEE Transactions on Power Systems, vol. 31, no. 4, pp. 2752–2763,

[41]
R. P. O'Neill, M. Williard, B. Wilkins, and R. Pike, “Mathematical
Programming Model for Allocation of Natural Gas.” Operations
Research, vol. 27, no. 5, pp. 857–873, 1979.

[42]
“AEMO.” [Online]. Available:

[43]
H. Ameli, M. Qadrdan, and G. Strbac, “Coordinated operation strategies for
natural gas and power systems in presence of gas‐related flexibilities,”
IET Energy Systems Integration, vol. 1, no. 1, pp. 3–13, 2019.

[44]
A. J. Osiadacz, Simulation and analysis of gas networks.1em
plus 0.5em minus 0.4emGulf Pub Co, 1987.

[45]
E. S. Menon, Gas Pipeline Hydraulics.1em plus 0.5em minus
0.4emTaylor & Francis Inc, 2005.

[46]
C. Borraz-Sánchez and R. Z. Ríos-Mercado, “Improving the
operation of pipeline systems on cyclic structures by tabu search,”
Computers and Chemical Engineering, vol. 33, no. 1, pp. 58–64, 2009.

[47]
M. A. Duran and I. E. Grossmann, “An outer-approximation algorithm for a
class of mixed-integer nonlinear programs,” Mathematical
Programming, vol. 36, no. 3, pp. 307–339, 1986.

[48]
“Gurobi Optimizer Reference Manual, Version 5.0,” www.gurobi.com,
pp. 1–786, 2019.

[49]
S. Mhanna, G. Verbič, and A. C. Chapman, “Tight LP Approximations for
the Optimal Power Flow Problem,” 19th Power Systems Computation
Conference, PSCC 2016, 2016.

[50]
E. W. Cheney and A. A. Goldstein, “Newton's method for convex programming and
Tchebycheff approximation,” Numerische Mathematik, vol. 1, no. 1,
pp. 253–268, 1959.

[51]
J. E. Kelley, Jr., “The Cutting-Plane Method for Solving Convex
Programs,” Journal of the Society for Industrial and Applied
Mathematics, vol. 8, no. 4, pp. 703–712, 1960.

[52]
R. H. Byrd, J. Nocedal, and R. A. Waltz, “Knitro: An integrated package for
nonlinear optimization,” in Large Scale Nonlinear Optimization,
35–59, 2006.1em plus 0.5em minus 0.4emSpringer Verlag, 2006,
pp. 35–59.

[53]
J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh
approach to numerical computing,” SIAM review, vol. 59, no. 1, pp.
65–98, 2017.

[54]
I. Dunning, J. Huchette, and M. Lubin, “Jump: A modeling language for
mathematical optimization,” SIAM Review, vol. 59, no. 2, pp.
295–320, 2017.

[55]
A. Wächter and L. T. Biegler, “On the implementation of an interior-point
filter line-search algorithm for large-scale nonlinear programming,”
Mathematical Programming, vol. 106, no. 1, pp. 25–57, 2006.

[56]
I. S. Duff, “MA57—a code for the solution of sparse symmetric definite and
indefinite systems,” ACM Transactions on Mathematical Software,
vol. 30, no. 2, pp. 118–144, jun 2004.

[57]
O. Kröger, C. Coffrin, H. Hijazi, and H. Nagarajan, “Juniper: An open-source
nonlinear branch-and-bound solver in julia,” in Integration of
Constraint Programming, Artificial Intelligence, and Operations
Research.1em plus 0.5em minus 0.4emSpringer International
Publishing, 2018, pp. 377–386.

[58]
“Future Fuels CRC.” [Online]. Available:

[59]
A. Xenophon and D. Hill, “Open grid model of Australia's national electricity
market allowing backtesting against historic data,” Scientific Data,
vol. 5, pp. 1–21, 2018.

[60]
K. Purchala, L. Meeus, D. Van Dommelen, and R. Belmans, “Usefulness of DC
power flow for active power flow analysis,” in Power Engineering
Society General Meeting, 2005. IEEE, June 2005, pp. 454–459 Vol. 1.

[]Sleiman Mhanna
(S'13–M'16) received the B.Eng. degree (with high distinction) from the Notre Dame University, Lebanon, and the M.Eng. degree from the American University of Beirut, Lebanon, in 2010 and 2012, respectively, both in electrical engineering. He received the Ph.D. degree from the School of Electrical and Information Engineering, Centre for Future Energy Networks, University of Sydney, Australia, in 2016. He is currently a Research Fellow at the Department of Electrical and Electronic Engineering, The University of Melbourne, Australia. His research interests include computational methods for integrated multi-energy systems, decomposition methods, and demand response.

[]Isam Saedi (S'20) obtained the master’s degree in electrical power systems engineering from The University of Manchester, U.K., in 2015. He is currently in pursuit of his Ph.D. degree at the Department of Electrical and Electronic Engineering, The University of Melbourne, Australia. His research focuses on the assessment of integrated electricity-gas-hydrogen systems in the presence of different coupling technologies and scenarios for different sectors.

[]Pierluigi Mancarella
(SM) is Chair Professor of Electrical Power Systems at The University of Melbourne, Melbourne, Australia, and Professor of Smart Energy Systems at The University of Manchester, Manchester, UK. His research interests include techno-economic modelling of integrated multi-energy systems; security, reliability and resilience of future networks; and energy infrastructure planning under uncertainty. Pierluigi is an Editor of the IEEE Transactions on Power Systems, IEEE Transactions on Smart Grid, and IEEE Systems Journal, and an IEEE Power and Energy Society Distinguished Lecturer.