N
 ‐1 security‐constrained coordinated scheduling of integrated electricity and natural gas system considering gas dynamics and wind power uncertainty

Funding information National Key R&D Program of China, Grant/Award Number: 2018YFB0905000; National Natural Science Foundation of China, Grant/Award Numbers: 51877189, 51761135015, U2066601; Young Elite Scientists Sponsorship Program by China Association of Science and Technology, Grant/Award Number: 2018QNRC001 Abstract This paper develops a novel N-1 security-constrained coordinated scheduling model of integrated electricity and natural gas system (IEGS) with consideration of gas dynamics and wind power uncertainty. Direct quantile regression method and Monte Carlo simulation are applied to model the uncertainty associated with wind power forecast in scenarios. Then, the bilinear functions of natural gas constraints are transformed to a set of linear constraints based on McCormick envelopes. To integrate gas dynamics with security analysis in IEGS, a linear dynamic gas flow model considering N-1 contingency of pipeline is established by means of big-M method and the fully implicit finite difference method. Finally, the proposed stochastic coordinated scheduling problem of IEGS is formulated as a mixed-integer linear programming model. Numerical experiments are implemented to verify the superior accuracy of the linear dynamic gas flow model and the effectiveness of the proposed stochastic scheduling model.


INTRODUCTION
With the advantages of flexible operation and convenient regulation, natural-gas fired units (NGUs) can effectively address the variability and uncertainty of renewable energy, and have become a superior choice for the generation expansion in modern power systems [1,2]. It is estimated that the proportion of power generation from natural gas will account for 24% of the global power generation by 2040, leading to an increasing dependence of power system on the natural gas system. Meanwhile, the rapid development of power-to-gas technology blazes a new way for energy storage accelerating a higher penetration of renewable energy [3,4]. Consequently, the coordinated operation of integrated electricity and natural gas systems (IEGSs) plays a more and more critical role in improving the energy utilization efficiency and ensuring the reliability of electricity supply [5].
In recent years, lots of works have been accomplished to study the modelling and scheduling of IEGS. Non-convex Weymouth equation is widely used in the steady-state model of natural gas system [6,9]. A novel mixed-integer linear programming (MILP) approach is formulated in [6] for short-term operation of hybrid natural gas and electricity supply system with consideration of the velocity and compressibility of gas flow. To overcome the challenges of the variability and uncertainty of wind power, a stochastic day-ahead scheduling model is developed with the utilization of Benders decomposition to ensure the feasibility of natural gas system [7]. To enhance computational performance, a modified second-order cone relaxation method and a Newton-Raphson based natural gas correction method are adopted to convexify the Weymouth equation [8]. A stochastic scheduling model of IEGS is proposed to consider the random contingencies of power system [9], where Weymouth gas flow equation is linearized by the piecewise linear approximation method. In addition, the randomness of wind power generation poses severe challenges to the secure operation of IEGS. Probabilistic forecasting is capable of quantifying the uncertainties of wind power generation and providing significant information to security-constrained coordinated scheduling of electricity and gas networks [10][11][12][13][14]. It has been proven that nonparametric probabilistic forecasting which does not rely on the certain probability distribution of forecasting errors is superior to the parametric probabilistic forecasting in reliability and accuracy [15].
The significant difference of travelling velocity between gas and electricity has non-ignorable impacts on unit commitment and dispatch decisions. Thus, it is critical to take into account dynamic characteristics of gas flow. To further study the nodal pressure violation of natural gas system caused by consumption perturbations of NGUs, Laplace transform is adopted to establish the equivalent dynamic model of gas flow [16]. A dynamic scheduling framework of electric-gas systems in different coordination scenarios is presented in [17], of which a novel continuous time optimal power flow model is developed to quantify the real-time gas delivery. By applying Wendroff difference method to approximate the partial differential equations (PDEs) of transient gas flow, a dynamic optimal energy flow model is obtained to achieve the optimal operation of IEGS [18]. In [19], Euler finite difference technique is used to establish a standardized procedure for the dynamics of natural gas transmission networks considering wind power uncertainty. The discretization and linearization process of PDEs has notable impacts on the accuracy of the dynamic gas flow model, which needs to be further studied.
As an important security criterion, N-1 criterion requires that the IEGS be robust enough to hedge against any single failure of components which may lead to the blackouts. Failure of single pipeline may result in the loss of generation because multiple large power generation plants are likely supplied by a single pipeline [20]. As a result, Federal Energy Regulatory Commission has called for necessary study on the unexpected contingency of pipelines in the secure operation of power system [21,22]. In [23], the N-1 contingency analysis for the natural gas system is conducted by using linear sensitivity factors to achieve fast calculation of possible violations. A day ahead economic scheduling model of IEGS is proposed by adopting spinning reserve to address the N-1 generator outages without involving N-1 transmission contingency [8].
This paper develops a MILP model for N-1 securityconstrained coordinated scheduling of IEGS considering gas dynamics and wind power uncertainty. The wind power uncertainty is modelled by a representative scenario set based on Monte Carlo simulation and scenario reduction method. To address the non-convexity and nonlinearity of the momentum equation, McCormick envelopes is utilized to transform the constraints with bilinear terms to a set of linear constraints. Compared with the conventional methods (e.g. ref. [18,28]), this linearization procedure can realize convex relaxation without the rough assumption of fixed gas flow velocity in the day-ahead coordinated scheduling problem. For the purpose of establishing an algebraic dynamic gas flow model, the large time step and limited pipeline subdivision should be adopted to reduce the scale of algebraic equations in the large pipeline networks calculation(e.g. ref. [18,19]), which has great negative impacts on the accuracy of the finite difference methods. To meet this challenge, the fully implicit finite difference (FIFD) method with high accuracy even for large time step and limited pipeline subdivision is adopted to discretize the PDEs in time and space. Subsequently, the linear dynamic gas flow model (LDGF) model is established to integrate N-1 criterion into dynamic natural gas flow model by using the big-M method. Finally, the proposed stochastic model for N-1 security-constrained scheduling problem is formulated as a MILP problem, which can be efficiently solved to obtain the optimal operation strategy of IEGS. The established coordinated scheduling model can be beneficial for integrated gas and electricity network operators such as grid operators to minimize the expected operation cost in the presence of wind power uncertainty [24].
The major contributions of this paper are summarized as follows: 1. A novel MILP model is developed for the N-1 securityconstrained coordinated scheduling of IEGS considering gas dynamics and wind power uncertainty. 2. An algebraic LDGF model with high accuracy is established based on the FIFD method and McCormick envelopes to explore the benefits of gas dynamics for power systems. 3. N-1 criterion of both electricity and natural gas system is integrated into the coordinated scheduling model to guarantee that IEGS is sufficiently robust to suffer from the outage of transmission line and gas pipeline. 4. Probabilistic information involved in wind power scenarios is generated by nonparametric probabilistic forecasting method and Monte Carlo simulation to deal with the renewable uncertainty.
The rest paper is organized as follows. The dynamic gas flow model considering N-1 contingency of gas pipeline is presented in Section 2. Section 3 presents the MILP formulation of the N-1 security-constrained coordinated scheduling model. Numerical simulation results are given in Section 4. Finally, this paper is concluded in Section 5.

