Development of an integrated dynamic model for supply security and resilience analysis of natural gas pipeline network systems

Hui Su ，Enrico Zio ，Zong-Ji Zhng ，Chng-Zhng Xiong ，Qin-Shng Di ，Qing-Wi Wu ，Jin-Jun Zhng ，*

a National Engineering Laboratory for Pipeline Safety/ MOE Key Laboratory of Petroleum Engineering/Beijing Key Laboratory of Urban Oil and Gas Distribution Technology，China University of Petroleum-Beijing，102249，Beijing，China

b Dipartimento di Energia，Politecnico di Milano，Via La Masa 34，20156，Milano，Italy

c MINES ParisTech，PSL Research University，CRC，Sophia，Antipolis，F(xiàn)rance

d PipeChina West East Gas Pipeline Company，Dongfushan Road 458，Pudong District，200122，Shanghai，China

e China Special Equipment Inspection and Research Institute，Heping Street 2，Chaoyang District，100013，Beijing，China

Keywords:Natural gas pipeline networks Dynamic modeling State space model Graph theory Resilience of natural gas supply

ABSTRACT An integrated dynamic model of natural gas pipeline networks is developed in this paper.Components for gas supply，e.g.，pipelines，junctions，compressor stations，LNG terminals，regulation stations and gas storage facilities are included in the model.These components are firstly modeled with respect to their properties and functions and，then，integrated at the system level by Graph Theory.The model can be used for simulating the system response in different scenarios of operation，and evaluate the consequences from the perspectives of supply security and resilience.A case study is considered to evaluate the accuracy of the model by benchmarking its results against those from literature and the software Pipeline Studio.Finally，the model is applied on a relatively complex natural gas pipeline network and the results are analyzed in detail from the supply security and resilience points of view.The main contributions of the paper are:firstly，a novel model of a complex gas pipeline network is proposed as a dynamic state-space model at system level;a method，based on the dynamic model，is proposed to analyze the security and resilience of supply from a system perspective.

1.Introduction

Natural gas plays a crucial role in the World energy portfolio.The International Energy Agency estimates that 25%of the World's energy will come from natural gas by 2030，while the global natural gas consumption will double the 2008 level (International Energy Agency，2008).Then it is not surprising that the governments around the World are becoming increasingly serious regarding the natural gas supply security and the pipeline network infrastructure which assures it by connecting the demands and the sources through large space distances (Zhu et al.，2017).

Natural gas pipeline networks are complex systems composed by a large number of units and sub-systems，varying from structures and functions.Because of the physical property of natural gas，the network system presents complex dynamic behaviors in time and space.These characteristics increase the difficulties to assess the supply security of the gas grid system.In general，natural gas flows from suppliers to customers under the driving force of the gas pressure.Because of the friction between the pipe inner surface and the gas，there are pressure drops during the process of gas flowing.The drops are compensated by compressor stations located at specific intervals in the pipeline network system.The transmission system operators (TSOs) control the pressure to balance the demands，the system transmission capacity and the capacity of gas suppliers.Besides，storage facilities，especially underground gas storages (UGSs)，are installed to increase the flexibility of the system in response to supply disruptions，demand peaks and congestions of the pipeline networks.

Many efforts have been performed in the development of models of natural gas pipelines and pipeline networks.Most of the works can be classified into two groups:steady models and dynamic models.In steady models，pressures and flows are hypothesized to be constant in time;on the contrast，the evolution of physic parameters in the pipelines is accounted for in dynamic models.

Steady models are mostly developed based on the balance of the input and output of the gas and are often applied to solve complex problems，e.g.，optimization problems and parameter analysis problems.For example，Szoplik (2016) studied the relationship between system performance and the temperature of air based on a steady state simulation;üster &Dilaverǒglu developed a steady model to optimize the cost of operations of gas pipeline systems(üster and Dilavero?lu，2014);Woldeyohannes &Majid(Woldeyohannes and Majid，2011) analyzed the influence of the compressor parameters setting on the performance of gas grids，under the steady assumption.In general，steady models play an effective role in design，operation optimization and performance analysis of natural gas pipeline systems.However，the consideration of transient flows in the gas pipelines cannot be analyzed from the steady perspective and supply security analysis requires a system-level dynamic model.

