An integrated mixed integer linear programming model for resilient and sustainable natural gas supply chain

ABSTRACT Disruptions in the natural gas supply chain result in reduced throughput and associated emissions and losses, causing significant economic, environmental, and social impacts. Therefore, it is crucial to design supply chains that are resilient and sustainable to prevent or reduce the effects of disruptions. This paper proposes a novel Mixed-Integer Linear Programming model, which optimizes the natural gas supply chain in terms of resilience and sustainability, by examining the impact of an additional workflow design (contingency pipeline) located between the shutdown inlet and outlet nodes in the transmission echelon. The model is applied to a “real world” case, using data collected from gas companies operating in Nigeria. Both steady and transient states of the system are examined in this study through a set of scenarios. The best final solution was found to yield 93.6% performance increase when compared to target throughput and 63% performance increase with the introduction of the contingency when compared with the baseline scenario.


Introduction
Climate change concerns have led Governments around the globe to adopt legally binding commitments to cut down their greenhouse gas (GHG) emissions to limit the global warming effects (UNFCCC 2015). Today's society confronts the significant challenge of providing sustainable energy to meet the growing global demand (British petroleum 2019). In 2019, the energy consumption increased by 1.3% in relation to 2018, representing less than half of the 2.8% growth in 2018; the increase in consumption was met mainly by Natural Gas (NG) and renewables (British Petroleum 2020). NG is a significant player in the energy mix and a reliable energy fuel that has been reported to bridge the gap between conventional and renewable sources (Hamedi et al. 2009;Mac, Brouwer, and Scott Samuelsen 2018). Projections of the future global energy mix suggest that 85% of total energy supply growth will be generated by NG and renewables, with renewables becoming the largest power generation by 2040 (British petroleum 2019). The NG-fired power plant is a relatively low-carbon, flexible technology, which can be used to enable a low carbon energy transition (Ríos-Mercado and Borraz-Sánchez 2015; Sayed, Wang, and Bi 2019). Although NG is associated with strong GHG emissions at different supply chain levels, it is argued that it generates less emissions compared to other fossil fuels and promotes cleaner production (Balcombe, Hawkes, and D 2018;Emenike and Falcone 2020;Hao et al. 2016;Zhang et al. 2020). Furthermore, flexible NG-fired power plants can satisfy changes in peak demand and network congestions, resulting from intermittent renewable energy (such as wind and solar) in the network, challenging the power system security (Hutagalung et al. 2017;Ioannou et al. 2019).
The NG supply chain consists of several interconnected nodes (see Figure 1), composed of supply, transmission, distribution, compression, storage, and production sectors (interconnected by physical and financial infrastructure, information sharing, and conveyance) rendering the system broad and minimization of the transportation tariff and the maximization of the transported gas volume, with gas flow and compressor stations (Alves, de Souza, and Luiz Hemerly Costa 2016). Apart from cost minimization and profit maximization, researchers have investigated key objectives such as portfolio diversification, flexible contracts, capacity planning for transportation, safety stocks, and system redundancy (Emenike and Falcone 2020). These objective functions can promote the flexibility and resilience of the NG supply chain process. Emenike and Falcone (2020) stated that optimizations have been carried out on different levels of the supply chain echelon due to the complexity of the NG supply chain. Generally, NG optimizations are mostly evaluated on the transportation level involving the pipeline and compressor nodes (Abahussain and Christie 2013;Bopp et al. 1996;Farzaneh-Gord and Rahbari 2016;Sabri and Beamon 2000;Vasconcelos et al. 2013).
The transportation pipeline network is a complex, large-scale optimization problem with inherent nonlinearities, including, for example, the flow pressure loss constraints, and other hydraulic equations (Humpola and Fügenschuh, 2015). Several authors have developed models to simulate the nonlinear dynamics of the gas networks. Zhou et al. (2021) developed a Mixed-Integer Nonlinear Programming (MINLP) model for an underground pipeline network of gas storage (UNGS) that aimed at minimizing the total cost of pipelines, platforms, and stations. Considering the flow conditions of injection and production, the algorithm finds the best pipe network layout, the topology of the platform and the central station, and the pipe diameters of the ground pipe network of the system. Humpola and Fügenschuh (2015) used a MINLP model to study a pipe network design problem, which considered nonlinear and nonconvex potential-flow coupling constraints to define the relationship between the flow rate at an arc and the pressure at the end nodes. Kabirian and Hemmati (2007) developed a strategic planning model for natural gas networks such that the optimization of the nonlinear model addresses the short-run development plan where the location of compressor stations, pipeline routes, and sources of natural gas was considered to reduce transmission network cost while meeting increasing energy demand. Also, to improve the active control for gas transmission systems Sukharev and Kosova (2017) considered the problem associated with technical parameter identification in an unsteady state using a nonlinear model. Similarly, Mikolajková et al. (2017) developed a model of a pipeline network for NG distribution through a multi-period MINLP formulation, using the overall system cost, including fuel costs, investment and operation costs as the objective function. Constraints of the model included mass and energy balance for the network nodes, pressure drop equations, and gas compression in compressor nodes.
Nevertheless, such design optimization problems can become very computationally intricate. To address this, Hong et al. (2020) developed an integrated MILP method to derive the optimal gathering pipeline network, by considering the minimization of total construction cost as the objective function and using a piecewise method to linearize the nonlinear hydraulic equations. Zhang et al. (2017) also developed an integrated MILP method to derive the optimal production well gathering pipeline network toward minimizing the total investment while considering terrain, obstacles, and other constraints.
Disruptions in the NG supply chain can cause significant economic, environmental, and social impacts (Emenike and Falcone 2020) and a number of scholars have studies the design of resilient supply chains to prevent or reduce the effects of disruptions. Carvalho et al. (2014) studied the design of the NG supply chain aiming at reducing the impact of disruption due to factors outside the control of the business. The researchers proposed a decentralized method toward NG network resilience to failures when disruption occurs. Three strategies to manage congestions from disruptions were identified: network expansion, implementation of congestion pricing to cap the consumption of heavy users that cause network bottlenecks, and the grouping of consumers with similar suppliers' dependencies. One of the pinch points is the pipeline as it limits the quantity of flow based on its capacity available for NG transportation. The pipeline can also be affected by leakage or shutdown at the compressor station. To adopt a performance measurement, a new index that evaluates the functionality of a NG distribution network was proposed. It involved the restoration process after the occurrence of an external disruption, using recovery time as the main factor (Cimellaro, Villa, and Bruneau 2015). Focusing on the restoration problem of an interdependent infrastructure network following the occurrence of a disruptive event and the different interruption scenarios, Yasser, Barker, and Albert (2019) devised a resilience-driven restoration model to tackle and restore the system to normal state. Four key factors to be minimized were considered in the model, namely, the restoration time and the total cost as the sum of flow cost, restoration cost, and disruption cost (i.e., unmet demand cost). A decentralized algorithm model that controlled congestion in the NG supply chain affected by disruption, was proposed by (Carvalho et al. 2014). The model ensured that the available network capacity is distributed to each route without reducing network throughput. The approach considered a wide range of scenarios on a country-by-country basis within Europe and adopted a mitigation strategy. Specific indices were introduced to verify the results, based on per capita throughput and coefficient of throughput variation. The supply chain's resilience has attracted the attention of industry experts and researchers due to its considered competitive advantage (Balcombe, Hawkes, and D 2018;Christopher and Peck 2004;Moslehi and Reddy Pourhejazy et al. 2017).
This study presents a Mixed-Integer Linear Programming (MILP) framework for optimizing the resilience and environmental sustainability of a NG supply chain. To this end, the proposed model examines the impact of a relief pipeline contingency employed to address the unplanned disruption to a network and its subsequent closure. Compared to existing literature related to the NG supply chain optimization, the novelty of this paper lies on, firstly, the modeling of the impact of a relief pipeline redundancy on the NG supply chain as part of the mitigation planning problem (MPP) and secondly on the integration of resilience and environmental sustainability objectives in the proposed Mixed Integer Linear Programming (MILP) model.
The rest of the paper is structured as follows. In Section 2, the definition of the problem, model assumptions, and the overview of the scenarios investigated are presented. Section 3 outlines the mathematical programming model. In Section 4, the case study, composition of the infrastructure and the overview of the NG-proposed workflow are described. Section 5 describes the results generated in the steady and transient states. Finally, Section 6 presents the concluding remarks and discussion.