Governing equations of natural gas flow
In natural gas system, gas flow driven by pressure travels along the pipeline. The gas states can be determined by the time-dependent continuity equation and momentum equation [17,19]. The momentum equation is expressed as where ρ, π and s represent the density, pressure and flow velocity of natural gas, respectively; g is the gravity acceleration; d and ε are the diameter and friction factor of the gas pipeline, respectively.
The continuity equation is given as where F and A is the mass flow rate and cross-sectional area of gas pipeline, respectively. The state equation can be expressed by where c is the sound speed calculated by where R is the gas constant; Z is the compressibility factor related to the temperature and density of gas; T is the gas temperature. Moreover, with proper assumptions and simplifications [19], the momentum Equation (1) can be simplified to the following form In addition, the mass flow rate F, nodal pressure π and flow velocity s also need to satisfy the following relationship

Linearization of the dynamic gas flow model
The direct application of dynamic gas flow model (2)-(6) will make the optimization problem non-convex, which is complicated to solve. Convex relaxation is an effective way to transform the non-convex functions into convex functions. However, the tightness of relaxation has a significant influence on the solutions of the convex problems. McCormick Envelope is a crucial relaxation method that not only guarantees convexity, but also has sufficiently tight bounds [25]. Thus, McCormick Envelope is employed to linearize the dynamic gas flow model (2)- (6) in this paper.
In Equations (5) and (6), there are two bilinear terms F ⋅ s ands ⋅ . To convexify and linearize Equations (5) and (6), a continuous variable is introduced to replace the bilinear term F ⋅ s in Equation (5), and the momentum Equation (5) becomes The lower bounds of the variable are given as follows, where F min and F max are the minimum and maximum mass flow rate of gas pipeline; s min and s max are the minimum and maximum flow velocity. The upper bounds of the variable are expressed as: By replacing the bilinear term s ⋅ in Equation (6) with another continuous variable , Equation (6) becomes Similar to the definition of variable , the lower and upper bounds of variable are given by where min and max are the minimum and maximum nodal pressure respectively. After the above mathematical manipulations, bilinear Equations (5) and (6) can be replaced by the linear constraints (7)- (16). However, the continuity Equation (2) and momentum Equation (7) are still PDEs which are extremely difficult to calculate.

