A Robust Formulation Model for Multi-Period Failure Restoration Problems in Integrated Energy Systems

: The risks faced by modern energy systems are increasing, primarily caused by natural disasters. As a new form of multi-level energy complimentary utilization, integrated energy systems are attracting more and more attention for their high-e ﬃ ciency and low-cost. However, due to the deep coupling relationship between systems, they are more susceptible to natural disasters, resulting in a cascading failure. To enhance the resilience of the integrated electricity-gas system, this paper proposes a failure restoration strategy after a natural disaster occurs. First, the temporal constraints of the dispatching model are considered, and the failure restoration problem is molded into a multi-period mixed-integer linear programme, aiming to recover the interrupted loads as much as possible. Second, since the uncertain output of distributed generation sources (DGs) such as wind turbines and photovoltaic systems will threat the reliability of restoration results, the robust formulation model is incorporated to cope with this problem. Third, we propose a new modeling method for radial topology constraints towards failure restoration. Moreover, the Column and Constraints Generation (C&CG) decomposition method is utilized to solve the robust model. Then, the piecewise linearization technique and the linear DistFlow equations are utilized to eliminate the nonlinear terms, providing a model that could be easily solved by an o ﬀ -shelf commercial solver. The obtained results include the sequence of line / pipeline switchgear actions, the time-series dispatching results of electricity storage system, gas storage system, and the coupling devices including the gas-ﬁred turbine, power to gas equipment. Finally, the e ﬀ ectiveness of the proposed restoration strategy is veriﬁed by numerical simulation on a 13-6 node integrated energy system.


Introduction
In recent years, with the development of the integration of renewable energy and the growing demand for environmentally-friendly operation, integrated energy systems are believed to be a promising way to satisfy energy customers' demand with high-reliability, high-efficiency and low-cost [1][2][3][4]. Among the different kinds of energy systems, the natural gas system and electric power system play essential roles as the most common options for complementary operation as an integrated electricity-gas system (IEGS) [5][6][7][8]. Natural gas consumption by the electric utilities in the U.S. increased by 60.3% from 2008 to 2018. At the same time, 35.1% of electricity generation came from natural gas in 2018 [9]. The co-optimization of the electric power system and the natural gas system has received more and more attention. There is a bi-directional energy flow between the electric network and natural gas network through the coupling devices, the gas-fired turbines (GFT) and the power-to-gas (P2G) [10][11][12]. In this way, the IEGS can achieve economical operation with high efficiency.
Especially in highly coupled IEGS, there are many electric lines, natural gas pipelines, DGs, etc. It is challenging to co-optimize these devices over a multi-period. Some literature has studied the sequential process of failure restoration in electric distribution networks. In [33], the sequential actions of switchgear is generated and several microgrids are formed to recover the interrupted loads in the unbalanced electric distribution network. Based on the genetic algorithm, [34] proposes a two-stage failure restoration method for getting the feasible topology in the first stage, generating the optimal sequence actions of switchgear in the second stage.
Meanwhile, the DGs in the distribution power system can provide emergency power supplies during the process of failure restoration [35][36][37] and improve the resilience [38]. However, due to the intermittency and uncertainty of DGs, it is challenging to co-optimize the multiple DGs during a failure restoration process. The information gap decision theory was utilized to cope with the uncertainty in failure restoration by Chen [39], and the obtained robust solution can guarantee that the restoration results are not worse than a specific value. A two-stage robust failure restoration model was established by Chen [40], and the restoration strategy can guarantee that the security constraints will not be violated in the worst scenario of uncertainty.
However, as far as the sequential failure restoration method for an IEGS considering the uncertainty of DGs and load demands is concerned, there are relatively few studies available. To fill this gap, this paper proposes a robust multi-period failure restoration method for the IEGS by co-optimizing DGs, battery storage system (BSS), gas storage systems (GSS), GFT, P2G and line/pipeline switchgear. To achieve a full-process, time-series failure restoration, the problem is formulated into a mixed-integer linear programming model which is tractable to solve since the linear DistFlow [41] equation and the piecewise linearization technique are utilized. The effectiveness of the proposed method is verified by simulation on a 13-6-node IEGS which is composed of a 13-node electric network and a 6-node natural gas network. The main contributions of this paper may be summarized as follows: (1) Proposing a three-stage failure restoration process, including the action sequence of switchgear, so in this way the feasibility of restoration results can be ensured; (2) Modeling the failure restoration problem of IEGS as a multi-period mixed integer linear problem, considering the temporal characteristics of various devices, including DGs, BSS, GSS and line/pipeline switchgear; (3) Considering the bi-directional coupling of IEGS, the GFT and P2G are coordinated to achieve the failure restoration, and the reliability of the IEGS is improved; (4) Proposing a new modeling method for radial topology constraints towards failure restoration; (5) By introducing the linearization technique, the nonlinear mathematical model of the failure restoration is formulated into a mixed-integer linear programming problem, which can be solved by an off-the-shelf commercial solver; (6) The Column & Constraint Generation (C&CG) [42] algorithm is used to solve the three-layer robust model. A modified 13-6-node IEGS is used as a case to verify the effectiveness of the proposed model.
The rest of the paper is organized as follows: Section 2 describes the deterministic model formulation of the failure restoration for an IEGS. A robust failure restoration model is then proposed in Section 3. Section 4 describes the case study of an IEGS composed of a 13-node electric network and a 6-node natural gas network. Finally, Section 5 summarizes the paper and discusses possible future work.

