Convex Relaxations of Maximal Load Delivery for Multi-contingency Analysis of Joint Electric Power and Natural Gas Transmission Networks

Recent increases in gas-fired power generation have engendered increased interdependencies between natural gas and power transmission systems. These interdependencies have amplified existing vulnerabilities to gas and power grids, where disruptions can require the curtailment of load in one or both systems. Although typically operated independently, coordination of these systems during severe disruptions can allow for targeted delivery to lifeline services, including gas delivery for residential heating and power delivery for critical facilities. To address the challenge of estimating maximum joint network capacities under such disruptions, we consider the task of determining feasible steady-state operating points for severely damaged systems while ensuring the maximal delivery of gas and power loads simultaneously, represented mathematically as the nonconvex joint Maximal Load Delivery (MLD) problem. To increase its tractability, we present a mixed-integer convex relaxation of the MLD problem. Then, to demonstrate the relaxation's effectiveness in determining bounds on network capacities, exact and relaxed MLD formulations are compared across various multi-contingency scenarios on nine joint networks ranging in size from 25 to 1,191 nodes. The relaxation-based methodology is observed to accurately and efficiently estimate the impacts of severe joint network disruptions, often converging to the relaxed MLD problem's globally optimal solution within ten seconds.

Components directed to i ∈ J Natural Gas Parameters f ij , f ij ∈ R Lower, upper mass flow bounds for (i, j) ∈ A p i , p i ≥ 0 Lower, upper pressure bounds for i ∈ J s i ≥ 0 Upper supply mass flow bound for i ∈ R d i ≥ 0 Upper demand mass flow bound for i ∈ D w ij ≥ 0 Resistance of pipe (i, j) ∈ P τ ij ≥ 0 Resistance of resistor (i, j) ∈ T α ij , α ij ≥ 0 Lower, upper scalars for (i, j) ∈ W ∪ C Natural Gas Variables Demand at delivery i ∈ D p i ≥ 0 Pressure at junction i ∈ J z ij ∈ {0, 1} Status of controllable element (i, j) ∈ V ∪ W Interdependency Modeling K Set of links between i ∈ D and j ∈ G h 1,2,3 ij Coefficients of heat rate curve for (i, j) ∈ K

I. INTRODUCTION
B ETWEEN 2012 and 2040, global electric power genera- tion capacity is predicted to increase from 21.6 million gigawatt-hours (GWh) to 36.5 million GWh.Of this, gasfired generation is expected to increase from 22% to 28%

Gas
Power load delivered (%) Fig. 1.A high-level illustration of natural gas and power transmission network responses to a hypothetical severe disruption.The shaded region indicates the points in the disruption and restoration timeline that are studied in this paper using an optimization-based assessment of damaged network capacities.[1].This growing dependence underscores the increasing sensitivity of power systems to upstream disruptions in gas pipelines.The most recent example is the February 2021 Texas power crisis, where the Electric Reliability Council of Texas experienced a loss of nearly 52.3 GW (48.6%) of its generation capacity.Nearly half of the loss was attributed to a lack of gas-fired power generation [2].Other examples include the 2014 polar vortex, where curtailments in gas delivery resulted in roughly 25% of generation outages throughout the Pennsylvania-New Jersey-Maryland Interconnection [3].Disruptions in the gas grid can also inhibit the transport of fuel required for residential heating.This begets an important tradeoff between gas delivery and power delivery during severe network disruptions.Understanding these interdependencies is critical for the resilience of gas and power delivery systems.
The contingency response measures considered in this paper are illustrated in Figure 1.Given a severe disruption, (i) gas and/or power load deliveries decrease as gas and/or power network elements are impaired and effects begin to cascade, and (ii) cascading effects subside, and a new stable operating point is realized.After (ii), load can gradually be restored via operational methods until (iii) network repairs begin.These restorative actions are performed until (iv) all gas and power loads can be delivered.Repairs continue until (v) all gas and power network components are again operational.Addressing all event types within Figure 1, however, is a substantial task.To make the scope more manageable, we focus more narrowly on ascertaining optimal steady-state operating points between events of types (ii) and (iv), i.e., decisions that maximize gas and power load delivery in the surviving gas-power system.
In this paper, this task is formalized as the steady-state joint Maximal Load Delivery (MLD) problem.The problem is informally stated as follows: given severely damaged gas and power networks in which multiple components have become nonoperational, maximize the amounts of prioritized gas and active power loads that can be served simultaneously in the damaged joint network, subject to steady-state natural gas and alternating current (AC) power network physics.The nonconvex physics and discrete nature of operations in the joint network (e.g., the opening and closing of valves in the gas network) render this a challenging mixed-integer nonlinear program (MINLP).To increase its tractability, we develop a mixed-integer convex programming (MICP) relaxation of the MLD problem.The MICP is found to be an effective means for bounding maximum total deliverable gas and power loads.
This paper expands upon existing MLD methods for independent gas [4] and power networks [5], as well as approaches from joint network modeling, to formulate and solve the joint gas-power MLD problem.Its contributions include • The first formulation of the gas-power MLD problem; • A reliable MICP relaxation of the MLD problem; • Proof-of-concept analyses of MLD gas-power tradeoffs.The remainder of this paper proceeds as follows: Section II reviews relevant gas, power, and joint steady-state optimization models that appear in the literature, then formulates the requirements for AC power and gas pipeline operational feasibility as an MINLP; Section III formulates the MLD problem as an MINLP, then proposes an MICP relaxation; Section IV rigorously benchmarks the MINLP and MICP formulations across multiple joint gas-power networks of various sizes, then provides proofs of concept for joint multi-contingency analysis using the MLD method; and Section V concludes the paper.