Differencing form of the dynamic gas flow model
To establish a solvable model of natural gas transmission network, the FIFD method [27] is adopted to transform the PDEs to algebraic equations. Compared with partially implicit finite difference method (e.g. the Crank-Nicolson method [29]), the FIFD method is demonstrated to be unconditionally stable by Von Neumann stability analysis [27]. The partial derivative with respect to time t in Equation (2) and (7) can be approximated by where H represents the variables including ρ, π, F and in Equations (2) and (7); Δt is the time resolution; the superscript t and t+1 denote the time level; the subscript i and j represent the inlet and outlet of each gas pipeline. The partial derivatives with respect to space x can be approximated by the centreddifference form, expressed as Moreover, the individual term H at the junction between two nodes is replaced by Replacing the spatial step Δx in Equation (18) by the length of the pipeline l i j and substituting Equations (17)- (19) into Equations (2) and (7), the following linear differencing forms (20) and (21) can be derived.

Boundary conditions of natural gas Network
In the large-scale natural gas network, there exist many types of non-pipe elements including gas source, compressor station, etc. At source nodes, the pressure is considered as a constant 0 s , shown as 0 s = s ⋅ t (22) where s is the coefficient matrix of gas source, the element sg in the coefficient matrix s equals to "1" if the gas source located at node g, else sg =0; t represents the pressure of all nodes at time t; 0 s is the constant pressure of gas source nodes. For sink nodes, the mass flow rate F is also a constant value, expressed as where l is the coefficient matrix of sink nodes, the element li in the coefficient matrix l equals to "1" if node i is a sink node, else li =0; F t denotes the mass flow rate of all nodes at time t; F 0 l is the mass flow rate of all sink nodes. For the junctions of pipelines, an incidence matrix B is introduced to describe the relationship among pipelines. Element "0" in incidence matrix B indicates the junction does not belong to the pipeline, while "−1" ("1") denotes the junction is the inlet (outlet) of the pipeline. Incidence matrix B can be represented as the sum of the two matrices B + and B -: where B + is all the positive elements in matrix B; Bis all the negative elements in matrix B. Then, the first Kirchhoff's law at the junctions of pipelines can be expressed as: where F inlet t and F outlet t represent the mass flow rate of all inlet nodes and outlet nodes at time t, respectively.