Introduction to Failure Restoration
First, to emphasize the importance of the temporal characteristics of failure restoration, the several stages that happen when the IEGS suffers from a fault need to be analyzed. Traditionally in an electric power system, failure restoration studies mainly focus on the two-stage failure restoration process, as shown in Figure 1a, including the faulty stage and the restored stage, but in an IEGS, due to the tightly coupled nature, the two-stage failure restoration process cannot represent the temporal characteristics of the topology, DGs, BSS and GSS, which may lead to unfeasible operation results in the restoration process. To this end, this paper applies the three-stage failure restoration process for the IEGS. As shown in Figure 1b, there is one more stage, the failover stage, which includes the dynamic process of topology, the output of DGs, BSS, and GSS. Second, the structure of the IEGS needs to be introduced. The electric network is bidirectionally coupled with the natural gas network through gas-fired turbine (GFT) and power-to-gas (P2G) devices. The energy management system (EMS) of IEGS can achieve the information acquisition and remote dispatching through the remote terminal unit (RTU), and the controllable equipment such as GFT, P2G, BSS, and GSS are dispatched to achieve the reliable and economical operation. The software function of failure restoration can be integrated into the energy management system of IEGS. After a large-scale outage accident (such as earthquakes, typhoons, and other harsh natural disasters) happens, neither the electric nor gas transmission systems can supply energy to the corresponding distribution systems. At the same time, the electric distribution network and the gas distribution network also suffer inevitable failures. By dispatching the output of GFT, P2G, BSS, GSS, and actions of switchgear in the electric network and the action of the pipeline valves in the natural gas network, the dispatching center can quickly restore the electric power supply and natural gas supply after the large-scale outage happens. Besides, the GFT can transform natural gas into the electricity for supporting the electric network. Conversely, P2G can supply natural gas networks by consuming electricity. Through both approaches, the purpose of maximizing the restored load can be achieved. Figure 2 shows the flowchart of the restoration process.

Objective of the Failure Restoration Model
In order to simplify the modeling of failure restoration and guarantee the feasibility of the problem, it is necessary to make the following assumptions: (1) The topology of the electric network and the natural gas network maintain a radial operation, in which the energy flow is bi-directional; (2) The transmission network cannot supply energy to the distributed electric and natural gas network. The IEGS needs to perform a local failure restoration. It is assumed that the IEGS can maintain the islanded operation before the fault is repaired. The goal of failure restoration is to maximize the interrupted load, including electrical and natural gas loads. It can be formulated by the following mathematical forms: where: the electrical load P D i,t and the natural gas load F D n,t are multiplied by the coefficients β L i and β L n , respectively, representing the priority of different loads. ∆t indicates the action period of each operation step, one hour is taken [31]. Subscript t denotes serial time. T is the set of time which belongs to the whole restoration period. Ω e and Ω g represent the set of nodes in the electric network and natural gas network, respectively. The above mathematical model is a mixed-integer nonlinear programming problem. In this paper, the nonlinear terms are handled with a linearization method. The constraints in (1) will be described separately in the following sections.