Problem definition
This study addresses the mitigation planning problem (MPP) of a natural gas supply chain by optimizing a relief pipeline contingency employed to address the unplanned disruption of the network. Disruption and subsequent shutdown of the compressor node results in emission losses, downtime, and shortage in throughput supply. Excess trapped NG in the mainline supplying the compressor station is usually emitted, inducing an environmental threat because of continuous flow from upstream. The NG network's initial performance describes the known operating state; the optimization model is then applied to study the impact of the disruption.
The main features and assumptions of the proposed MILP model are summarized as follows (the list of symbols can be found in the Nomenclature): • The given planning horizon is divided into equal time intervals tPT over a period of 30 months. • The inlet nodes can include the supplier ðiÞ, the processing plant jÞ ð , the compressor ðkÞ, and the city gate ðgÞ nodes. • Outlet nodes can include the processing plant ðjÞ, the compressor ðkÞ, the city gate ðgÞ and the consumer ðqÞ nodes. • Shutdown and startup periods t 2 T are defined for every node k. During shutdown period, t, there is a loss of gas represented as Z E t ð Þ . The startup is defined as the time the plant requires to start running, featured as Y k;t ð Þ Shutdown is denoted as R k;t ð Þ , while the operating time is represented as X k;t ð Þ . • The demand volume for consumer, d a m;t ð Þ , is associated with a dedicated capacity, rc max m;t ð Þ . • During the shutdown, R k;t ð Þ the loss Z E t ð Þ through emissions is recorded for a time duration.
• The impact of disruption on the gas flow to the plant is bounded by the minimum S min k ð Þ and maximum S max k ð Þ mass flow rates. • The nodes from the supplier to the consumers are interconnected. There are no dedicated storage units. Pipelines provide limited temporary storage for every given time. • It is assumed that no more than two plants are simultaneously shutdown. • Disruption to the network nodes is identified as the main cause for the shortfall in supply. In steady state, the pipelines are operating at constant flow rates with a constant pressure profile, which may differ from point to point, but does not change with time (Menon 2005). In this work, The General Flow Equation (GFE) for the steady state isothermal flow in a gas pipeline is introduced, which is the basic equation of the pressure drop as a function of the flow rate. If the inlet pressure at the upstream is constant in the steady state, the flow rate will increase if the downstream outlet pressure is reduced. The pressure drop was modeled following Menon's (2005) equation as follows: where f=friction factor P b =base pressure P 1 =upstream pressure P 2 =downstream pressure G=gas gravity, (air = 1.00) T f =average gas flowing temperature (°R (460+°F)) L=pipe segment length Z=gas compressibility factor at the flowing temperature (dimensionless) D=pipe inside diameter • The initial state of the network is static except for the introduced relief pipe, while the dynamic state of the system is also studied. • The different statuses of the plant nodes include the following: operational, shutdown or start up.