II. BACKGROUND FOR NETWORK MODELING
The past decade has seen remarkable theoretical and algorithmic advances in the independent fields of power and natural gas network optimization.A recent survey of relaxations and approximations used in power system optimization is presented by [6].The study in power most related to this paper is by [5], who introduce the AC MLD problem and propose various relaxations to increase its tractability.The method was later extended by [7] to identify the k components that maximize network disruption, as well as by [8] to identify multicontingency scenarios that would benefit from more detailed cascading analyses.The MLD problem was also exploited by [9], who applied it within a bilevel optimization for balancing wildfire risk and power outages.Finally, an implementation of the power MLD problem is presented by [10], who provide formulations via the POWERMODELSRESTORATION package.Their implementation, in fact, serves as a computational foundation for the power system modeling portion of this paper.
As with power, the growing utilization of gas networks has led to a variety of optimization studies.A summary of recent work related to the optimization-based assessment of gas network capacities is provided by [11].Steady-state models and approximations of gas network components amenable to optimization applications are provided by [12].However, the study in gas most related to this paper is by [4], who develop the steady-state gas MLD problem and an MICP relaxation.Another related study is by [13], who examine the the problem of identifying the k components of a gas network whose simultaneous failure maximizes disruption to the network.
An even more recent body of literature has examined the optimal coordination of gas and power infrastructures.A review of joint gas and power planning is given by [14].Other studies have focused on market coordination and energy pricing problems [15], [16].Many studies have assumed the

Model 1 Power Network Modeling Requirements
networks to be fully coordinated, examining optimal scheduling of generator dispatching and gas compressor operations [17].Recent studies have expanded upon these earlier joint "optimal gas-power flow" problems, developing specialized formulations and algorithms for related applications [18]- [20].
A smaller number of studies have considered joint problems related to restoration, e.g., scheduling of general large-scale interdependent infrastructures in [21] and last-mile restoration of joint gas and power systems in [22].The remaining subsections build upon previous studies to define the requirements for steady-state operation of a damaged joint gas-power network.

A. Power Transmission Network Modeling
a) Notation for Sets: A power network is represented by an arbitrarily directed graph (N , E ∪E R ), where N is the set of buses, E is the set of forward-directed branches (or lines), and E R is the set of branches in their reverse orientation.The set of generators (producers), loads (consumers), and shunts are denoted by G, L, and H, respectively, which are attached to existing buses i ∈ N .We let the subset of these components attached to i ∈ N be denoted by G i , L i , and H i .We next define the decision variables and constraints required to model a damaged AC power network's steady-state operations.
b) Power Network Modeling Requirements: The MINLP formulation for AC power network feasibility, as defined for AC MLD analysis, is presented in Model 1 and detailed by [5].Here, Constraints (1a) and (1b) model Ohm's law for lines, where S ij ∈ C denotes the variable power along each line; Y ij ∈ C and Y c ij ∈ C are constants denoting the line admittance and line charging; V i ∈ C denotes the variable voltage at bus i ∈ N ; and T ij ∈ C denotes constant transformer properties.Constraints (1c) model power balances from Kirchhoff's current law for each bus, where S g k ∈ C denotes the variable power supplied by generator k ∈ G; S d k ∈ C denotes the maximum power that can be delivered at load k ∈ L; and Y s k denotes the admittance of bus shunt k ∈ H.Note that z d k , k ∈ L i allows each load to vary between zero and its predefined maximum, and z s k allows for shedding fixed bus shunts from the network.These modifications ensure that power balance constraints are satisfied in damaged networks.
Constraints (1d)-(1i) impose engineering limits and variable bounds.Constraints (1d) bound the apparent power flow on each line, representing thermal limits.Constraints (1e) ensure that each voltage phase angle difference is limited by predefined lower and upper bounds, θ ∆ ij and θ ∆ ij , respectively.Constraints (1f) bound the voltage magnitude at each bus, where V i and V i denote lower and upper bounds, respectively.Here, z v i ∈ {0, 1} is a discrete variable that allows each bus to become de-energized when isolated from load or generation.Similarly, Constraints (1g) bound power generation, where S g i and S g i denote lower and upper bounds, respectively, and z g i ∈ {0, 1} allows for each generator to become uncommitted when required to satisfy Constraints (1a) and (1b).