Many works are concerned with the transient response of pipelines and pipeline networks，mostly focusing on the improvement of the numerical algorithms used to solve the non-linear hyperbolic partial differential equation (PDE) system that describe them mathematically.Various of numerical methods，e.g.，characteristics method(Elaoud et al.，2017;Sun et al.，2000)，finite difference method (Reddy et al.，2006;Uilhoorn，2017)，finite volume (Xiaojing and Weiguo，2011;Zhang，2016) and finite element method (Bisgaard et al.，1987;Durgut and Leblebiciǒglu，2016)，have been adopted for modeling the transient flows in pipelines.Herrán-González et al.(2009a)developed a method to simulate the transient behaviors of a gas pipeline network by combining the implicit Crank-Nicholson method and the characteristics method.Zhang (2016) used the finite volume method for solving the transient flow problem in a pipeline network.Pambour et al.(2016a)developed an integrated transient hydraulic model to simulate the dynamic behaviors in gas pipeline network systems;the model includes the sub-systems which are important to natural gas supply.

Some unconventional methods have also been developed，such as Control Theory(Alamian et al.，2012;Xiaojing and Weiguo，2011;Zecchin et al.，2009)，analogy method (Ke and Ti，2000) and intelligence algorithms(Madoliat et al.，2016，2017).Tao&Ti(Tao and Ti，1998) developed a dynamic model for gas grids based on the transformation of the gas network to an analogous power grid model including current，voltage，capacitance and resistance.Reddy et al.(2006) developed a transfer function model based on control theory to simulate the transient process in the pipelines.Alamian et al.(2012) developed a state-space model for a simple triangle pipeline network.They used Laplace transformation to simplify the PDEs as a set of transfer functions，which were further transformed to a state space model.Madoliat et al.(2016) developed an intelligent-based method，based on particle swarm optimization gravitational search algorithm (PSOGSA)，to perform the transient simulation in pipeline networks;the accuracy and efficiency have been confirmed in their case study.

The evaluation of complex gas pipeline networks in terms of security of supply requires a system model to simulate the system response under possible scenarios (Kuznetsova et al.，2014).The above works，however，have not given attention to the analysis of the dynamic properties of the gas grids from an overall system perspective.For accounting of the complexities of the system(E Zio，2007) and the constraints from multiple aspects，e.g.，physical properties，technology limitations and market uncertainties，managers and operators need to acquire a comprehensive and clear picture of the global gas pipeline network systems from a supply security perspective.Possible consequences of potential threats(and their combinations)，e.g.，drop of gas productions and sudden increase of demands，should be cautiously studied，and the effects of emergency strategies should be examined to evaluate the resilience of the gas pipeline network.

The requirements of the system dynamic model are stated in detail as follows:

a.The model should be able to describe the responses of the gas pipeline network system under single or multiple disturbances，with acceptable computational burden.

b.The components affecting the capacities of transmission and supply should be comprehensively considered in the model;the important properties of these components should be included.

c.In the model，the constraints of the components should be implemented.

d.The model should be as lean as possible，at the desired level of accuracy and provide the necessary system perspective.

The contribution of this paper is to develop an integrated，statespace formed dynamic model of natural gas pipeline networks.The model includes the components and sub-systems related to the system capacities of supply and transmission，e.g.，pipelines，compressor stations，valves，suppliers，natural gas storages，LNG terminals and customers of gas.Their constraints of operation and the physical limitations are also considered in detail.The individual models of the components of the complex pipeline network integrated via graph theory into one integrated linear state space model.The resulting state-space dynamic model allows reproducing the transient responses of the systems and analyze the consequences of disturbances(e.g.，demand fluctuations)from the system perspectives of stability，observability and controllability.Finally，a supply resilience analysis is developed based on this integrated dynamic model which is used to evaluate the ability of natural gas pipeline networks to maintain stable supply of gas for customers.

2.Model development

The development of the integrated dynamic model of natural gas pipeline network systems is presented in detail in this section(Fig.1).

2.1.Fundamental equations and simplifying assumptions

The gas flow in pipelines is governed by the law of conservation of mass，the law of conservation of momentum，the law of conservation of energy and the real gas law.Considering an infinitesimal control volume in a pipeline with an infinitesimal length dx and a constant cross-sectional area S (shown in Fig.2)，these laws yield the following partial differential equations system (Eqs.(1)-(4))with the assumption that parameters of the gas flow along the pipe are averaged over the cross section area S.