Dynamic gas flow model considering N-1 pipeline contingency
In USA, there are about 516,599 km onshore and offshore gas pipelines, it is reported that the incident rates of natural gas pipeline range from 0.015% to 0.3% (per km year) [33]. Both natural disasters (e.g. earthquakes, extreme frosts etc.) and human factors (e.g. bad maintenance operations, land excavations etc.) can cause the failure of gas pipelines [23]. The redistribution of nodal pressure caused by the outage of gas pipeline may lead to abnormal disconnection of gas-fired power plants and blackout. However, the slow dynamics of gas flow can effectively mitigate the impacts of the various contingencies and fluctuations in energy demands as well as renewable energy, which can enlarge the feasible operation region of IEGS as shown in Figure 1.
To integrate N-1 contingency of gas pipeline into the optimization model, a contingency matrix n is introduced to represent the operation states of all gas pipelines in natural gas networks. The element n t p in n is utilized to indicate whether the gas pipeline is in service, where n t p = 1 represents the normal operation state of gas pipeline p, and n t p = 0 represents the contingency state of gas pipeline p. The first column K n indicates that all gas pipeline are operating in a normal state, and the column K c means that one of the gas pipelines is in the contingency state.
The big-M method is employed to equivalently transform the continuity Equation (20) and the momentum Equation (21) from algebraic equations to inequality constraints. The continuity Equation (20) is formulated as The momentum Equation (21) is transformed to where LM G is the "big M" value. When pipeline p is under normal operation state (n t +1 p = 1), inequalities (27)-(30) degrade to equality constraints; when pipeline p is under contingency operation state (n t +1 p = 0), LM G is large enough to guarantee (27)- (30) are satisfied regardless of the difference of mass flow rate and pressure between two nodes. In addition, mass flow rate F t +1 i of node i should also subject to the capacity limits of pipeline p given by where F p,min and F p,max are the minimum and maximum mass flow rate of pipeline p, respectively. In addition, gas pressure m of node m should subject to where min m and max m are the minimum and maximum pressure of node m, respectively.
Since the proposed LDGF model is related to both time t and space x, the steady-state values of gas pressure , mass flow rate F and flow velocity s at t = 0 calculated by Weymouth equation [18] are considered as the initial condition.

3.1
Objective function Power transmission system and natural gas transmission system can be unified as an integrated energy system, in which multi energy infrastructures are jointly operated. The objective of the N-1 security-constrained coordinated scheduling problem is to minimize the expected total operation cost, expressed as where i,t and p j ,t denote the cost coefficients of gas sources and cost coefficients of coal-fired generators, respectively; z represents the probability of occurrence of wind power generation scenario z. GSU z k,t and GSD z k,t are the start-up and shutdown binary variables for the generation unit i at time t, respectively.

Operation constraints of electricity network considering N-1 transmission contingency
In the electricity network, the power demand is satisfied by coalfired, gas-fired, and renewable energy generation, expressed as (34) where P P2G,z b,t is the power consumption of P2G unit b; P c,z i,t and P g,z j ,t denote the power generation of coal-fired generator i and NGU j respectively; P r,z e,t is the power output of renewable resource e.
Failure of transmission line is also a major cause of blackout. It is reported that the values of power transmission line failure rates range from 0.1% to 0.9% (per km year) in Canada [32]. Thus, it is significant to include N-1 transmission contingency in the operation of power system. To take into account N-1 transmission contingency, binary variable where LM p is a number with big value to ensure that the inequalities (35) and (36) (50) and (51) [9].
where P c i,min and P c i,max are the minimum and maximum outputs of coal-fired generation unit i respectively; U In deterministic approach, spinning reserve is required in power system to withstand the uncertainty associated with wind power generation [24]. The capacity of reserve is restricted by where R t is the total scheduled reserve; R min is the minimum reserve requirement.

Modelling of wind power uncertainty
The variability and uncertainty of wind power pose severe challenges to the secure operation of power systems. Probabilistic forecasting taking the form of predictive quantiles is capable of quantifying the uncertainties of wind power generation [15]. Denote the probability density function and cumulative distribution function of wind power generation by w t and W t respectively. The corresponding quantile q t is defined by (53)-(54).
Pr(x t ≤ q t ) = (53) where Pr() represents the probability operator; x t is the random wind power generation; α is the nominal proportion of the quantile q t . Then, the probability density function w t and cumulative distribution function W t of wind power generation can be approximated by a set of predictive quantiles, expressed aŝ whereq i t +k|t denotes the estimation of actual quantile q t +k|t ; W t +k|t is the predictive quantile series with proportions from zero to one.
The predictive quantile series W t +k|t of wind power generation can be efficiently obtained via the direct quantile regression method [15]. The direct quantile regression method that takes the advantages of both quantile regression and extreme learning machine can formulate the probabilistic forecasting as a linear programming problem.
Given the predictive quantile series W t +k|t , numerous scenarios are generated to model the uncertainty of wind power via Monte Carlo simulation [31]. The scenario-based stochastic optimization approach is adopted to take into account the wind power forecasting uncertainty in the coordinated scheduling of IEGS. However, scenario-based stochastic scheduling with a large number of scenarios suffers from heavy computational burdens. It is essential to reduce the complete scenario set to several representative forecast scenarios via an effective scenario reduction method [26]. In this paper, SCENRED provided by the General Algebraic Modelling System is utilized as a tool for the scenarios reduction [30].