Modelling of the Electric Network
The linear DistFlow equation has been widely used in the power flow of radial distribution networks and microgrid systems, as shown in Figure 3. The equations are relatively easy to solve because the non-linear terms in the original power flow equation are neglected. Based on the DistFlow equation, this paper introduces the binary variables of each line, and the modified power flow constraints of the electric network are obtained. For each node in the electric network, the power flowing in and out should be equal: where P ki,t represents the active power flow in the line (k,i) which connects the nodes k and i. P dg i,t is the generated active power of DGs in the node i at time t. P D i,t denotes the active power consumption in the node i at time t. v i,t stands for the status of restoration in node i at time t, node is restored when v i,t = 1, otherwise v i,t = 0. Ω elc l denotes the set of lines in the electric network. (3) is the corresponding reactive power equation.
For each line in the electric network, the nodal voltages of both sides have a relationship with the active and reactive power flow of the corresponding line. The restoration status variables are also introduced: where Equation (4) represents the relationship between the voltage drop and the power flow in each line. V i,t denotes the amplitude of voltage in node i. R ij , X ij are the resistance, reactance of line (i,j), respectively. x ij,t represents the restoration status of line (i,j). b ij,t is the auxiliary variable of line (i,j). When the line (i,j) is energized, x ij,t = 1 From the constraint (5) we know the corresponding b ij,t = 0, thus the voltage equation (4) is feasible. On the contrary, when the line For each line in the electric network, there is a capacity limit: where P max ij is the limit of active power flow in the line (i,j). The sign of negative and positive denotes the bi-direction of power flow. When line (i,j) is energized, x ij,t = 1, the capacity constraints of lines are active. Otherwise, the power flow of lines is zero. (7) is the corresponding reactive constraint of line (i,j).
For each nodal bus in the electric network, there is a voltage limit: where V min , V max are the minimum and maximum limit of the voltage of each node. 0.9 and 1.1 are taken, respectively [43]. A battery storage system (BSS) can assist in handling the intermittent problem of DGs and load demands which has an inter-temporal dispatching requirement. The main concern of BSS dispatching is to control the power of charging or discharging, and the remaining SOC. The mathematical model of the BSS adopted in this paper does not consider the specific battery model and parameters, but only the general modeling, the mathematical form is as follows: where α c , α d are the charging and discharging status variable of BSS, respectively. α c = 1, α d = 1 indicates that the BSS is in charging and discharging, respectively. P c bss,t , P d bss,t denotes the active power of charging and discharging, respectively. ϕ bss t is the status of charging of BSS at time t. η c , η d are the efficiency of charging and discharging, respectively. For simplicity, it is assumed that the efficiency is not affected by the charging power. Besides, enough reactive capacity is assumed.

Modelling of the Natural Gas Network
The natural gas network comprises a gas source, pipelines, compressor, storage and load ( Figure 4). For simplicity, the steady natural gas flow formulation is utilized, the pipeline storage and dynamic characters are ignored [44]. Meanwhile, the radial operation topology is assumed. With these assumptions, the constraints of nodal balance and pipeline flow in a natural gas network could be formulated as follows: where F ki,t represents the gas flow of pipeline (k,i) at time t. F D i,t is the gas consumption of node i at time t. F P2G i,t stands for the output amount of gas of P2G. F P2G i,t stands for the gas consumption of GFT. Ω gas l is the set of pipelines in the natural gas network, Ω g is the set of nodes in the natural gas network. The natural gas pipelines can be divided into two categories: pipelines with compressors and pipelines without compressors. For the pipelines with compressors, the gas flow constraints can be expressed as follows: where π i,t represents the pressure of node i at time t. x ij,t is the binary variable, representing the status of restoration in the pipeline (i,j). The pipeline is restored when x ij,t = 1, otherwise x ij,t = 0. λ ij is the parameter of the compressor in the pipeline (i,j), corresponding to different operating conditions. The compressor could be illustrated as Figure 5, energy consumption of compressor is ignored just for simplification. At the same time, for natural gas pipelines without compressors, the gas flow has the following relationship with node pressure: where Φ indicates the comprehensive parameter of the pipeline [45]. sgn ij,t (π i,t , π j,t ) represents the direction of the gas flow in the pipelines. Because of the existence of quadratic terms in the Equation (17), there is a need to eliminate these nonlinear terms. Squaring both sides of (17), letting µ i,t = π 2 i,t . Then the auxiliary variables δ, η are introduced to do the incremental linearization [46,47]: Equation (19) is a piecewise linear form of a square function, (20) is a piecewise linear form of independent variables. Equation (21) are constraints of the introduced variable corresponding to piecewise positions and quantities, k is piecewise number, and N is the total piecewise quantities. After applied with the above incremental linearization, the non-linear term in the (17) can be eliminated.
The GSS can provide standby natural gas supply after system failure to support short-term failure restoration. Similar to BSS, the model of GSS can be formulated as follows: where S i,t indicates the remaining gas capacity of GSS in node i at time t. S min i , S max i represent the minimum and maximum gas capacity limit of GSS, respectively. F in i,t , F out i,t indicate the input and output gas flow of GSS in node i at time t.
The nodal pressure of natural gas network needs to be within a reasonable range in order to ensure safe operation. The corresponding constraints are as follows: where π min ,π max are the minimum and maximum limit of pressure in the natural gas network, respectively. IEGS achieves bi-directional coupling through GFT and P2G. By opening the P2G device, surplus electricity can be converted into natural gas and injected into the natural gas network for storage and transportation. The GFT can consume the natural gas for generating the electricity into the electric network: where F GFT i,t represents the gas consumption of GFT in node i at time t. P GFT i,t is the output active power of GFT in node i at time t. G represents the low calorific value of natural gas (35590 kJ/m 3 ) is taken from [48]. Ω GFTU denotes the set of GFT units.
where F P2G i,t represents the output natural gas of P2G in node i at time t. P P2G i,t is the active power consumption of P2G in node i at time t. Ω P2GU denotes the set of P2G units.