Continuity equation:

Momentum equation:

Energy equation:

Fig.1.The flowchart of the model development.

Fig.2.A control volume in a general gas pipeline.

State equation:

The momentum equation consists of the inertia term，the convective term，the pressure force term，the shear force term and the force of gravity term.In the shear force term，the frictional shear stress τwis related with the dynamic pressure ρv2/2 by the Darcy-Weisbach relation as follows:

The friction factor f is usually estimated by the empirical Coolebrook-white correlation in the condition of turbulent flow(Finnemore and Franzini，2002):

The Reynolds number (Re)，which is the ratio of inertia and frictional forces，is defined as follows:

According to Eq.(6)，f can only be solved iteratively because it is implicit in the equation.In this paper，one explicit approximation，which is valid for turbulent flow in pipelines，is applied to simplify the calculation (Finnemore and Franzini，2002):

Considering the influence of the curvature on effective transportation capacities of the pipelines，the effective friction factor feis calculated based on the efficiency factor φeas follows:

According to the above derivations，the equation system describing the gas flow of pipelines is rather complex and can only be solved by numerical methods.Also，a natural gas pipeline network can contain hundreds of pipelines，nodes and other components:it is impossible to include all the details in the equations system，because of the unaffordable computational burden.Hence，the equations system should be simplified by neglecting some terms，while maintaining the physical accuracy at an appropriate level.Other assumptions generally accepted in previous researches(Pambour et al.，2016b)，are here adopted:

A.The influence of temperature changes on gas flow is negligible.The gas temperature，equal to the environment temperature in which a pipeline is located，is constant in space and time.The environment temperature is correlated to burial depth and ambient temperature.

B.The convective term is ignored on account of its negligible influence compared to the other terms in the momentum conservation equation (Herrán-González et al.，2009b;Pambour et al.，2016a).

2.2.Pipeline modeling

Generally，pipeline flow transient models are developed based on both the mass equation and the momentum equation，and coupling the individual pipeline models by Kirchhoff's first law at the nodes (junctions) (Farzaneh-Gord and Rahbari，2016);this increases the number of equations of the system model and the computation cost.Hence，in this work，only the mass equation and the momentum equation are used to model the nodes and the pipelines.This can significantly reduce the scale of the system models，especially for complex gas grids，and maintain an appropriate accuracy.

On the basis of the assumption A and B in Section 2.1，Eqs.(2)-(4) are simplified as Eq.(10):

Then Eq.(11) is discretized by the finite differential method:

2.3.Node modeling

In this section，we develop the dynamic model of the nodes based on the mass conservation law.The mass equation (Eq.(2))and the real gas equation (Eq.(4)) are simplified based on the assumptions in Section 2.1:

Equation(19)reflects the relation between the pressure change rate and the change of flow.At a node，we assume Δx equals to zero.Therefore，the mass equation at the node j is transformed as follows:

where Qj，nrepresents the gas flow into node j from the connecting pipeline n;when gas flows from the pipeline n to the node j，Qj，nis positive;otherwise Qj，nis negative.Ljrepresents the gas uploaded to (Lj＜0) or downloaded from (Lj＞0) the pipeline network;when the node represents a junction，Lj=0。

Combining Eq.(19)and Eq.(20)，the dynamic model of the node j is as follows:

In order to keep consistent with the pipeline dynamic model in Eq.(15)，Eq.(21) is further transformed as:

wherepj0Qj，n0，Lj0represent the constant values of the corresponding variables at the equilibrium condition of the system.

Finally，the dynamic model of the nodes is:

In the gas pipeline networks，the nodes can mainly be divided in three groups:junctions，suppliers and demand sites.The suppliers can be further divided into gas plants，LNG terminals and underground storages (UGS).These nodes have different characteristics and we need to set constraints on the variables in Eq.(23)according to these differences.

2.4.Control and mathematical modeling