Coupling constraints of IEGS
In IEGS, power system and natural gas system are interconnected by NGUs and P2Gs. Natural gas can be converted into electricity by NGUs, and electricity can be converted into natural gas by P2Gs. In this study, the energy curves of NGUs and P2Gs are described by two linear functions, expressed as: where g f and P2G represent the energy conversion efficiency coefficients of NGUs and P2Gs respectiveely.
Meanwhile, the capacity limitation of NGUs and P2Gs are given as where P g max and P P2G max are the maximum rated power of NGUs and P2Gs respectively.

Summary of the stochastic scheduling problem
The flowchart of the stochastic N-1 security-constrained scheduling problem is presented in Figure 2. The 24-h ahead predictive quantiles of wind power generation can be produced via direct quantile regression method. Then, numerous scenarios are generated by Monte Carlo simulation to account for the uncertainty of wind power and the scenario reduction method is employed to reduce those scenarios to several representatives. Based on the generated representative scenarios, the probability information of wind power forecasting uncertainty can be integrated into the coordinated scheduling model. Then, the PDEs governing the dynamic gas flow are transformed into a group of algebraic linear equations and inequalities with the aid of FIFD method and McCormick envelopes. Coupling the DC power flow model with the established LDGF model, this coordinated scheduling problem of IEGS can be formulated as a MILP model.
In the established MILP model, N-1 contingencies of power transmission line and gas pipeline are considered, and gas pipelines are assumed to completely out of service under contingency state. Despite of large-scale variables and constraints, the proposed MILP model can be easily solved by mature optimization software, such as CPLEX.

System description
In the study, the test IEGS is constructed by combining the IEEE 24-bus power system and a 25-node natural gas system. As shown in Figure 3, the test system includes 7 coal-fired generators, 2 NGUs, 1 P2G, 1 wind farm, 35 transmission lines, 2 gas sources and 26 pipelines. The hourly forecast data of electricity load, natural gas load and renewable power is depicted in Figure 4. The maximum active power flow of each electricity transmission line is 250MW and the minimum pressure of the gas node is 4Mpa. Big-M value is set to 10 7 . The test data of natural gas system is given in [34]. Direct quantile regression approach [15] is employed to produce the 24-hour ahead predictive quantiles of wind power generation. All experiments are