Topology Constraints
In the whole process of failure restoration, the electric network involves the operation sequence of line switches, load switches, BSS, etc. The natural gas network involves the operation sequence of pipeline valves, GFT and P2G. Topology constraints ensure that the IEGS can maintain the radial topology in each step of the entire failure restoration. It should be emphasized that the essential difference between the reconfiguration and failure restoration problem modeling is how to formulate the radial constraints. There are a considerable number of studies about distribution network reconfiguration, the explanation of which is given here to make this paper self-contained. The reconfiguration of a distribution network is essentially a minimum spanning tree problem with an objective of loss-reduction. Generally, the single substation node scenario is considered in the reconfiguration problem. Thus, to ensure the radial topology after the completion of reconfiguration, the following concise constraint could be adopted to the model: Equation (29) can be interpreted in the following logic. The topology graph consists of vertexes and edges. The sufficient and necessary condition of radial graph is that the number of vertexes equals the number of vertexes minus one. The number one on the right-hand side of (29) indicates the number of substations (i.e., root node). Figure 6 shows a possible reconfiguration process in which the constraint (29) always holds. The operation nature restricts that the topology must be radial, i.e., it must be a single spanning tree. However, things change in the failure restoration case. Because there are potentially two kinds of nodes after the change of topology completed, i.e., energized nodes and non-energized nodes which can form multiple not-connected subgraphs. In this case, we call this subgraph a spanning forest consisting of multiple spanning trees. We will give a simple example to demonstrate this case in Figure 7. The branch S1 is assumed to be faulty, the downstream area, including nodes 2,4,5 is de-energized. After the service restoration is executed, the tie-branch S7 is closed, branch S4 is open. Thus, the original graph is divided into two parts, i.e., two spanning trees. In this case, the value of left-hand side of (29) is 5, but the right-hand side equals to 6. It is obvious that (29) is valid to make topology radial only when the second term of right-hand side equals to the numbers of formed trees, in Figure 7 this number is 2. However, we cannot be aware of how many trees will be formed in advance. Therefore, the radial constraint (29) in reconfiguration cannot be applied to the failure restoration straightly.  In this paper, we will introduce new radial constraints for the failure restoration. Since we do not know the exact numbers of the possible trees that would be formed, the idea is to connect all the remaining trees except for the restored part. Thus, the fictitious node is introduced to be the root node of the de-energized area, the fictitious branches connecting the real load nodes and fictitious node are generated correspondingly. In this way, there will always be two trees after the failure restoration executed. The following mathematical constraints are introduced: Ω l f is the set of fictitious branch. H ij is the fictitious flow of branch (i,j).
are fictitious generation and demand in node i, respectively. Each de-energized node (i.e., those y i equals 0) requires a little fictitious load demand, and only the fictitious node N f supplies the fictitious generation. Thus, Equation (30) ensures that each de-energized node connects the fictitious root node N f . After applying these constraints, the graph after the failure restoration is a 2-tree, the original minimum spanning forest problem has been transformed into a 2-tree spanning problem. Back to Equation (29), since we already know there must be two trees after the failure restoration, the corresponding radial constraint could be modified into following form: However, there is another situation we should take into consideration which the existing studies have not mentioned before. That is, the de-energized node can potentially be contained in the energized tree. The reason is that the load priority w i has been included in the objective function. In general, if we do not account for w i , the failure restoration would restore the loads which have a shorter path to the substation node. Because the shorter path has less possibility to violate the voltage and capacity constraints. Figure 8b shows this common situation in which the nodes 2,4 could not get restored while node five get energized again. However, if taking the load priority into consideration, let's assume that node 2 has a higher priority than node 5 (i.e., w 2 > w 5 ), the restored topology could be the Figure 8c in which the de-energized node 5 is included into the energized area, the corresponding Equation (33) is not valid, the left-hand side of which is 7 while the right-hand side is 6. Furthermore, we will introduce another indicator variable to exclude the above invalidity. Intuitively, if we know the number of nodes which are de-energized nodes embodied in the energized area, then Equation (33) can be modified to valid. These nodes have an obvious characteristic that they have at least one child node which has a status of being restored. In Figure 8c, the child node of 5 is node 2 which is restored. Therefore, we can exclude these nodes from the connection with the fictitious network by force the corresponding fictitious branch to be open (branch N f -5 in Figure 8c). In this way, (33) can be valid. The general mathematical formulation of above analysis can be represented by the following logical constraints: If a node has at least one restored node, i.e., jk∈Ω elc l Then the fictitious branch which connects that node should be open, i.e., x ij = 0. If a node does not have restored node, i.e., Then the fictitious branch which connects that node could be open or closed, i.e., x ij is binary. We then combine the above logical constraints to the following compact linear form: The parameter M in (36) has a similar idea as the Big-M method, but in this case the M represents the maximum numbers of children nodes of one specific node. Now, we have finished the processing of radial constraints in failure restoration. Then, we can obtain the deterministic failure restoration model which consists of objective function (1), operation constraints (2)-(28), fictitious network constraints (30)-(32), and radial constraints (33)(36).