The dynamic behavior and the supply capacity of the natural gas pipeline network systems are controlled mainly by adjusting the operation parameters of the regulation stations (RS)，compressor stations(CS)，gas suppliers(GS)，demand sites(D)and underground storages(UGS).From the mathematical perspective，it is impossible to solve a set of equations when the number of unknowns is larger than the number of equations.Hence，besides the models of pipelines and nodes，the system dynamic model requires additional independent linear equations to close the entire problem.The regulation stations and the compressor stations are controlled by automatically maintaining the operation parameters(e.g.flow rate，inlet pressure and outlet pressure) at the desired set point within the ranges of the constraints.At the points of gas entry and exit，the flow rates and the pressures are typically controlled at desired set points.The additional equations describing the control modes and the constraints of different facilities are listed in Table 1:

Table 1Control modes and constrains of the active components in natural gas pipeline network systems.

Table 2Pipeline data for the triangle pipeline network.

Table 3Input information for the dynamic simulation of the triangle pipeline network.

Besides，we should also notice the specific operation laws of LNG terminals and UGSs:the actual operation of the natural gas storage is carried out by monitoring the operating point (withdrawal/injection rate-working gas inventory)and ensuring that it lies within the storage envelope(Kopustinskas and Praks，2014).The operating region of the LNG terminal is limited by the LNG regasification capacity and working inventory (Ghorbani et al.，2016;McGuire et al.，1986).

2.5.The integrated model

To develop the systematic dynamic mode of the gas pipeline networks，the key is appropriately representing the componentsand sub-systems in one framework and integrating their dynamic models together.Kirchhoff's first law is typically applied to couple the pipelines and the nodes，which increases the model complexity and the computation cost.Considering the network structure of the gas transmission systems，graph theory，which has been applied in power grids (Cadini et al.，2017;Correa and Yusta，2013)，supply chains (Soni et al.，2014)，transportation systems (Mattsson and Jenelius，2015)，etc，is used here as the framework of gas pipeline network modeling.

In this work，a natural gas pipeline network is modeled as a graph consisting of directed edges and nodes.Pipelines，compressor stations and regulation stations are described as edges and junctions，suppliers，UGSs and demands are described as nodes.The edges are further divided into passive components(pipelines) and active components (regulation stations and compressor stations)，because the behaviors of the passive components are described by physical laws and the active components are controlled externally.

The graph topology of the pipeline network is modeled by the following node-branch incidence matrix:

Combining Eq.(23) and the incidence matrix，all the dynamic modes of the nodes and the pipelines are integrated as Eq.(25):

Based on the definition of AI，we need to use the transposed matrix of AI to integrate the edges.Considering the differences between the active elements and the passive elements，the transposed matrix of AI is divided into two parts:

Based on Eqs.(15)-(18) and Eq.(31)，the integrated dynamic model of the pipeline is as follows:

where Kpis the integration of the results of Eqs.(16)and(17)on the basis of BIpipe，i.e.replacing the corresponding elements in the BIpipeby the results of the equations.

In Eq.(32)，the elements corresponding to active components are assigned a null value.

Based on the control modes and the linear equations in Section 2.4，the active components in the system are integrated as follows:

where Cpcontains the coefficients of the linear equations representing the given control modes and the matrix BInonpipe，i.e.replacing the corresponding elements in BInonpipeby the results of the coefficients.The term diag represents a unit diagonal matrix.The elements in the matrix Valuesetare the desired operation values of the active components:

In Eq.(34)，the elements corresponding to active components are assigned a null value.

The control modes of the suppliers，the UGSs and the demands can be directly modeled by adjusting the node dynamic models.When the control mode is selected as the flow control，we only need to adjust ΔL to the desired value;when the pressure control mode is chosen，Δ ˙p is kept equal to 0 at all the time and the initial pressure is the set value.

Fig.3.Flowchart of dynamic simulation.

Fig.4.Layout of the typical triangle gas pipeline network.

Finally，the integrated，systematic gas pipeline dynamic value is developed by combining Eqs.(25)，(32) and (34) as follows:

where m represents the number of edges;n represents the number of nodes;the elements corresponding to the active components in Δ ˙Q and ΔQ are assigned null values.Hence，the system dynamic model is an algebraic-differential equation system.From the Control Theory perspective，the model is a state-space model(Chen and Chen，1984):

where A is the state matrix;X is the state vector;B is the control matrix;U is the control vector.

Several dynamic properties of the gas pipeline network can also be directly derived from the state-space model，such as stability and controllability，allowing for different perspectives of analysis of the complex gas network systems.

2.6.Algorithm