B. Natural Gas Transmission Network Modeling a) Notation for Sets:
A gas pipeline network is modeled using a directed graph (J , A), where J is the set of nodes (i.e., junctions) and A is the set of components that connect two nodes.The sets of receipts (producers) and deliveries (consumers) are denoted by R and D, respectively.These components are considered to be attached to junctions i ∈ J .The subset of receipts attached to i ∈ J is denoted by R i and the subset of deliveries by D i .The sets of horizontal and short pipes are denoted by P ⊂ A and S ⊂ A, respectively; the set of resistors by T ⊂ A; the set of valves and pressurereducing regulators by V ⊂ A and W ⊂ A, respectively; and the set of compressors by C ⊂ A. Additionally, the set of node-connecting components incident to i ∈ J where i is the tail (respectively, head) of the arc is denoted by δ + i := {(i, j) ∈ A} (respectively, δ − i := {(j, i) ∈ A}).We next define the decision variables and constraints required to model a damaged gas network's steady-state operations.
b) Gas Network Modeling Requirements: The MINLP formulation for gas network feasibility, as defined for MLD analysis, is presented in Model 2 and detailed by [4].First, Constraints (2a) model nodal physics, i.e., mass flow conservation at junctions i ∈ J .Here, f ij ∈ R denotes the variable mass flow along each node-connecting component; s k ∈ R + denotes the variable supply at receipt k ∈ R; and d k ∈ R + denotes the variable demand (or load) at delivery k ∈ D.
Constraints (2b)-(2u) model the physics of node-connecting components.Constraints (2b) model the Weymouth relationship for steady-state flow in a gas pipeline for each horizontal pipe (i, j) ∈ P. Here, p i ∈ R + denotes the variable pressure at junction i ∈ J , and w ij ∈ R + denotes the constant mass flow resistance of the pipe.These constraints are the most frequent sources of nonconvex nonlinearity in modeling the gas system.
Constraints (2c) model short pipes in the network, which provide resistanceless mass transport between two junctions.Constraints (2d) model resistors in the network, which act as surrogate components capable of modeling pressure losses elsewhere from pipes.Here, pressure loss is modeled according to the Darcy-Weisbach equation, where τ ij ∈ R + is the resistance, which is a function of the resistor's unitless drag Model 2 Gas Network Modeling Requirements factor and (possibly artificial) diameter.Note that like Constraints (2b), these constraints are also nonconvex nonlinear.Constraints (2e)-(2g) model valves in the network.Here, the operating status of each valve (i, j) ∈ V is modeled using a discrete variable z ij ∈ {0, 1}, where z ij = 1 indicates an open valve and z ij = 0 indicates a closed valve.Constraints (2e) prohibit flow across each valve when z ij = 0. Constraints (2f) and (2g) model, when a valve is open, that the pressures at connecting junctions are equal.They also model the decoupling of junction pressures when the valve is closed.
Constraints (2h)-(2k) model regulators (i.e., pressurereducing valves) in the network.Similar to valves, the status of each regulator is modeled using a discrete variable z ij ∈ {0, 1}, where z ij = 1 and z ij = 0 indicate active and inactive statuses, respectively.Constraints (2h) prohibit mass flow across each regulator when z ij = 0. Constraints (2i) ensure that mass flow across each regulator is in the same direction as the loss in pressure.Constraints (2j) and (2k) model the remaining pressure dynamics.Here, each regulator has a corresponding scaling factor, α ij .This factor models the Natural Gas System Electric Power System Fig. 2. Diagrammatic representation of a small joint gas-power network.Here, D 1 contributes to the objective term η G (•) and L 1 , L 2 contribute to η P (•).
Finally, the linkage between gas and power systems occurs at relationship between junction pressures when the regulator is active, i.e., α ij p i = p j .The factor is further limited by the bounds Constraints (2j) and (2k) require that, when a regulator is active, pressures are defined according to the scaling relationship.Otherwise, the pressures at the junctions connected by the regulator are decoupled.Constraints (2l)-(2s) model compressors in the network.Each compressor (i, j) ∈ C models an increase in pressure at junction j ∈ J by a variable scalar α ij .Without loss of generality, bidirectional compression is not considered, although each compressor may allow for uncompressed flow in the opposite direction.These different behaviors of compressors are modeled by employing three different sets of constraints.The first are Constraints (2l) for compressors that prohibit reverse flow, where α ij and α ij are minimum and maximum pressure ratios.The second are Constraints (2m) and (2n) for compressors where reverse flow is allowed and α ij = 1.Note that here, if f ij < 0, then p i = p j .Finally, Constraints (2o)-(2s) model compressors where uncompressed reverse flow is allowed and α ij = 1.In this case, the behavior of each compressor is disjunctive in its flow direction.To model this disjunction, discrete variables y ij ∈ {0, 1} are introduced in Constraints (2o) to model the direction of flow through each compressor.Here, y ij = 1 indicates flow from i to j, and y ij = 0 indicates flow from j to i. Constraints (2p)-(2s) model the pressures and pressure differences between junctions as per the specified flow direction and compression ratio bounds.
The remaining Constraints (2t)-(2v) are variable bounds.Constraints (2t) are mass flow bounds, Constraints (2u) are pressure bounds, and Constraints (2v) are receipt and delivery bounds.Note that Constraints (2v) differ from the typical assumption of fixed supply and demand.These modifications ensure mass conservation is satisfied in damaged networks.