Overview of the NG network model and alternative pathway (relief pipeline) component
An overview of the NG supply chain infrastructure with the additional workflow (pipeline redundancy) is shown in Figure 2. A detailed description of the system can be found in Emenike and Falcone (2020). The relief pipeline is located in the midstream between the gas processing plant and the distribution center (red colored line). The city gate station (CGS) regulates the natural gas (NG) pressure by using expansion valves (Sheikhnejad, Simões, and Martins 2020).
When the emergency shutdown occurs, the gas flows through emergency outlet between the valve and the compressor station. Without this relief pipeline, the gas is emitted through a relief valve to the environment, resulting in throughput losses and pollution to the environment. This relief pipe gradually flows the gas to a sale line or to another compressor station depending on its proximity to the sale line.
The introduction of the relief pipe is to ensure continuous flow and help the network withstand the impact of the disruption. The relief pipeline takes the excess flow during shutdown so that the pressure at the end does not increase excessively, leading to a reverse flow (where the pressure at the end of the pipe is greater than that at the inlet). The mathematical model describing the necessary conditions for the optimization process is presented in Section 3.
If the flow is constant in time and the pipeline is straight and horizontal, then according to Tomasgard et al. (2007), it can be assumed that the system is in steady state, and the time resolution is strategic or tactical and not operational. Table 1 summarizes all optimization scenarios implemented in Section 5. The baseline scenario (BS) is used to benchmark the compressor's performance and throughput without introducing the relief pipeline (emergency outlet). The mean throughput of the BS constitutes a key performance indicator used to compare the scenarios. The redundancy (relief pipeline) is introduced in the other scenarios, and it is set to operate only when the disruption occurs. In scenario 1, the redundancy is introduced through the opening of the alternative pathway valve, allowing the flow into the alternative pathway in the steady state. Scenario 2 introduces flow constraints in an extended time-series to investigate the impact of flow constraint when the redundancy is operating. Scenario 3 comprises the extended time-Figure 2. Case study system layout with additional workflow (adopted from (Emenike and Falcone 2020)). series without the flow constraint. The extended time-series in scenarios 2 and 3 is required to forecast future data points in the planning horizon and to eliminate possible deviance in the data. Two scenarios are analyzed under the transient state. Scenario 4 analyses pressure variation resulting in pressure surge from the plant closure, while scenario 5 analyses the prolonged closure in both the inlet and outlet nodes on the mainline. The scenarios in both steady and transient states are summarized in Table 1.