From Eq.(36)，the system dynamic model is an algebraicdifferential equation system.Considering the stiffness of the mathematic problem resulting from the large differences between the absolute changes of flow rates and the pressures，we apply an implicit algorithm with variable-step and variable-order to solve the model and perform the simulation.

The simulation process of the developed system dynamic model is shown in the flowchart of Fig.3:

The space discretization is performed based on the criterion by(Kralik et al.，1988)，as follows:

which means that the elevation between the inlet and the outlet must be less than 200 m and the ratio between the length of the pipeline and its diameter must less than 3 × 104.

3.Case study

3.1.The classical triangle pipeline network

The accuracy of the developed integrated model is verified by benchmarking against the results from other models applied to a typical triangle network(Fig.4).The example pipeline network has been used in many papers(Ke and Ti，2000;Osiadacz and Pienkosz，1988;Pambour et al.，2016b) and three of them are chosen as the benchmarks in this paper.Besides these，the result of TGNET (a professional software)(“PipelineStudio|Energy Solutions，”n.d.)is also used as a benchmark.

The information of the triangle pipeline network is shown in Table 2.The properties are selected from the references in order to compare the simulation results.The initial conditions are calculated based on the steady state by the 4th order Runge-Kutta algorithm.The supplier (Node 1) is pressure-controlled with a constant value of 5 MPa (Fig.5) and the demands are flowcontrolled.The boundary conditions are shown in Fig.6.The dynamic simulation is performed by the variable-step，variable-order method and the computation time is 0.539 s to obtain a converged solution.We also use the 4th order Runge-Kutta algorithm，Adam-Bathforth-Moulton algorithm and Rosenbrock algorithm to perform the simulation，and the computation times are respectively 5.493 s，4.411 s and 5.176 s，which shows that the variablestep，variable-order method algorithm is more efficient for this problem than the others.The computation environment is Inter(R)Xeon(R) (CPU E3-1505M v5 @2.80 GHz) (see Table 3).

The simulation results(from the developed model，the literature references and the TGNET software)of the pressures at the demand sites are compared in Fig.7.The results from the developed integrated model is close to the results of TGNET.There are deviations between the integrated results and the literature references’ results，which could be due to the different method to estimate the compressibility factor and the friction factor.

Fig.5.Pressure condition at supply node (Node 1).

Fig.6.Demands at node 2 &node 3.

The average flow rates in every pipeline and the gas flow from the supply node are compared with the TGNET results.From Figs.8 and 9.，the results of the developed model and the commercial software are very similar.

3.2.Model application to a complex gas pipeline network for the analysis of supply security and resilience

The integrated model is applied for analyzing the system behavior from the supply security perspective.Specially，the model is applied to a relatively complex gas pipeline network，which constitutes a part of a real-world pipeline network.The target gas pipeline network comprises 37 pipelines (total length of approx.1100 km，diameters ranging from 950 mm to 1014 mm)，two pipeline importers，23 demand sites(including city gates，factories，power plants and export stations)，two compressor stations(pressure ratios ranging from 1.02 to 1.18)，seven regulation stations，one UGS and one LNG terminal.The control modes of the two pipeline importers are pressure-controlled while the LNG terminal and the UGS are set at flow-control mode，to analyze their response under different control modes.The regulation stations are set as inactive modes.

In the following，the steady state conditions are used as the initial conditions.As generally crisis of natural gas result from decreases in the capacity of sources and/or abnormal increases of demands，three scenarios(Rodríguez-Gómez et al.，2016;Zeniewski and Bolado-Lavin，2012)(normal condition，supply decrease of the UGS and demand increase of two customers) are simulated and their results are analyzed.In this part，these critical conditions，which can significantly impact the supply security of this natural gas pipeline network，is considered as typical scenarios for verifying the effectiveness of the developed model.

Scenario 1:the system is operated in normal conditions.For convenience of the analysis，we assume that demands fluctuations only occur at two city gates(Customer 15 and Customer 17)and the trends of their demands are presented in Fig.11.This scenario is treated as the benchmark (normal condtions).

Scenario 2:Based on the normal conditions，the capacity of UGS suddenly reduces to 50% at the 8th hour.

Scenario 3:Based on the normal conditions，the demands of Customer 5 and Customer 4 suddenly increase to 150% at the 8th hour.