Robust Model Formulation
Wind turbines and photovoltaic systems are uncontrollable distributed generations, and their outputs have substantial uncertainty. Meanwhile, the load demands change all the time. Despite the fact there are many studies on methods for forecasting the output of distributed generations and load demands, errors in forecasting are unavoidable. The focus of this paper is not on the method of forecasting, but is based on the assumption that the forecasting range of distributed generations and load demands have been obtained. Based on this, the uncertainty set in the form of a box is used to describe the forecasted error: where: P i is the forecasted error. g i is the uncertainty coefficient. (39) represents the set of uncertainty in robust optimization theory. The value of τ is called the uncertainty budget. The conservativeness of the robust model can be controlled by adjusting the value of τ. When τ = 0, the uncertainty set is empty, and the robust model degenerates into the traditional deterministic model. When τ = |Ω e |, the uncertainty set is the largest corresponding to the worst scenario which means the maximum deviation from the forecasted value. By choosing the appropriate uncertainty budget τ, the extreme cases can be avoided.
With the uncertainty has been modeled, the robust failure restoration model can be formulated as the following three-layer mixed-integer linear programming form: Cy (41) indicates the relationship between topology and energy flow, Equation (42) is the topology constraints. The electric and natural gas flow constraints and remaining operation constraints are represented by inequalities (43). A, C, D, b, f are the constant matrix or vector related to the corresponding constraint.
The first-layer determines the sequential topology and switch actions. For the uncertainty, the min problem of the second-layer looks for the worst operating scenario in which the minimal restored load demands would be achieved. The third-layer aims to find the optimal dispatching strategy to recover load demands as much as possible.
The three-layer robust optimization problem is difficult to solve directly. For such a problem, the C&CG algorithm [42] and the benders decomposition algorithm [49] are often adopted. Both algorithms use an iterative framework which decomposes the original problem into two sub-problems.
The difference is that the original tangent planes and generated variables are added to the main problem in the C&CG algorithm. Thus, the number of iterations is significantly reduced compared to the benders decomposition method. Therefore, this paper uses C&CG algorithm to solve (41)-(43).

Case Study
In this section, an IEGS composed of a 13-node electric network and a 6-node natural gas network is utilized as a case study to verify the effectiveness of the proposed method. The simulation environment is under MATLAB. The optimization problem is modeled by YALMIP [50] and solved by CPLEX [51].