C. Interdependency Modeling
As in [23] and [15], gas and power systems are connected via heat rate curve models for gas-fired power generators, i.e., i:(i,j)∈K Each constraint links the real power generated at possibly multiple generators with a single gas delivery.Here, h 1,2,3 i are coefficients of the heat rate curve for i ∈ G, and K is the set of linkages between gas-fired generators in G and their corresponding gas delivery points in D G ⊂ D. Furthermore, h 1 i ≥ 0 for all (i, j) ∈ K, and thus the left-hand side is always a convex function.However, note that Constraint (3) is nonlinear nonconvex when h 1 i = 0. Finally, the presence of h 3 i z g i ensures that when z i = 0, the intercept of the heat rate curve, and thus both generation and gas required, will be zero when a generator is uncommitted from the dispatch scenario.
A diagramatic illustration of the joint network model is illustrated in Figure 2. Here, gas and power systems are linked by the single interdependency K 1 , which relates the delivery D 2 to the generator G 1 .Contributions to the gas and power delivery objectives, which are later described in Section III-A, are represented by green and red colored nodes, respectively.

D. Challenges
Although independent gas and power MLD models were explored by [4] and [5], respectively, the joint MLD problem that includes Constraints (1)-( 3) is more challenging.Most importantly, the nonlinear nonconvexities that appear in Models 1 and 2 arise primarily from different sources: Model 1 includes many nonlinear equations with bilinear variable products, whereas Model 2 includes more manageable quadratic nonlinear equations.To model them exactly, Model 1 must be formulated as a highly challenging MINLP, but Model 2 can be written as a more tractable mixed-integer nonconvex quadratic program.These differences suggest potentially incompatible numerical methods and solving technologies in practice.
Although omitted here for brevity, this work also employed a number of important preprocessing steps used to ensure the construction of feasible damaged joint networks that satisfy Constraints (1)-(3).Finally, Models 1 and 2 constrain each system using steady-state physical assumptions.In practice, modeling the transient dynamics of the gas system could be crucial.However, as will be shown in subsequent sections, even the steady-state variant considered in this paper is computationally challenging.This work is an important first step toward building MLD techniques that also consider transients.

III. MAXIMAL LOAD DELIVERY FORMULATIONS
This section derives the joint gas-power MLD formulations used throughout the remainder of this paper.First, Section III-A defines the competing objectives of the joint MLD problem.Section III-B poses lexicographic and weighted MLD formulations that prioritize the gas-power delivery tradeoff in different ways.Section III-C derives MICP relaxations of the MINLP MLD formulations.Finally, Section III-D summarizes the naming conventions used for these various MLD formulations, which are then empirically compared in Section IV.

A. Objectives of the Maximum Load Delivery Problem
The objective of the MLD problem is to maximize the amount of nongeneration gas and active power load delivered simultaneously under a multi-contingency scenario.Note that the maximization of nongeneration gas, specifically, allows the model to decouple practical objectives of the gas system (e.g., delivery of fuel for residential heating) from practical objectives of the power system.However, because the delivery of nongeneration gas load can inhibit the amount of active power generation, and thus active power delivered, there exists an important tradeoff between these two objectives.For notational ease, we first write the two objective functions as Here, Equation (4a) denotes the normalized sum of all prioritized nongeneration gas demand, where D ′ := D \ {j : (i, j) ∈ K} (i.e., the set of all nongeneration gas deliveries), and β i ∈ R + is a predefined restoration priority for delivery i ∈ D ′ .Similarly, Equation (4b) denotes the normalized sum of all prioritized active power loads.Note that for all of the experiments considered in this study, The tradeoff between nongeneration gas and active power load naturally lends the MLD problem to the broader category of multi-objective optimization.A thorough survey of multiobjective optimization methods in engineering is presented by [24] and describes a number of techniques for specifying preferences among multiple objective functions.These include weighted sum, weighted product, lexicographic, and bounded objective optimization methods.In Section III-B, we define lexicographic and weighted sum variants of the MLD problem.