Effectiveness of LDGF model
In this section, a single pipeline of natural gas network with 80 km in length and 0.6 m in diameter is utilized to verify the effectiveness of the proposed LDGF model. The inlet pressure is set as 0.2 MPa. The outlet pressure and the mass flow rate of inlet are variables to be determined. Time resolution is set as 5 min.
1) Accuracy analysis of linearization method: To verify the accuracy of the proposed linearization method (7)-(16), a relaxation error V o is introduced to represent the quality of solutions, defined as where o represents the values of mass flow rate F and pressure π estimated by the proposed linearization method;ō is the values of mass flow rate F and pressure π assessed by the nonlinear model (2)-(6) accordingly.
1) The proposed linearization method is compared with the average velocity method widely used in [18] and [28]. The maximum linearization errors in solutions obtained from different methods are listed in Table 1. The maximum errors of pressure and mass flow rate assessed by the proposed linearization method are 3.26% and 2.11%, which are the Wendorff difference-based gas model [18] and Euler difference-based gas model [19] are employed as benchmarks. Figure 5 compares the outlet pressure estimated by FIFD method and benchmarks with the solutions of PDEs and steady-state model. It can be seen that the performance of PDEs is significantly different from that of steady-state model without considering gas dynamics. Meanwhile, the LDGF model has better performance with respect to accuracy than benchmarks and its maximum error is only 3.7% compared to the PDEs, demonstrating its superior approximation to PDEs. 3) Sensitivity analysis of time resolution to LDGF model: Sensitivity analysis of the time resolution to the LDGF model is presented in Figure 6. It can be easily found the proposed LDGF model performs well and obtains the minimum error when the time resolution is selected as 5 min. With the time resolution increasing from 5 to 60 min, the performance of LDGF model becomes slightly worse, but the solution is still acceptable due to the satisfactory approximation of most critical dynamic characteristics of the PDE simulation curve. Therefore, the dynamic characteristics of gas flow can be well described when the proposed LDGF model is applied to the coordinated scheduling problems with time resolution of 1 hour.
1) Impacts of gas dynamics on IEGS scheduling: The response time of gas flow within a pipeline ranges from tens of minutes to hours, which is determined by the parameters of pipelines including length, topology and diameter [27]. The response time difference of power and gas flow has significant impacts on the coordinated scheduling results of IEGS.
The scheduling results of different cases are listed in Table 2. It can be seen that the total generation cost of Case 1 is higher than Case 2. Considering the gas dynamics effectively expands the feasible operation region of coordinated scheduling, and more conservative scheduling results are obtained with steady-state based gas flow model. The mass flow rates of source node 1 under Case 1 and Case 2 are compared in Figure 7. It can be seen that the drastic changes in mass flow rate are significantly alleviated at hours 3-5, 6-9 and 11-14. Benefiting from the linepack characteristic of gas pipeline, the natural gas is stored in the pipeline during valley periods and released during peak periods, which plays the role of peak clipping and valley filling. Thus, gas dynamic is of great significance to improve the operation flexibility of IEGS. The nodal pressures of NGU 2 in different cases are depicted in Figure 8. With the consideration of N-1 contingency, the normal operation pressure of NGU 2 is 2) Effectiveness of scenario-based stochastic approach: Due to the randomness and variability of wind power, the system becomes more vulnerable to various contingencies. Stochastic conditions are introduced to Case 3 for the comparison between the deterministic and the stochastic solutions. Based on the predictive quantile series of wind power, 2000 stochastic scenarios are accordingly generated by Monte Carlo simulation, and the scenario reduction method is utilized to reduce those scenarios to 10 representatives. The weights and total operation costs of each representative scenario are listed in Table 3. The generation cost varies in scenarios because each representative scenario after reduction indicates a combination of random wind power outputs within 24 hours. The expected generation cost of Case 3 is $ 1,459,486, which is lower than that of Case 2. With probability information involved in wind power scenarios provided to the coordinated scheduling model, the operation strategy with enhanced economics can be obtained.
The computation time of each case is given in Table 4. Compared with the steady-state gas flow model, the dynamic gas flow model includes not only spacedependent variables, but also time-dependent variables. By considering gas dynamics in the coordinated scheduling problem, the computation time of Case 2 is longer than that of Case 1 scheduling with steady state gas flow constraints. Moreover, by employing the scenario-based stochastic approach to deal with wind power uncertainty, variables associated with scenarios are introduced into the scheduling problem and the computation time increases prominently.
3) N-1 contingency analysis: To further verify the validity of the proposed model against possible N-1 transmission contingency, the following three additional cases are considered.
Case 4: No contingency is considered in Case 3. Case 5: Only N-1 gas pipeline contingency is considered in Case 3. Case 6: Only N-1 power transmission line contingency is considered in Case 3. Table 5 lists the scheduling results of additional Cases 4-6. Obviously, the operation cost of IEGS increases with the integration of N-1 criterion into the coordinated scheduling model, because more operation cost will be produced to deal with the vulnerability of the IEGS to N-1 contingencies. By considering both N-1 contingencies, the operation strategy of Case 3 is the most conservative and generates the highest operation cost.
The scheduling results of NGU 1 with different contingency constraints are compared in Figure 9. For the N-1 criterion of power and natural gas system is not violated during valley periods of demand, the active power output of NGU 1 remains unchanged under different contingency scenarios. Compared with the scheduling results with no contingency, the output of