Simulation Parameters
The topology of the IEGS is shown in Figure 9. In order to distinguish the different nodes of two networks, the letters G and E are used to represent natural gas and electric nodes, respectively. Moreover, the brackets are used to represent lines and pipelines. The dotted lines represent the tie-lines which are open under normal condition. Two networks are coupled by one GFT at node E4 and one P2G at node E9. The original three-phase unbalanced IEEE 13-node distribution network [52] has been modified into a balanced network. The line parameters of IEEE 13-node distribution network are shown in Table 1, and the base electric loads are shown in Table 2, which will fluctuate by percentage over time, as shown in Figure 10. The weight coefficients are randomly generated from 1 to 5 with a uniform distribution which could be modified by the operators on the basis of real-life condition. The pipeline parameters of the natural gas network are shown in Table 3, the loads data of the natural gas network is shown in Table 4. There are two compressors in the pipeline (G3, G5) and (G1, G2), Table 5 shows the parameters. The specific parameters of these two devices are shown in Table 6. Two energy storage devices are embedded in the system, Table 7 shows the specific parameters.

Failure Restoration Process Analysis
Firstly, the objective function value (i.e., total restored load with a weighted coefficient, including electric and natural gas load) are compared in the three cases. As shown in Table 8, the objective value of three cases increase in turn, and the computation time is about 1 second. For Case 1, the excess natural gas in the GSS cannot be converted into electricity because of the lack of existence of an electricity-gas coupling device, so the downstream of the failed pipeline gets no energy supply. Thus, the total restored load is relatively low, just only 55.51% of Case 3. For Case 2, only the GFT coupling is considered, and the excess natural gas can be transformed into electricity, resulting in a higher restored electric load compared to Case 1. However, the downstream of the failure pipeline in the natural gas network cannot be restored either, so the natural gas parts in Case 1 and 2 are the same. The objective function value of Case 2 is about 86.45% of Case 3. However, due to the failure to restore the natural gas load, the electric load is higher than Case 3. In Case 3, the extra electricity generated by the photovoltaic system can be converted into natural gas by the P2G, thus further restoring the natural gas load. Among the three cases, Case 3 has the maximum objective function value 39,874.69, which means that the maximum electric and gas loads are restored, the effectiveness of the proposed two-way coupling failure restoration is demonstrated. Secondly, the restoration process in each case would be analyzed in detail. The restoration actions of lines/pipelines and load in three cases are shown in Table 9. In each step of the restoration process, all constraints in (1) need to be satisfied. In Case 3, compared with Case 1 and Case 2, six extra loads (electric load: E4, E9; gas load: G3, G4, G5, G6) and four loads (gas load: G3, G4, G5, G6) are restored, respectively. - In Case 1, because no extra natural gas is transformed into electricity, the load in the electric network is entirely supplied by the BSS and DG. Thus, the DG located at node E6 and E10 are restored first, the downstream load E8 and E12 are restored through closing the tie-lines (E6, E11), (E12, E8) as shown in Figure 11. To emphasis the change of topology, DGs and BSS, ESS are not drawn in the figure.
In Case 2, due to the existence of GFT at node E4, the electricity it generates which assists the restoration of the electric network. The downstream nodes E9 and E8 are restored through closing the tie-line (E4, E9). Meanwhile, because of the local electricity injection at node E4, the output of BSS at node E1 is reduced, and the extra load E7 is restored. Besides, with the support of DG at node E6 and E10, all electric loads are restored. The restored topology is shown in Figure 12.  In the restoration process of Case 3, the natural gas loads can be restored because the electricity can be transformed into natural gas through the P2G at node E9, as shown in Figure 13. The restoration path includes DG at nodes E6 and E10 so that the clean energy generated by them could be transferred to natural gas through P2G. However, due to the long restoration path, the excessive load will lead to the nodal voltage violating the constraints, so the electric load at E5 and E7 would not be restored. Although the amount of restored electric load is less than Case 2, the total restored load in Case 3 is higher than Case 2 because of the extra restoration of the natural gas load.
Through the co-operation of line switch actions, BSS, GSS, P2G and GFT, the restoration process of Case 3 achieves the multi-period failure restoration of the IEGS the maximum interrupted loads are restored, and also the clean energy's penetration is promoted. This process embodies the complementary advantages of IEGS which increase the reliability of energy supply.
The BSS and coupling devices play a key role in the process of failure restoration. Figure 14 shows the active output of BSS, GFT and the active power consumption of P2G. In general, with the temporal change of the network topology and load, the output of coupling devices changes in real-time in order to make the operation constraints hold. In the fourth and seventh step of the restoration process, DGs are connected to the grid to generate electricity, resulting in the reduction of the output of BSS. The third and eighth step corresponds to the restart-up of GFT and P2G, respectively.