B. Lexicographic and Weighted MLD Formulations
To explore the gas-power tradeoff, we introduce three MLD models that prioritize gas and power delivery in different ways.The first is a lexicographic formulation that maximizes the amount of nongeneration gas load delivered first.This situation is representative of common contractual requirements for gas grid operators.Here, the MLD is written as the program where η G (d * ) is the optimal objective when maximizing gas delivery alone.The second MLD is a similar formulation that maximizes the amount of active power load delivered first, i.e., maximize η G (d) Constraints (1)−(3). (MLD-P) The last is a single-level formulation that weights normalized sums of nongeneration gas and active power delivery, i.e., maximize λη where 0 < λ < 1 is a weighting parameter for the objective.Note that (MLD-G), (MLD-P), and (MLD-W) are mixedinteger nonlinear, nonconvex programs.The nonconvexities arise from three sources: (i) discrete operations of controllable components (e.g., z g i for generator commitment); (ii) bilinear products that appear in both gas and power network physics (e.g., V i V * j in Ohm's law); and (iii) nonlinear equations used for satisfying physical relationships (e.g., the Weymouth equation for pipes).In the following sections, we leverage a number of relaxations to render these problems more tractable.

C. Relaxation of Bilinear Products and Nonlinear Equations a) Convexification of Power Physics:
The primary sources of nonconvexity in Model 1 are the bilinear products that appear in Constraints (1a)-(1c) (e.g., V i V * j ).A large body of literature has developed relaxations of similar terms, and for a comprehensive review, we refer the reader to [6].In this paper, we develop a model based on a second-order cone (SOC) relaxation of the AC power flow equations, first presented by [25] and used for power MLD analysis in [5].
The primary insight of the SOC formulation is that variable products (|V i | 2 and V i V * j ) can be lifted into a higherdimensional variable space (W ii and W ij , respectively).This renders terms involving these products linear, and the relaxation in the new W -space is ultimately strengthened via This is an SOC constraint, lending the formulation its name.b) Convexification of Gas Physics: Many nonconvexities in Model 2 appear in the form of nonlinear equations (e.g., Constraints (2b)) and bilinear variable products (e.g., Constraints (2i)).To resolve both, direction variables y ij ∈ {0, 1} are first introduced for each node-connecting component (i, j) ∈ A. We also introduce variables π i ∈ R + to denote squared pressures p 2 i for i ∈ J .This first allows for a partial linearization of the Weymouth equations for pipelines, i.e., Then, variables ℓ ij for (i, j) ∈ P are introduced to model the difference in squared pressures across each pipe.The introduction of y, π, and ℓ, as well as convexly relaxing the equalities in Constraints (6), give rise to the convex relaxation Note that Constraint (7d) is the primary physical relaxation, i.e., the Weymouth equation need not be satisfied with equality.Convexification of the remaining nonlinear nonconvex terms in Model 2 is accomplished in a similar manner to the above.Here, for brevity, we omit the derivation of the full mixedinteger convex relaxation used throughout the remainder of this study.For a complete derivation and description of the relaxed mixed-integer convex model, we defer to [4].
c) Convexification of Gas-fired Generation: Constraints (3) are linear when h 1 i = 0 but nonconvex when h 1 i > 0. In the latter case, Constraints (3) can be convexly relaxed as where However, in our experiments, all h 1 i are zero, and the relaxation is not required.
2) (MLD-*-R): Formulations where power and gas constraints use SOC and MICP relaxations, respectively.These formulations provide different tradeoffs between model accuracy and computational performance.An empirical evaluation of both allows us to quantify the effects of the relaxations, as well as to guide our subsequent MLD analyses.

IV. COMPUTATIONAL EVALUATION
In the following, Section IV-A describes the networks, computational resources, and parameters used throughout the computational experiments; Section IV-B compares the efficacy of exact and relaxed MLD formulations on randomized N −k multi-contingency scenarios; Section IV-C evaluates the runtime performance of formulations over the same experimental sets; Section IV-D provides a proof-of-concept MLD analysis across the same experimental sets, illustrating the tradeoffs when lexicographically maximizing gas and power load delivery; and Section IV-E provides a proof-of-concept Pareto analysis of load delivery on a single joint network.

A. Benchmark Datasets and Experimental Setup
The computational experiments in this paper consider gas and power networks of various sizes that appear in the literature or have been derived by subject matter experts.These networks are summarized in Table I.The networks in this table are named according to the number of junctions in the natural gas network (e.g., NG11) and the number of buses in the electric power network (e.g., EP14).The references from which the gas, power, and/or joint network properties are derived appear in the second column of this table.The numbers of natural gas and electric power system components of the joint networks vary substantially and are specified in the second and third delineated portions of Table I, respectively.For networks that reference [27], heavily loaded variants of the corresponding electric power network datasets are used.
Joint gas-power network properties are summarized in the last column of Table I.Here, NG25-EP14 uses the linking and heat rate properties of the joint network instance developed by [29], and NG146-EP36 uses the properties of the instance developed by [23].Linkages within the NG247-EP240 network were derived from open data, and heat rate curves were estimated in a manner similar to [23].The remaining networks combine instances from GASLIB and PGLIB-OPF to create new joint networks of various sizes.The purpose of these new instances is twofold: (i) to explore the tractability of joint MLD instances as network sizes grow and (ii) to explore the tradeoffs involved in maximizing gas versus power delivery.In these new instances, the number of gas-fired generators, |K|, was estimated to be near min{0.25|D|,0.4|G|}, i.e., ≈ 25% of all gas deliveries or ≈ 40% of all generators.After determining the total number of gas-fired generators, the largest-capacity generators in each power network were then assumed to be linked to the smallest-withdrawal delivery points in the gas network.The heat rate at each gas-fired generator was then assumed to be equal to the proportion between the maximum withdrawal at the delivery point and the maximum power at the generator.Note that these networks thus use synthetically generated linkages between GASLIB and PGLIB-OPF instances, and these linkages are not necessarily reflective of real-world datasets.They are, however, instances where gas and power interdependencies are consequential, which in turn allows for a meaningful computational exploration of the MLD method.All of the MLD formulations considered in this paper were implemented in the JULIA programming language using the mathematical modeling layer JUMP, version 0.21 [31]; version 0.9 of GASMODELS, a package for steady-state and transient natural gas network optimization [32]; version 0.18 of POWERMODELS, a package for steady-state power network optimization [33]; and version 0.4 of GASPOWERMODELS, a package for joint steady-state gas-power network optimization [34].Furthermore, for the exact nonconvex nonlinear representation of Model 1 in (MLD-*), the polar form of the AC power flow equations, introduced by [35] and implemented by [33], was used.Similarly, for the exact representation of Model 2, the mixed-integer nonlinear nonconvex formulation described by [4] and implemented by [32] was leveraged.
Each optimization experiment was prescribed a wall-clock limit of one hour on a node containing two Intel Xeon E5-2695 v4 processors, each with 18 cores @2.10 GHz, and 125 GB of memory.For solutions of (MLD-W), version 0.7 of the open source JUNIPER MINLP solver was used [36].Within JUNIPER, IPOPT 3.12 was leveraged as the nonlinear programming solver, using a feasibility tolerance of 10 −6 and the underlying linear system solver MA57, as recommended by [37] for nonlinear network problems.Note that JUNIPER does not provide global optimality guarantees for (MLD-W), and feasible solutions obtained from the solver serve only as lower bounds on the true amount of maximum deliverable load.For solutions of (MLD-*-R), GUROBI 9.1 was used with its default parameterization.Here, since (MLD-*-R) is mixed-integer convex, globally optimal solutions are obtained.However, since (MLD-*-R) is a relaxation, a globally optimal solution corresponds only to an upper bound on (MLD-*).