FIGURE 9
Scheduling results of NGU 1 with different contingency constraints

FIGURE 10
Active power flow through transmission lines in the different N-1 contingencies at T = 9h NGU 1 is restricted by power N-1 criterion during power peak at hours 9-11 and 18-20, while the generation of NGU 1 is limited by gas N-1 criterion during gas peak at hours 8-10 and 12-13. By considering both N-1 contingency, the scheduling problem has the most security constraints and NGU 1 has the minimum total power generation. Due to the tight coupling of power and gas systems, both the N-1 criterion in power and gas systems can significantly affect the scheduling results of NGUs. Figure 10 shows the active power flow through each transmission line in the different N-1 contingency scenarios. Compared with the normal operation state of the power system, the active power flow of each line is redistributed due to the outage of transmission line. Moreover, in the single contingency scenarios of line L11 and line L26, the active power flows of line L12 and line L32 approach to its upper and lower limits respectively, which demonstrates the effectiveness of the proposed MILP model against N-1 transmission contingency in power system. However, the N-1 contingency in natural gas system is more complicated than that in power system due to the distinct response time of power and gas flow. It is assumed that the con-  tingency takes place at T = 8 h and is overhauled at T = 12 h. The pressure of node 5 under the contingency scenario of pipeline 2-3 with and without considering the dynamic operation constraints of gas network are compared and depicted in Figure 11. It can be observed that the failure of pipeline 2-3 leads to a significant drop in node pressure, which seriously threatens the secure operation of gas system. Moreover, the pressure of node 5 involving steady-state network constraints surpasses the lower limit, while the nodal pressure is always within the security region under dynamic natural gas constraints. The proposed dynamic gas flow model can accurately model the redistribution of nodal pressure after the failure of gas pipeline and avoid excessive defence against the contingency.

Verification on larger system
To further validate the effectiveness of the proposed stochastic N-1 security-constrained coordinated scheduling model, a larger case [8] consists of modified IEEE 118-bus system with 40-node gas network is tested. The topology of the large test system is shown in Figure 12. The total operation costs of deterministic cases are presented in Table 6. With gas dynamics considered in Case 2, the production cost significantly reduced from $ 4,954,788 to $ 4,922,309. The compressibility and storability of natural gas make security constraints of IEGS less stringent and provide more flexibility to coordinated operation. In this test system, the capacity of NGUs accounts for 17% of total power generation capacity. Available flexibility will increase with the growth of installed capacity of natural gas fired units.
Occurrence probabilities and scheduling results for different scenarios in Case 3 are listed in Table 7. The expected total operation cost for the ten scenarios in Case 3 is $ 4,907,064, which is lower than that of Cases 1 and 2. The computational time of different cases is proposed in Table 8. By introducing stochastic conditions, the IEGS can operate more economically, but computation time of the scheduling model increases accordingly.
By considering N-1 criterion, the IEGS can continue to operate even if one component is out of service. In the test system, the outage of pipeline 4-14 cannot impose huge impacts on IEGS, but the failure of pipeline 31-32 can result in severe redistribution of mass flow rate and pressure. Therefore, identifying the most critical contingencies and decreasing the number of monitored contingencies are of great significance to the economics and reliability of IEGS.

CONCLUSION
The increasing dependence between power and gas systems significantly highlights the importance of the IEGS co-optimization issues. With consideration of gas dynamics and wind uncertainty, this paper proposes a novel stochastic MILP model to solve the N-1 security-constrained coordinated scheduling of IEGS and obtain the optimal operation strategy. Numerical simulation results validate the merits of the proposed scenario-based stochastic approach in dealing with uncertainties involved in wind power. The transient process of the gas flow can be well modelled by the established LDGF model to evaluate the operation status of IEGS more accurately. After considering the N-1 contingency in both power and natural gas systems, the operation strategy of IEGS is more defensive and can better prevent the huge impacts caused by component failures. Off time of unit i for period t R t The total scheduled reserve in power system for period t