Comparison of Different Radial Topology Constraints
In order to verify the effectiveness of the proposed radial topology constraints, we design another test to compare the existed modeling method with our model. Assuming there is a fault happens at the line (E2, E7), the topology of the 13-node distribution network has been represented by the graph consisting of vertexes and edges, as shown in Figure 15. For highlighting the superiority of the proposed radial topology constraints, we will focus on the electric distribution network, and the natural gas network is ignored in this test. Once the fault detection and isolation accomplished, the de-energized area can be identified which is set of load buses including E7, E8, E9, E10, E11, E12, E13 in Figure 15. Then, the two restoration modeling methods are applied to recover the interrupted loads. Figure 16 shows the restoration results using the traditional radial topology constraints [52] modeling method, the load buses E8, E9, E10, E11, E13 get restored. The E7 and E12 do not get restored because the radial topology constraints enforce that there must be two trees after the restoration. Then, we utilize the radial topology constraints proposed in this paper to make the simulation, the restoration results have a better performance on the amount of restored load demands which is 24.9% more than modeling method in [53]. As shown in Figure 17, the restored loads include E7, E8, E10, E12, E13, the unrestored loads are E9 and E11. The reason for getting the different restoration results is that the proposed radial topology constraints allow the formation of flexible trees after the restoration.
In this case, only one tree is formed, the E11 and E9 are shed other than E7 and E12 because of higher priority coefficient of the loads which makes the better restoration results.

Comparison of Deterministic and Robust Model
Then, in order to prove the effectiveness of the proposed method under uncertainty, Case 4 is designed for simulation. In Case 4, the worst scenario of uncertainty is considered. Meanwhile, the deterministic model of failure restoration is utilized. As shown in Table 10, the objective function value of Case 4 is 92.03% of Case 3 which means that the deterministic model affected by the uncertainty. There is a need to shed more load under the worst scenario which cause the less objective function value, i.e., the restored energy demands, including the electrical and gas loads.  Table 11 illustrates the restoration results under different DGs penetration. We can observe that with the increase of the DGs penetration, the total restored loads get reduced which means the more loads are shed. The reason behind this is that the proposed model is based on robust optimization, the uncertainty has a great impact on the restoration results. The robust model makes all the restoration solutions feasible under the worst scenario of the uncertainty, thus it needs to shed more loads to handle the potential variation of the DGs' output. Meanwhile, we further carry out a test to compare the results obtained by the deterministic model and the proposed robust model under normal uncertain scenarios. We randomly generate 2000 fluctuation scenarios of the uncertain output of DGs using the Monte Carlo sampling method. The restoration results are presented in Table 12. The success rate represents the portion of scenarios in which the restoration results do not violate the nodal voltage or line capacity constraints. It can be seen that the proposed robust model can get a valid restoration solution in all scenarios. Moreover, the average restored power of robust model is greater than the deterministic model. Thus, the proposed robust model has a superiority than the deterministic model under the condition of uncertainty.

Conclusions
Most of the existing research on failure-related aspects of integrated electricity-gas systems focuses on failure prevention and reliability assessment. From the post-failure perspective, this paper proposes a multi-period failure restoration method considering two-way coupling for IEGS. The problem is formulated as a mixed-integer linear programming model, utilizing linearization technology and linear Distflow equations. On the one hand, the proposed failure restoration strategy considering the bidirectional coupling of the electric and natural gas networks has the best result compared to other kinds of coupling; on the other hand, the proposed failure restoration strategy could coordinate the sequential actions of switchgear, BSS, GSS, GFT, and P2G, guaranteeing the safe operation during the whole failure restoration process. The proposed radial topology constraints modeling method can get a better restoration result than the traditional modeling method. Besides, the case study shows that the robust model can handle the uncertainty problem, leading to a more reliable restoration result than the deterministic model. Based on this paper, follow-up works can focus on dynamic natural gas network modeling, and heating network coupling.