Mathematical programming model
A mixed integer linear mathematical programming model is adopted in this work. GAMS mathematical programming system for optimization supports interfaces with several optimization algorithms or solvers. The GAMS programming model is considered a reliable optimization tool for mathematical modeling of the supply chain, where the run time varies based on the objective to be achieved.

Model Formulation
In this section, the MILP optimization model describing the gas supply chain with the disruption and loss elements is introduced, along with the alternative pathway that can serve as a capacity for expansion. The optimized model formulation aims to make the supply chain more resilient and sustainable. The detailed list of symbols used in the paper can be found in the Nomenclature.

Objective function
The objective of the optimization problem is to maximize the resilience of the supply of NG using flow volume flexibility from supplier to consumer nodes and minimize the associated emissions during plant shutdown. In the context of this study, the resilience has been approximated through the minimization of the throughput losses of the supply chain system. The flexibility of the supply chain nodes will help achieve the targeted resilience. For simplification, the multi-objective optimization problem has been formulated as a single-objective function, which is composed of the aggregate volume flexibility as a function of flow (represented as Z1), as well as the losses and emission savings, the operating status of the plant, and the additional flowline (Z2). The same weight was applied to all individual items of the objective function. As such, the objective function can be expressed as: zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl }|ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl { Where, • Supply node volume flexibility (SVF), which corresponds to the maximum inlet flow from the supplier to the processing plant multiplied by the mass flow rate (from node i to j in time t). This ensures that the amount of product from the supply node is within the capability of the processing plant node. • Processing node volume flexibility (PVF), which corresponds to the maximum inlet flow from the processing plant to the compressor multiplied by the mass flow rate (from node jto kin time t). This ensures that the amount of product transported from the processing plant node is within the capacity of the compressor node. • Transmission node volume flexibility (TVF), which corresponds to the maximum inlet flow from the compressor plant to the consumer node multiplied by the mass flow rate. This is a function of node k, m and g multiplied by maximum pressure in the compressor. This ensures that the amount of product from the compressor node is within the capacity of the consumer node. • Operating status (OS), which corresponds to the operating status of the plant with respect to the minimum run time after startup and the minimum shutdown time. • Additional flowline (AF), which corresponds to the flow to relief pipe from the main line during shutdown multiplied by the maximum flow rate. This describes the redundancy introduced to mitigate the impact of the shutdown.

Supplier capacity constraint
The supplier node represents the gas fields, which is often owned by multiple parties with production rights to produce in commercial quantities. Constraint (2) ensures that the supplied NG to the production plant is less than or equal to the supply capacity and the production plant capacity (Zhou, 2021). The total NG volume from all related gas wells does not exceed the maximum production capacity of gas fields in the supply nodei, as shown in constraint (3):

Compressor capacity constraint
The compressor is one of the vital nodes in the gas network system. It is used throughout the natural gas network to move gas from the upstream to the midstream and finally to the downstream at different pressures. It helps to exert pressure that has been lost due to friction in gas pipelines (Menon 2005) and to reduce volume by providing the necessary force to move the gas along the pipeline.
Constraints (4) and (5) ensure that NG flow from the processing plant node, j; to the compressor node,k, minus the gas losses and emissions at time, t, do not exceed the compressor capacity. The shutdown of plant node k is taken into consideration when there is a flow from plant jto plant k to account for the losses during plant disruption. Constraint (6) ensures that if node k supplies to the consumers, capacity must not exceed that of plant node k. The constraint is represented below:

Power plant capacity constraint
In constraint (7), the power plant is being supplied directly from plant node k. In this constraint, it is expected that NG flow from the plant node k to the transmission pipeline does not exceed the power plant capacity. It is assumed that the shutdown of plant node k affects primarily the supply of NG to the power plant.

City gate station capacity constraint
The city gate is a point at which a local gas utility receives gas from a transmission system. It supplies gas to the customers in the city at required consumption pressure. When the city gate station is opened, all NG flow from compressor node, k; in the transmission pipeline should not exceed the city gate capacity. This is ensured in constraint (8).

Consumer demand constraint
Based on the contractual agreement at every time period, consumers' demand should be satisfied by the suppliers as shown in constraints (9) Storage capacity constraint Constraints (12) -(14) suggest that the upper and lower limits of the temporary storage capacity should not be surpassed. In constraints (12) and (13), the gas sent to the pipeline for storage should be less than or equal to the capacity of the line packing storage capacity. Constraint (14) represents the minimum and maximum inventory storage levels, i.e., it indicates that the gas storage must fall between its minimum and maximum limits. The parameter y S w;t ð Þ represents the initial inventory in the storage, while the variables YM þ k;w;t ð Þ and YM À w;k;t ð Þ are the inflows and outflows from and to the compressor and the storage.

Material input-output balance constraint
In constraint (15), as there is no mass buildup in any node of the network, each node will be constrained according to the mass balance law (Zhou et al., 2021). Thus, for every node of the network system: P flow in ¼ P flow out. NG supplied from the gas fields to the processing plant equals NG transmitted from the processing plant to the compressor: Furthermore, the total NG supplied from the processing plant to the compressor and from the compressor should equal the NG from the compressor to the power plant and from the compressor to the city gate station, as shown in constraint (16).
The local distribution company (LDC) oversees the operation of the city gate station. Constraint (17) shows that the NG supplied from the compressor to the gas company equals the NG supplied from the gas company to the industrial and domestic consumers: Pressure equality constraint A simple maximum flow restriction is introduced. It is assumed that the distance between the pipeline nodes is limited in length. In the steady state, the pressure in and out of a node remains constant with time (the pipeline is assumed to be straight). The flow is assumed to be isothermal, namelyΔT ¼ 0. Constraint (18) refers to the transmission pipeline, from pipeline to compressor and from compressor to pipeline as represented below: Pressure inequality constraint Constraint (19) ensures that the nodal pressure does not exceed the maximum level at time t (Zhou, 2021). Also, pressure in the outlet node from processing plant j to compressor should not exceed the maximum pressure, so that in Constraint (20), pressure is kept within acceptable levels: P in j;k;t ð Þ � P out j;k;t ð Þ ; "j 2 J; k 2 K; t 2 T The pressure movement at both the mainline and the alternative pathway in time-series during the shutdown period in the transient state is represented in constraints (21) and (22). The inlet and outlet pressures are multiplied by the mass flow rates to derive the pressure variation.
The introduced redundancy measure operates within maximum and minimum pressure bounds for each period according to constraint (23). Constraint (24) displays the time when there is a flow from nodes j to k and from k to z during shutdown such that if no flow is recorded from either node at a defined time, it does not affect the pressure balance. The bigM z;t ð Þ represents a large number.

Capacity expansion constraint
A lower and upper bound for the cumulated capacity for expansion is introduced in constraint (28). By introducing the compression factor, the redundancy is modified in constraint (29). Then, the proportional capacity for expansion in constraints (30) and (31) is to ensure that the proportional capacity is not more than the capacity before expansion multiplied by the maximum proportional capacity for expansion and is not less than the capacity before expansion multiplied by the minimum proportional capacity for expansion, proposed by (Liu, 2013). An additional constraint is introduced when the redundancy is fully operational, and the trapped gas is rerouted as demonstrated in constraint (32).

Flow constraint
Constraint (33) ensures that the total losses do not exceed the gas flow supplied to the relief line during the shutdown. A corresponding upper and lower bound of the flow before and during the shutdown is introduced in constraint (34). Constraint (35) ensures emission losses during the shutdown of the compressor plant do not exceed the capacity of the alternative pathway, this constraint should be ignored if flow to alternative pathway does not occur and where binary for the alternative pathway is R (k,t) = 0.

Startup and shutdown constraint
The emergency shutdown is temporally bound by the assumption that all state variables are constant and that there is no mass accumulation in the system. Constraints (36)-(38) of the model indicate the start-up and the shutdown of the compressor station through:The following binaries are introduced relating to startup and shutdown actions of a plant node.
• X k;t ð Þ = 0 if there is no start up for plant node k at the beginning of the period; otherwise, 1. These constraints show that if plant node k is already running (i.e., X k; t ð Þ ¼ 1) at the beginning of the planning horizon then shutdown R k; t ð Þ ¼ 1. If plant node k starts operating at the start of the planning horizon Y k; t ð Þ ¼ 1. In constraint (36), the emission that occurs during shutdown is taken into consideration. However, simultaneous startup Y k; t ð Þ and shutdown R k; t ð Þ is disallowed.
Constraints (39) and (40) introduce the minimum online time for plant node k after its startup. It is expected that the plant will operate for a given period after its startup. The total period that plant node k has been operating continuously since its last startup is greater than the minimum online time.
Constraints (41) and (42) ensure that the total time that plant node k has been shutdown continuously is less than the minimum offline time.
In constraints (43) and (44), the maximum idle time is the maximum time that plant k is switched off continuously after its last shutdown, which is expected to be longer than plant shutdown time (R k;t ð Þ = 1).

Case study input data
A typical network has been selected as a case study (key components are shown in Figure 2). The case study comprises three gas streams that converge into a single supplier node (iÞ. The system includes a single processing plant node j1 ð Þ, four compressors in each compressor station (k1; undefined, undefined), two main transmission pipelinesðp1; p2) and one natural gas company represented by a city gate station (q1Þ. Although the model is able to consider different types of consumers (such as domestic, power production and industrial consumers), only the demand of the power plant consumers (m) is assumed in this study. The contractual volume obligation for power plant consumers is 360 mmscfd. It is assumed that this demand should be met for the entire planning horizon. The case study is applied to a "real world" case, using data collected from gas companies operating in Nigeria.
Each section of the supply chain is represented with a node, with the performance of prior nodes affecting the activity in subsequent nodes. The contractual volume obligation for power plant consumers is 360 million standard cubic feet per day (mmscfd). It is expected that the demand should be met monthly for the entire planning horizon. The values of the case study are analyzed over a 30month planning horizon. For the period under review, it is assumed that the disruption occurs in three different periods at times t8; t19; and t27 respectively, over the planning horizon. The problem of the study does not include dedicated storage. Reference parameters used within the time horizon are summarized in Table 2.

Results
The examined optimization problem was modeled using the General Algebraic Modeling System (GAMS) 26.14 with the CPLEX solver 12 in an intel ® core ™ i7, under standard configurations. The optimal solution was found within reasonable solution time with a zero-optimality gap. Data obtained is used to compute the baseline flow rate where monthly shortages were identified. In the BS, the mean flow rate at the beginning of the planning horizon provides a reference point to assess and compare the other scenarios.

Steady state
In a steady state, the values representing the gas flow of the system are independent of time. The computation is made in a deterministic environment where all parameters, constraints, and objective function are known. As such, the performance of the supply chain under the steady state (with no pressure variation with time) including the redundancy component can be determined.

Scenario 0: Baseline scenario
Under the baseline scenario, the mean flow rate is displayed in Figure 3. The mean flow rate amounts to 200.38 mmscfd, which means that the target mean flow of 360 mmscfd is not achieved and there is a shortfall of 159.62 mmscfd. To optimize the supply chain, the topology redundancy is introduced as a mitigation strategy in Scenario 1.
The performance level of each compressor, k, with respect to the corresponding minimum mass flow rate when X k;t ð Þ ¼ 1 and Y k;t ð Þ ¼ 1 is displayed in Figure 4. This is calculated by multiplying the minimum mass flow rate with the operating time of the compressor. Improved performance of the compressors is seen toward the end of the planning horizon with k2 and k3; outperforming k1, k4.

Scenario 1: Shutdown with redundancy
This scenario describes the plant shutdown with redundancy, and it operates when the disruption binary (R k;t ð Þ ¼ 1Þ is activated. The main pipeline node closes and the trapped gas between the inlet and outlet nodes flows through the alternative pathway. Figures 5 and label Figure 6. label display the throughputs (flow rates) following the occurrence of a disruption at the specific times in the planning horizon with the presence of the relief pipeline.
illustrates the gas flow rates for a standard mass flow rate in node k, while in Figure 6, the throughput reflects a varying mass flow rate introduced in node k to allow the adjustment of the flow to the disruption. For the latter case, the mean flow rate is increased from 200.38 (BS) to 327.67 mmscfd.
The outputs illustrated in the graphs correspond to the flow rates with no pressure drop with time between the inlet and outlet nodes of the pipeline. Each node in the network is bound within lower and  upper pressure limits. The mass balance constraint is applied to all relevant parameters, subtracting the emission losses resulting from the compressor plant shutdown. The improved flow rate is obtained by relaxing the disruption period such that the shutdown time is defined. The shutdown of the compressor station means that at least one compressor plant among k1 to k4 in the mainline is not operating during the entire shutdown period. To achieve the results in Figures 5 and label Figure 6. label, node k2 must be operating at the start of the planning horizon, such that the scheduling of supply to consumers comes from one to two compressors at any given time.
The performance level (throughput) of each node k is illustrated in Figure 7, where there is gas flow from the mainline to the relief pipe during the shutdown at time t R F k;z;t ð Þ ¼ 1 À � , and when X k;t ð Þ ¼ 0 and Y k;t ð Þ ¼ 1. To assess the performance of the node (k1 to k4), the minimum mass flow rates is  multiplied by the operating time of the compressor. While k3 and k4 remained unchanged, the performance of k1 increased and k2 varied in relation to the baseline scenario.

Scenario 2: Flow Constraint in Extended Time
In scenario 2, the upper and lower bound flow constraints (34) and (35), are introduced. The flow constraint introduces the upper and lower bound of the flow before and during the shutdown to ensure the total flow entering the inlet node and the total flow leaving the outlet node are within bound limits. This scenario is a combination of extended time sequence at successively equally spaced points and flow constraints when the redundancy is operating. To capture all possible changes that may affect the throughput, the input series is extended twice the original length to the end of the corresponding period by halving the number of each time steps throughout the planning horizon. The extended time-series is introduced where the impact of the flow constraint on the throughput is investigated such that the flow from the processing plant to the compressor is subject to the min/ max mass flow rate of the operating status of the plant. Additional operating time of the plant does not affect all the compressors at the same time. As displayed in Figure 8(a), the only exceptions are the shutdown times of t8, t19 and t27. For the normal period (Figure 8(b)), the no operating time is at t9 where only k3 is affected, while for the extended time (Figure 8(c)), the no operating time is found at t8.5 where k3 is not operating, and t19:5 where k1; k2 and k3 are not operating. The average throughput resulting under Scenario 2 conditions amounts to 336.078 mmscfd.

Scenario 3: No Flow Constraint in Extended Time
Scenario 3 is introduced to evaluate the impact on the overall system performance when the flow constraint is removed from the optimization problem. The extended time-series is also introduced without a corresponding upper and lower bound flow constraint. The flow rate performance for normal, extended and cumulative time in this scenario is shown in Figure 9(a-c). The defined shutdown times (undefined, undefined, and undefined) performed at optimal which is seen in the average throughput. In fact, the average throughput improved from 336.078 mmscfd (with the flow constraint) to 336.900 mmscfd (without the flow constraint).

Transient state
The transient effects of time-varying demand for natural gas affect the compressor and pipeline operations mainly just ahead of the delivery point. Studying the transient condition is restricted to the mainline transmission node with an extended observation time until the opening of the relief valve. The mainline closure during the disruption produces a pressure buildup. The gas accumulation suggests an expected pressure rise with time after the closure at these nodes. The gas compressibility allows for continuous pumping of gas from the upstream over a period, which eventually increases the line pack in the midstream and downstream echelons.

Scenario 4: Pressure surge from plant closure
In scenario 3, the transient state is investigated during a cumulated period of 24hrs. Each hour is divided into 6-min intervals, making a total of 10 points for every hour. The operating status of the plant node between the gas plant and the compressor is multiplied by the binary variable X k;t ð Þ À � for operation on lower and upper bounds of the flow to determine the variation. The pressure interaction in the outlet node is determined by multiplying the binary variable R k;t ð Þ À � with the maximum pressure upper and lower bound limit of the disruption in the mainline. Figure 10 shows the mainline pressure variation time series at maximum mass flow rate. In Figure 11, the mainline pressure variation with time is determined using the minimum mass flow rate. During the shutdown, the pressure time-series in the mainline in both cases rises to maximum pressure between time points 0.48 to 1.67 at approximately 13:50 hrs before dropping to a stable pressure rate of 0.765 and 1.15 psia, respectively, from point 168 to 222 at approximately 10:00 hrs.
Also, the pressure interaction is examined when the alternative pathway is operating. Ignoring the bound limit of the inlet and outlet pressure while introducing the mass flow rate for the alternative pathway, the time-series at the point of variation is separated between the normal and extended time to examine the pressure interaction. As shown in Figure 12, there is close interaction between the pressure variation in both the normal and extended times displayed when the compressibility factor equals 1. This means that the deviation of the real gas from the ideal gas is insignificant and therefore does not affect the throughput in the alternative pathway. If the compressibility factor is less than 1 (z < 1), the pressure interaction in the redundancy is then shown in Figure 13. A reduced mass flow rate from 400 mmscfd in the mainline to 120 mmscfd in the relief pipe during the plant shutdown without changing the originating pressure will cause a pressure rise, which will be compensated as the gas enters the sale line. At this point, the relief pipe and the sale line are operating simultaneously.

Scenario 5: Prolonged inlet and outlet nodes closure
Scenario 4 explains the unsteady condition caused by prolonged closure of both the inlet and outlet nodes on the mainline. If the closure of the nodes is within an allowable time based on the pipeline interim storage capacity, the gas will continue to flow from the upstream, which is temporarily stored until the inlet and outlet nodes are opened. On the contrary, if the disruption period on node k exceeds the projected allowable time, then the accumulated gas is emitted because of continuous pressure rise.  The pressure behavior is examined over a 24hrs period at 6 mins intervals in the mainline and the relief pipeline represented as inlet and outlet nodes. The pressure in the mainline also known as the inlet pressure as shown in Table 3 starts to increase at approximately 11:54hrs triggering the outlet pressure to decrease once the alternative pathway is opened. At the point when the control valve in the alternative pathway is opened, the outlet pressure is relatively stable but changes slightly over time. The variation becomes more evident over time as gas continues to enter the alternative pathway. The pressure then normalizes as gas is fed into the sale line from the alternative pathway. Assuming the mainline inlet valve is re-opened, pressure begins to increase at approximately 17:48 hrs, which is offset as the gas begins to flow into the relief line.
The flow rate when the relief valve is opened (operating alternative pathway) and when the compression factor equals 1 with extended time is shown in Figure 14. This indicates that the effect of the shutdown is minimized. The optimized flow rate in Figure 15 is obtained when the pressure bound limit is introduced. The impact of the pressure change on the flow rate is given as 321.17 mmscfd, which corresponds to an improved flow rate. In Figure 15, the pressure bound limit is also introduced but the compression factor is larger than 1 resulting in a flow rate output of 327.03 mmscfd.
In Table 4, the results from all scenarios in both steady and transient states are presented. The best optimal result in comparison to alternatives is compared across all feasible alternatives.
The throughput performance of the system improved by over 63% when comparing scenario 5 to the baseline throughput performance of 200.38 mmscfd. Moreover, 93.6% of the target throughput performance of 360 mmscfd has been achieved when considering scenario 3. The best optimal solution was found in the steady state scenarios, when the flow constraint is removed. Figure 13. Outlet pressure when compressibility factor is lower than 1.    Optimized flow obtained if the capacity of node k is the same for all period in the planning horizon.

327.67
Optimized flow is obtained if the capacity of node k varies at different rates in the planning horizon.

276.38
Optimized flow is obtained when the flow constraint is introduced. 336.078 Optimized flow is obtained when the flow constraint is removed. 336.90 Transient Pressure variation in extended time.
200.38 When compression factor equals 1 with an extended time.
321.17 When pressure bound limit is introduced in extended time.
323.37 The compression factor is increased in extended time. 327.03

Concluding remarks and discussion
A novel MILP model has been presented to study the impact of a proposed workflow on a natural gas supply chain when interrupted. The proposed optimization framework investigates both steady and transient state scenarios, taking into consideration relevant constraints. The system's interruption is introduced at specified time instances and the relief pipe is activated when the inlet and outlet mainline nodes are closed. The relief pipe serves as an alternative pathway to mitigate the disturbance. The results have shown an improvement to the resilience of case study NG supply chain, following the introduction of the contingency. More in specific, the performance of the system was improved achieving an additional throughput of 127.29 mmscfd with no pressure drops when the contingency pipeline was introduced. An even better performance was obtained with the introduction of the flow constraint (Scenario 2) with a throughput of 336.078 mmscfd and a slight improvement of 336.900 mmscfd when Figure 15 time-series data is extended, and the flow constraint is removed (Scenario 3).
In the transient state, findings show that as the downstream pressure is reduced, the flow rate will increase to keep the upstream pressure constant. The additional pathway can remain open even after the mainline valves are re-opened, providing a two-way simultaneous flow to compensate for shortages pending when supply is improved. The increased compression factor in the extended time produced a flow rate of 327.03 mmscfd, which is an optimal solution in the transient state.
The developed model can be adopted in other supply chain systems after appropriate modifications are introduced based on the peculiarity of the system under review. Although the case study uses a relatively simple pipeline network, the model is versatile and able to simulate more complex problems with a larger number of nodes and pipes (e.g., system process integration) by exhibiting global optimum with well-defined solutions. The deployment of the alternative route for gas flow during the plant's shutdown results in economic cost implication, which has not been considered in this work. To the best of the authors' knowledge, this is the first work that addresses the optimization of the NG supply chain focusing on resilience aspects and in specific the optimization of throughput and emission losses, while considering an additional redundancy pipeline design aiming at mitigating the disruption effects. This research can be further developed to introduce more scenarios and perform a cost-benefit and lifecycle analysis to assess the cost implication compared to environmental and economic benefits.
Future studies in this area could further investigate the following aspects: • Identification of the most suitable location to introduce the redundancy depending on the network's need over the planning horizon to satisfy demand and loss reduction is a critical issue, as a wrong location may induce significant challenges. • Further research optimization modeling should consider savings on downtime and how downtime minimization will lead to profit maximization for the system operators. • A similar study could be carried out in a stochastic environment where logical consideration of uncertainty can help estimate future expectations, calculate likely returns, and estimate associated risks.

Indices/Sets
I Set of all suppliers, i 2 I J Set of processing plant producers, j 2 J K Set of compressors, k 2 K G Set of city gate stations, g 2 G M Set of power plant consumers, m 2 M W Set of gas storage stations, w 2 W Q Set of industrial consumers, q 2 Q T Set of time periods, t 2 T P Set of pipelines, p 2 P Z Set of relief pipelines, z 2 Z