B. Multi-contingency Damage Scenarios
This section examines the robustness and accuracy of the exact and relaxed weighted MLD formulations, (MLD-W) and (MLD-W-R), respectively, with λ = 0.5.Specifically, it studies these properties on large sets of randomized multicontingency or N −k scenarios, where k indicates the number of components simultaneously removed from the joint gaspower network.These scenarios are intended to capture the effects of severe multimodal network outages across joint systems.In each scenario, a random selection of 15% nodeconnecting components were assumed to be damaged (i.e., k ≈ 0.15N ).Through a parameter sensitivity study, we observed that this proportion of outages appeared to generate challenging MLD scenarios while providing interesting gas and power delivery tradeoffs among the coupled networks.For each network, one thousand such scenarios were generated.
Table II compares statistics of solver termination statuses across all N −k scenarios for each network and formulation.Here, "Conv."corresponds to the percentage of cases where the solver converged, "Lim." to cases where the solver time or other solver limit was reached, and "Inf." to cases that were classified as infeasible by the solver.Although both formulations are typically capable of converging on cases containing tens of nodes, for larger networks, (MLD-W-R) clearly outperforms (MLD-W), solving nearly all N −k instances.The results are especially dramatic for the three largest networks, where only one of three thousand (MLD-W) cases converges but all (MLD-W-R) cases converge.Note that three (MLD-W-R) cases are classified as infeasible due to numerical difficulties, but many more (MLD-W) cases are classified as infeasible due to the MINLP formulation and solver's greater tendency to converge to locally infeasible points.
Whereas Table II tion quality of relaxed formulations with feasible lower bounds obtained from (MLD-W).Here, "# Compared" corresponds to the number of cases used in each comparison, "Mean Obj." is the mean objective value obtained by (MLD-W) over all compared instances, "Mean" is the mean relative gap between (MLD-W) and (MLD-W-R) objective values, and "Median" is the median relative gap between objective values.In each such measurement, the relative gap is computed as where η is the objective value obtained when solving (MLD-W-R) and η is the objective value when solving (MLD-W).We note that, for NG25-EP30, five instances were excluded in the comparison: the three infeasible (MLD-W-R) instances and two instances that implied a negative relative gap.Proceeding with the analysis, the mean objective values for all sets of feasible solutions indicate that between around 50% and 75% of gas and power loads are being delivered across all multicontingency scenarios.Second, the mean relative gap between feasible solutions obtained by (MLD-W) and the upper bounds obtained by (MLD-W-R) are sometimes large, with the largest being 52.63% across all NG134-EP162 damage scenarios.
These extreme gaps have only two sources from which they can arise.First, a feasible solution obtained by JUNIPER for an (MLD-W) instance is not guaranteed to be near the globally optimal solution.That is, the globally optimal (MLD-W) objective value is potentially much larger than what JUNIPER reports at termination.Second, since (MLD-W-R) is a relaxation, it upper-bounds the globally optimal objective value of (MLD-W).The median column in Table III reports measures of centrality without the outliers that are likely arising from the first source of discrepancy.Through these measurements, (MLD-W-R) is observed to often provide reliable and tight bounds on the optimal objective of (MLD-W), with relative gaps most often ranging from nearly zero to less than 1.35%.This indicates that the relaxation is capable of providing tight upper bounds on maximum capacities of damaged networks.