The impact of the unexpected changes to the normal conditions are analyzed by comparing the flow rate profiles of the pipeline importers (Fig.12) and the pressure of the LNG terminal (Fig.13).The results show that，due to the supply reduction of the UGS and the sudden increase of demands，the pipeline importers increase their supply of gas to maintain the consumptions of gas.The difference in the reaction times，because of the physical flow time of gas，is easily seen in Figs.12 and 13.Furthermore，comparing the behaviors in the scenarios 2 and 3，one can observe the different impacts of the two unexpected events:the gas importer 1 reacts significantly to cover the loss caused by the supply drop of the UGS，whereas the gas importer 2 contributes more to fulfill the sudden increase of the demands because of the configuration of the overall natural gas sources.Although the LNG terminal is far away from the source of disterbances，it also has to adjust its pressure to the changes in the system，to maintain the specified amount of gas supply.

Fig.7.Simulated pressure at Node 2 and Node 3，compared to the results from literature and TGNET (Ke and Ti，2000;Osiadacz and Pienkosz，1988;Pambour et al.，2016b).

Fig.8.Simulated average flow rate in the three pipelines compared to the results from TGNET.

Fig.9.Simulated flow rate of the supply node compared to the results from TGNET.

Figure 14 shows the changes of flow rate and pressure of multijunction node (labeled in Fig.10) and a critical pipeline (the pipeline between the connection node and Customer 11)connecting the gas importer 2 with the other suppliers.These Figures allow analyzing the propagation of the impacts of the disturbance events in the overall pipeline network system and the different behaviors of the components can be understood.For example，in Scenario 2，this multi-junction node suffers a significant drop of operation pressure after the sudden reduction of the UGS output.This pressure drop continues nearly 10 h until the “gap” is fulfilled by the other sources.

Another property which is important for these critical infrastructure is their resilience，i.e.ability to withstand the sudden“strains”and recover to a new equilibrium state(Cadini et al.，2017;Fang et al.，2016;bib_Zio_2016Zio，2016a，b;Zio，2016，2016).In a strained condition，the components of the gas pipeline network will respond.Because of the compressibility of natural gas，the effectiveness of the responses will take some time to become effective for the customers.During this time，the line pack gas plays a critical role of buffering for the robustness of supply.As a disturbance event occurs，the influenced customers lose their normal balance state and take gas from the line pack storage;the consumption slows down to achieve a new balance.Obviously，supply will be assured if the line pack is sufficient for every customer.The rate of line pack(LP)consumption in Eq.(38)，which equals to the load balance(LB)(flow-out gas minus flow-in gas)，is used here as indicator of system recoverability from a degraded state to a new safe supply(new system balance)state.A high level of line pack consumption rate means a severe deviation from the balance:

Fig.10.The complex gas pipeline network system.

Fig.11.Assumed demand fluctuations at Customer 15 and Customer 17.

The resilience analysis in emergency conditions of gas pipeline network is here performed based on the simulation by the integrated dynamic model.The results are shown in Figs.16-19.Besides the three scenarios，another scenario is supplemented in this part to give a clearer illustration.In the 4th scenario，the output pressure of compressor 1# is reduced to 69 bar from 70 bar.

Figure 15 illustrates the load balance evolution of the overall system under different scenarios.From Fig.15，we can observe a significant difference between the profiles of the load balances.The area between the normal and disturbed conditions represents the amount of line pack gas consumed:a large area means a large“l(fā)oss”from the supply security perspective;the pressure reduction event of Scenario 4 turns out to be more severe than the others and should be preferentially prevented.For example，the system resilience in Scenario 4 is much worse than the other scenarios.The load balance in Scenario 4 decreases to nearly -3 Nm3/s3，which indicates the delivery pressure will hit the lower limitation bound soon.

From the Figure，we can also estimate the time required by the system to achieve the new balance.Scenario 4 obviously takes more time to recover than Scenarios 2-3，because Supply 1 is the main source of the system and the pressure reduction significantly reduces its supply ability.

Figures 16-18 show the resilience of different customers in the gas network.From these Figures，we observe that the shapes of the profiles，which represent the recover properties and robustness of the customers，are different.Hence，the securities of different customers are different and we need to compare their resilience to unexpected events and pay more attention to the vulnerable ones.This kind of vulnerability comes from the inherent deficiencies of system design，considering the resources configuration，pipeline capacity and market planning.Besides，considering the different scales of the customers，they cannot be directly compared through their load balance.

Fig.12.Load evolution at the two pipeline importers in the three scenarios.

Fig.13.Pressure evolution at the LNG terminal in the three scenarios.

To overcome this problem，a transformation is performed on Eq.(39).This transformation aims to eliminate the influence of the differences of scales of gas consumptions at different demand sites.The load balance rate (in Eq.(39)) represents the line pack consumption rate per unit of gas demanded to maintain the service level:

The comparisons of the resilience of different customers are shown in Figs.19-21.From Fig.19 we can conclude that Customer 14 is most resilient while Customer 4 is the least.This kind of difference can be observed for most of the customers in the pipeline network.This kind of observations can be caused by different reasons，such as the topological property of the nodes and the configurations of the natural gas sources.By comparing the three Figures，we observe that the resilience performances of the customers also change with the changes of conditions.This is because the consuming rates of the relevant line-pack capacities are quite different under different conditions.Hence，in different emergency conditions，managers and operators should focus on different customers and take appropriate actions to maintain a reliable gas supply.

Fig.14.Pressure and flow rate evolution at the connection node and connection pipeline.

Fig.15.The load balance profile of the overall system in the four scenarios.

Fig.16.The load balance profile of Customer 4 in the four scenarios.

4.Conclusion

Fig.17.The load balance profile of Customer 6 in the four scenarios.

Fig.18.The load balance profile of Customer 14 in the four scenarios.

In this paper，an integrated dynamic model is developed to simulate the system operation and analyze the behaviors of components in complex natural gas pipeline networks，from the supply security and resilience perspectives.In the model，the basic governing equations of pipeline flows and the control mode equations of the important components，e.g.，compressor stations，gas suppliers，UGS and LNG terminals are considered with their technological constraints.The integrated dynamic model is developed by firstly adapting and simplifying these equations through appropriate assumptions and numerical transformation methods，and then integrating them in the network system by graph theory.The developed model is numerically solved by the implicit method with variable-step and variable-order for efficiently addressing the stiffness problem.

The main original contributions of this work are:

Fig.19.The load balance rates in Scenario 2.

Fig.20.The load balance rates in Scenario 3.

Fig.21.The load balance rates in Scenario 4.

(1) We provide a novel method to model a complex natural gas pipeline network as a dynamic state-space model，which allows analyzing system responses to disturbances from the perspective of supply security.The developed state-space model also allows for different perspectives of analysis of the complex gas network systems，such as stability and controllability.

(2) We introduce a method to analyze the supply resilience of complex natural gas networks.The method allows evaluating the capacities of the overall system and customers to withstand disturbances and recover to an equilibrium state.The results of this type of analysis can help managers and operators to find out vulnerable points and take appropriate actions to maintain a reliable gas supply.

For validation，the model has been firstly applied to a typical triangle gas pipeline network，which has been used in several works，and the results have been respectively compared to those of literature and of the commercial software TGNET.Finally，the developed model was applied to a relatively complex gas pipeline network for supply security and resilience analysis.

In future work，the system-level dynamic optimization will be performed based on this integrated simulation model to improve system supply security.Besides，the state-space model will be used to analyze the stability and supply security of complex natural gas network systems from the control theory perspective.

Acknowledgement

This work is supported by National Natural Science Foundation of China[grant number 51904316]，and the research fund provided by China University of Petroleum，Beijing [grant number 2462021YJRC013，2462020YXZZ045].

Appendix.Nomenclature

AI incidence matrix Q gas flow rate in pipeline

BI transposed matrix of AI R natural gas constant

Cp&CQmatrix of control parameters Re Reynold's number

c speed of sound S cross-sectional area

cvisochoric heat capacity t time

D inner diameter T temperature

DS discrete space step Tccritical temperature

diag unit diagonal matrix v natural gas velocity

f friction factor x pipeline coordinate

feeffective friction factor △x pipe segment length

g acceleration of gravity Z compressibility factor

H elevation α inclination

Kp&Kqfactor of Taylor's formular φeefficiency factor

Kp&KQmatrices of Kp&Kqη dynamic viscosity

L loads of nodes (vector)ρ natural gas density

L load of node τwshear stress

l length of pipeline ˙Ω heat exchange rate

p pressure of node