C. Computational Performance
This section compares the performance of (MLD-W) and (MLD-W-R) using the instances described in Section IV-B.The performance profiles for these cases are depicted in Figure 3  tens of nodes; (M) networks containing hundreds of nodes; and (L) networks containing more than a thousand nodes (i.e., NG603-EP588).In all such categories, it is shown that the (MLD-W-R) formulation is able to solve substantially more problems than (MLD-W) in significantly shorter amounts of time.For joint networks with tens of nodes, both formulations are able to solve many instances within the one hour time limit.For networks with hundreds of nodes, (MLD-W-R) is capable of solving most instances within ten seconds, while (MLD-W) requires hundreds or thousands of seconds to solve only a small proportion.For networks with thousands of nodes, (MLD-W-R) solves all instances within ten seconds, whereas (MLD-W) does not solve any.The efficiency of (MLD-W-R) compared to (MLD-W) highlights its applicability to (i) realtime multi-contingency analysis and (ii) analyses that would require distributions of many multi-contingency scenarios.

D. Proof-of-concept Maximum Load Delivery Analysis
Whereas Sections IV-B and IV-C study the computational and accuracy tradeoffs between (MLD-W) and (MLD-W-R), this section provides a proof-of-concept MLD analysis using the (MLD-G-R) and (MLD-P-R) formulations on the same set of N −k damage scenarios.Figure 4 displays 18 histograms that evaluate the proportions of gas and power loads delivered across solved damage scenarios for the nine joint networks while using the two problem specifications.Here, green bars correspond to histogram frequencies obtained from analyzing results of (MLD-G-R) solutions (i.e., gas prioritization) and red bars correspond to (MLD-P-R) solutions (i.e., power prioritization).Brown, overlapping bars correspond to frequencies that appear in both (MLD-P-R) and (MLD-G-R) histograms.These results indicate qualitative differences in the hypothetical robustness of each joint network.They also display the extremal tradeoffs between prioritizing gas versus power delivery in the presence of extreme outages.Finally, they indicate the sensitivity of each gas or power network to the interdependencies that link them.These histograms serve as basic proofs of concept for real-world MLD analyses.The left half of Figure 4 displays histograms of maximum gas load delivered in the presence of severe N −k outages.First, note that these histograms display a variety of load distributions across the cases and networks considered.Some networks, e.g., NG25-EP30, NG247-EP240, and NG603-EP588 suggest gas grids that are highly sensitive to the outages considered, with large proportions of damaged networks often incapable of delivering more than 50% of gas load.Other networks, e.g., NG40-EP39, NG146-EP36, and NG135-EP179 show less severe but still substantial sensitivities to these outages.The remaining networks display gas network sensitivities somewhere between these two extremes.
The overlapping histograms also display the tradeoffs encountered when prioritizing gas versus power delivery.In the three joint networks NG25-EP14, NG134-EP162, and NG603-EP588, gas and power interdependencies are mostly inconsequential, and prioritizing either gas or power barely affects the maximum gas capacity.This is likely a result of excess generation capacity in the corresponding power networks.Other networks, e.g., NG40-EP39, NG146-EP36, and NG135-EP179 show more interesting tradeoffs, where prioritizing either gas or power results in substantial changes in the overall maximum load distributions.The remaining networks show less interesting tradeoffs, although NG11-EP14 displays large tradeoffs, likely due to the drastic effects that even minor outages can have on the relatively small network.
The right half of Figure 4 displays histograms of maximum active power delivered in the presence of the N −k outages.First, the four networks NG25-EP30, NG134-EP162, NG247-EP240, and NG603-EP588 appear robust to outages in the joint network and are often capable of delivering more than 75% of the original power load.The remaining networks see a greater variety in their maximum load distributions.Whereas some networks, e.g., NG25-EP30, NG134-EP162, and NG603-EP588, appear less reliant on gas-fired power generators, the remaining networks exhibit more drastic changes when prioritizing gas versus power delivery.The most extreme example appears to be NG146-EP36, which is often capable of delivering a large amount of power across all N −k cases when power is prioritized but also often loses more than 25% capacity when gas delivery is prioritized.
We remark that, to solve (MLD-G-R) and (MLD-P-R), inner-and outer-level problems of the lexicographic maximization are solved sequentially.For example, to solve (MLD-G-R), (i) the inner level problem maximizing η G (d) is solved, yielding a solution d * , then (ii) η P (z d ) is maximized, subject to Constraints (1)-( 3) and η G (d) ≥ η G (d * ) − ǫ.The latter ensures that nongeneration gas load delivered in the outerlevel is at least that of the inner level, minus some feasibility tolerance ǫ, taken in this study to be 10 −7 .A similar algorithm is used for (MLD-P-R).We note that the general algorithm is not as numerically reliable as (MLD-W-R) and does not solve 469 of the 18,000 N −k cases considered in this subsection.This could be alleviated with a larger ǫ or direct use of lexicographic features available in some solvers (e.g., GUROBI).

E. Proof-of-concept Pareto Analysis
Together, (MLD-G-R), (MLD-P-R), and (MLD-W-R) allow for a variety of prioritizations of gas versus power load.As such, they serve as powerful tools for exploring the wide range of possible MLD solutions based on the relative importance of gas versus power delivery.This can provide gas and power grid managers with best-case capacity estimates depending on the type of coordination between the two systems.In turn, this enables a better understanding of the extremely complex yet practically important tradeoffs encountered during the operation of a damaged joint network.Whereas Sections IV-B through IV-D focus on analyzing performance and qualitative aspects of MLD analyses across a large number of joint networks, this section focuses on providing a proof-of-concept Pareto analysis on a single joint network, NG146-EP36.
Figure 5 shows a linearly-interpolated approximation of the Pareto front for mean active power versus gas delivery across the same set of N −k scenarios considered in previous sections.Here, the upper-left and lower-right endpoints correspond to means obtained from the (MLD-P-R) and (MLD-G-R) problem formulations, respectively.Interior data points correspond to means obtained from the (MLD-W-R) formulation, where the tradeoff parameter λ was varied to determine interesting and distinct points on the Pareto front.
First, note that when prioritizing power delivery, on average, 96% of active power is delivered but less than 45% of nongeneration gas is delivered.When gas delivery is prioritized, 88% of power is delivered, while 50% of gas is delivered.Between these two extremes, the amount of gas and power increases and decreases, respectively, with increases in λ.For λ 0.5, active power decreases more slowly as a function of λ, and for λ 0.5, the rate of decrease appears larger.In this case, λ ≈ 0.5 happens to represent a value where (MLD-W-R) begins to prefer maximization of gas delivery over power delivery.Thus, in practice, a point near this value of λ could be one which maximizes simultaneous delivery of the two quantities while having practically equal prioritizations.

V. CONCLUSION
Recent increases in gas-fired power generation have amplified interdependencies between natural gas and power transmission systems.These interdependencies have engendered greater vulnerabilities to gas and power grids, where natural or man-made disruptions can require the curtailment of load in one or both systems.To address the challenge of estimating maximum joint network capacities under these disruptions, this study considered the task of determining feasible steadystate operating points for severely damaged joint networks while ensuring the maximal delivery of gas and power loads simultaneously.Mathematically, this task was represented as the mixed-integer, nonlinear nonconvex joint MLD problem.
Three variants of the MLD problem were formulated: one that prioritizes gas delivery, one that prioritizes power delivery, and one that assumes a linear tradeoff between the two objectives.To increase the tractability of these problems, a mixed-integer convex relaxation of the joint network's physical constraints was proposed.To demonstrate the relaxation's effectiveness, exact and relaxed MLD formulations were computationally compared across a variety of N −k scenarios.The relaxation was found to be a fast and reliable means for determining bounds on capacities of damaged networks.
Two proofs of concept were then provided to showcase the analytical power of the relaxed MLD problems.The first provided comparisons between prioritizing gas versus power delivery in an MLD analysis.These examples showcased the sometimes substantial tradeoffs that should be considered in extreme outage scenarios.The second proof of concept provided a Pareto front approximation of gas versus power delivery across N −k scenarios using a single joint network.These proofs of concept highlight that the efficacy of the relaxation-based MLD method makes it a potentially valuable tool for complex real-world decision support applications.
Future work will focus on extending the MLD approaches developed in this paper.First, additional gas and power relaxations should be considered to more accurately and efficiently scale to joint networks containing many thousands of nodes.Preprocessing routines, such as optimization-based bound tightening, may also aid in improving existing relaxations.Second, the current problem assumes the full coordination between gas and power systems when deciding operations that maximize load delivery.The modeling of bidding mechanisms that drive both systems could provide more accurate joint capacity estimates.Finally, capturing transient dynamics in gas networks is sometimes crucial for understanding the effects of network disruptions, which may only be realized long after the disruption occurs.Future work should consider these transient effects when modeling load delivery in the gas network.
Disruptive event begins.(ii) Event and cascading effects cease.Load restoration processes begin without network repairs.(iii) Network repairs commence.(iv) Load restoration is complete.(v) Network repairs are complete.

Fig. 3 .
Fig.3.Performance profiles comparing the efficiency of (MLD-W) and (MLD-W-R) over the N −k instances described in Section IV-B.Here, the performance profiles are partitioned into three categories for (S) networks containing tens of nodes; (M) networks containing hundreds of nodes; and (L) networks containing more than a thousand nodes (i.e., NG603-EP588).

Fig. 4 .
Fig.4.Histograms comparing the proportion of total gas and power load delivered across all solved N −k scenarios for (MLD-G-R) and (MLD-P-R) variants.The x-axis indicates the proportion of load delivered, and the y-axis indicates the proportion of solved damage cases that deliver load within an interval.

Fig. 5 .
Fig. 5. Pareto front approximation of total active power load versus nongeneration gas load delivered over one thousand NG146-EP36 N −k scenarios.

TABLE I SUMMARY
OF JOINT NETWORK DATASETS USED IN THIS STUDY.

TABLE II COMPARISON
OF SOLVER TERMINATION STATUSES OVER WEIGHTED OBJECTIVE MLD N −k MULTI-CONTINGENCY DAMAGE SCENARIOS.
measures the numerical reliability of exact and relaxed MLD formulations, TableIIIcompares the solu-

TABLE III COMPARISON
OF SOLUTION QUALITY FOR EXACT AND RELAXED JOINT MLD FORMULATIONS OVER ALL N −k CONTINGENCY SCENARIOS.