1 Introduction

In recent years, renewable energy generation gains rising attention due to the lack of traditional resources. However, the electricity quality and the reliability of power grids are significantly affected by the intermittency and instability of renewable energy resources such as wind and solar. By contrast, natural gas is a more stable and reliable sort of resource which can provide continuous energy for both gas and electricity loads by gas-fired generators [1].

Due to the clean, efficient and high-quality characteristics of natural gas, it has been widely used as the main energy resource in some areas and the coordinated operation of natural gas and electricity system has been researched in many previous studies [2]. A basic model for integrating natural gas and electricity networks is presented in [3], which shows the fundamentals of natural gas network and describes the constraints for the energy transmission between electricity and gas systems. In [4], a steady state power flow model is presented for solving the combined optimization operating problem of different energy facilities based on the new concept of energy hubs. A decomposition method is applied to solve the security-based model proposed in [5] for the solution of SCUC problem considering natural gas transmission system. In [6], an mixed integer linear program (MILP) method is presented to formulate the optimal power flow in multi-carrier energy systems. The non-convex constraints of natural gas transmission system are linearized, so that the problem is reformulated as an MILP problem which can be solved by traditional optimization methods.

However, the operating stability and reliability are still affected by increasing uncertainties in energy systems. Disruptions in an energy system are sometimes inevitable, uncontrollable and unpredictable [7]. As service industry, the energy system must guarantee the continuity of energy supply for customers. Therefore, improving resilience of energy systems is of vital importance. Resilience is defined as the ability to provide and maintain an acceptable level of service in the face of faults and challenges to normal operation [8]. The contingencies and challenges for services range from natural disaster to terrorist attacks. As a role of defender, the system itself will take preventive measures for attacks before disruptions occur and respond to the damage after attacks. The resiliency analysis in electricity system has been studied by several researchers. In [9], a resiliency-oriented microgrid optimal scheduling model is proposed for minimizing the load curtailments when the service of main grid is interrupted. A non-cooperative game-theoretic framework is presented in [10] to study the strategic behavior of microgrids. The framework incorporates economic factors and stability and efficiency of microgrids, which is solved by fully distributed phasor measurement unit (PMU)-enabled algorithm to ensure the resiliency of the proposed method. In [11], a resilient distribution network planning problem is presented and formulated as a two-stage robust optimization model to minimizing the system damage by coordinating the hardening and distributed resource allocation. In [12], a distribution system operating method by microgrid formation after natural disaster is proposed to restore important loads from power outage. A tri-level optimization model for electric power system defense is presented in [13] which can identify critical elements in power grid to defend against unpredictable attacks. A column-and-constraint generation algorithm is applied to solve a two-stage robust optimization problem in [14]. In [7], a risk assessment model is proposed to determine potential vulnerabilities of power system and provide feasible plans for enhanced protections according to the budgets for power grid construction.

There are few work which has been done to analyze resilience of integrated energy system which includes electricity, gas and other forms of energy. As the operational feature of natural gas system is different from that in power grids and the operating conditions have impacts on electricity networks, the resiliency analysis is an essential topic to be researched. In [15], a methodology is proposed to locate the most vulnerable components to make sure the resilient operation of multiple energy carrier microgrids when terrorists attack the network. The model is formulated as a bi-level optimization problem to solve the optimal operation for multi-energy microgrids in consideration of security and resiliency. In [16], a methodology is proposed to identify and protect vulnerable components of integrated electric and gas infrastructures. The resilience is guaranteed by solving a tri-level optimization problem. A mixed integer linear programming and nested column-and-constraint generation algorithm is applied to solve the proposed model. A novel mixed integer linear programming for security-constrained power and gas flow is presented in [17]. The proposed model allows the integrated system operates in both normal and contingency conditions with the least violations. In [18], three models are proposed for identifying optimal energy flow solvability to ensure secure operating conditions with corrective controls.

In this paper, a coordinated operation model of energy system which integrates electricity and natural gas infrastructures is formulated. The model includes constraints for both electricity and natural gas transmission networks. Then, a three-stage robust optimization algorithm is presented to solve the defender-attacker-defender problem of the integrated energy system when contingencies occur. The constraints and operational features of natural gas system and the coupling of electricity and natural gas networks make it more complicated to settle the plans for both attackers and defenders. In the first stage, as defender, the integrated energy system must make plans for network enhancement to minimize the damage caused by unpredictable attacks. In the second stage, attackers will attack vulnerable components of electricity and natural gas networks to cause maximum damage to the entire energy system. In the last stage, defender responds to the results of disruptions, which is formulated as optimal energy flow of integrated electricity and natural gas system. As some of the constraints of electricity and natural gas transmission networks are nonlinear, Taylor series expansion algorithm is applied to realize the linearization. Then, the proposed three-stage robust optimization problem is reformulated into a two-stage optimization problem by the application of duality theory. The two-stage problem can be solved as an MILP by decomposition algorithms.

The main contributions of our paper are:

  1. 1)

    A nested Benders decomposition algorithm is applied to solve the proposed defender-attacker-defender problem. It is more effective to be applied for large-scale problems.

  2. 2)

    Resiliency of integrated energy systems is analyzed in the presence of energy storage system. By adding proper electricity and gas storage facilities, reliability of integrated energy system is re-analyzed.

Compared with [15], we incorporate the hardening plans before attacks, which the resilience of the whole system can be better improved. That is to say, a defender-attacker defender model is presented in comparison with the attacker-defender model in [15]. Compared with [16], we apply a nested Benders decomposition algorithm to solve the nested two-stage proposed model. Moreover, Taylor series expansion is applied to linearize the quadratic polynomials of cost functions and the gas flow square for gas pipelines and compressors. In addition, energy storage systems are also considered in the end of this paper to analyze the resilience of integrated energy systems.

The remainder of this paper is organized as follows. Section 2 describes the mathematical formulations for the optimal energy flow and the three-stage robust optimization model of the integrated electricity and natural gas system. Section 3 provides the solution methodologies to linearize and decompose the presented three-stage problem. Section 4 presents and discusses the numerical results of the proposed method and analyze the impact of the presence of energy storage devices. Finally, the conclusion is given in Section 5.

2 Mathematical formulation

In this section, a power flow model for coordinated electricity and gas network is presented. Based on this integrated system, a robust optimization model is applied for resiliency analysis. Any contingency occurred in either electricity system or gas network may cause coordination problems or safety and stability problems of both systems. Models are presented to figure out the optimal defense and operation plans against contingencies.

2.1 Integrated electricity and gas network optimal flow model

Figure 1 shows a general network of natural gas transmission system. Gas suppliers, compressors, storage systems, pipelines and gas loads are essential components in gas networks. Natural gas is transferred from producers to end customers through pipelines and compressors. Unlike electricity, natural gas can be injected into certain gas storage facilities during off-peak periods and used during high-demand periods [19]. Therefore, steady flow over pipelines can be easily maintained. Gas in pipelines has much lower speed than electricity and can be stored in pipelines in a short period of time due to its compressibility [20]. So the dynamic process is a key feature of gas network, and usage of steady-state gas flow models could result in sub-optimal results or even insecure operation decision. For simplicity, in this paper, we only analyze steady state behavior for integrated electricity and gas system. Therefore, for modeling natural gas system, we assume that the natural gas system operates in steady states and the line pack is ignored [16].

Fig. 1
figure 1

General network of natural gas transmission system

Full size image

For pipelines, there are two major variables which are gas pressure \(p_{m}\) at each node \(m\) and gas flow \(f_{mn}\) between two end nodes \(m\) and \(n\). Similar to the voltage in power system, gas pressure insures that gas can transfer from one gas node to another. In general, gas can only be delivered from higher pressure nodes to lower ones. For each node, there are gas flow injection and gas loads. Mathematically, the gas flow balance equation is expressed as:

$$s_{m} + \sum\limits_{(m,n) \in \varPhi } {f_{mn} } + f_{sm} = d_{m} + e_{m} (P_{Gm} )$$
(1)

where \(s_{m}\) is the gas supply from producers at node m; \(f_{sm}\) is the gas storage at node m; \(d_{m}\) is the non-electrical gas demand at node m; \(e_{m}\) is the gas demand for gas-fired generator at node m; \(P_{Gm}\) is the power output of the generator node m.

For each pipeline, the amount of gas flow is associated with end nodes pressures and properties of each pipeline [5], which is presented in (2).

$$f_{mn} = {\text{sgn}}(p_{m} - p_{n} ) \cdot C_{mn} \sqrt {\left| {p_{m}^{2} - p_{n}^{2} } \right|}$$
(2)

where \(\text{sgn} ()\) is the sign function; \(C_{mn}\) is a constant which depends on properties of pipelines such as length, diameter, temperature, friction and gas composition.

For pipelines which have compressors, natural gas flows from the lower pressure node to higher ones because of the existence of compressors. Inevitable friction between gas and pipelines will result in pressure loss which can be compensated by compressors [5]. The gas flows through compressors can be calculated in (3).

$$f_{mn}^{cp} = {\text{sgn}}(p_{m} - p_{n} )\frac{{H_{mn} }}{{k_{mn2} - k_{mn1} \left[ {\frac{{{{\rm max}} (p_{m} ,p_{n} )}}{{{{\rm min}} (p_{m} ,p_{n} )}}} \right]^{\alpha } }}$$
(3)

where \(H_{mn}\) is the power input of compressor between nodes m and n; \(k_{mn1} ,\;k_{mn2} ,\;\alpha\) are empirical parameters that depend on properties of compressors. Gas consumed by compressors is treated as transmission “loss” of natural gas network, which is a quadratic function of \(H_{mn}\):

$$f_{{H_{mn} }} (H_{mn} ) = a_{mn} H_{mn}^{2} + b_{mn} H_{mn} + c_{mn}$$
(4)

It is a simplified form of expression which can be found in [5].

The optimal flow for integrated electricity and gas system takes the minimum total cost of operation as objective function:

$${\text{min }}\sum\limits_{{i \in N_{G} }} {(a_{i} P_{Gi}^{2} + b_{i} P_{Gi} + c_{i} )} + \sum\limits_{{m \in N_{S} }} {\varLambda_{m} s_{m} }$$
(5)

where \(N_{G}\) is the set of non-gas-fired generator nodes; \(N_{S}\) is the set of natural gas nodes; \(\varLambda_{m}\) is the price of gas at node m.

The constraints of electricity system are shown in (6)–(12).

$$P_{Gi} - \text{Re} \left( {\sum {S_{ij} (V,\theta )} } \right) = P_{Li}$$
(6)
$$Q_{Gi} - \text{Im} \left( {\sum {S_{ij} (V,\theta )} } \right) = Q_{Li}$$
(7)
$$P_{Gi,{{\rm min}} } \le P_{Gi} \le P_{Gi,{{\rm max}} }$$
(8)
$$Q_{Gi,{{\rm min}} } \le Q_{Gi} \le Q_{Gi,{{\rm max}} }$$
(9)
$$V_{i,{{\rm min}} } \le V_{i} \le V_{i,{{\rm max}} }$$
(10)
$$- \theta_{i,{{\rm max}} } \le \theta_{i} \le \theta_{i,{{\rm max}} }$$
(11)
$$\left| {P_{ij} (V,\theta )} \right| \le P_{ij,{{\rm max}} }$$
(12)

Equations (6) and (7) are the power balance formulation; (8)–(12) represent the constraints of node voltage, generator output and lines transmission limits.

Apart from (1)–(5), the other constraints for natural gas system are presented in (13)–(15). λ is the dual variables for the constraints.

$$s_{m,{{\rm min}} } \le s_{m} \le s_{m,{{\rm max}} } \quad \lambda_{m}^{ 2} ,\lambda_{m}^{ 3}$$
(13)
$$p_{m,{{\rm min}} } \le p_{m} \le p_{m,{{\rm max}} } \quad \lambda_{m}^{\text{4}} \text{,}\lambda_{m}^{\text{5}}$$
(14)
$$H_{mn,{{\rm min}} } \le H_{mn} \le H_{mn,{{\rm max}} } \quad \lambda_{m}^{6} ,\lambda_{m}^{7}$$
(15)

Equations (13) and (14) are respectively the gas supply and gas pressure limits.

2.2 Robust optimal operation model for resilient energy systems

As the attacks on integrated energy system are random and unpredictable, the complex nature of this problem makes robust optimization be the most suitable method that can take account of inherent randomness and uncertainties [11]. This robust optimization problem aims to find the optimal enhancement and dispatch plans according to the worst case caused by attackers.

As described in the first section, resiliency analysis for integrated energy system can be formulated in three stages, which is also known as defender-attacker-defender game model. The first stage is the hardening of transmission lines by defender to minimize the damage caused by contingencies. In the second stage, attackers disrupt the energy system to produce a most damaging attack. In the third stage, defender will respond to the damage caused by attacks to optimize the operation of integrated energy system. Therefore, the three-stage problem can be formulated as a min-max-min problem which is shown as follows:

$$\mathop {\text{min}}\limits_{h \in H} \,\mathop {\text{max}}\limits_{a \in A} \,\mathop {\text{min}}\limits_{r \in R(H, A)} (f_{C} (r) + f_{LC} (r))$$
(16)

where h and a are respectively the hardening and attacking scenarios; H and A are respectively the feasible set for line hardening plans and the uncertainty set for attacking plans; r represents the optimal electricity and gas flow variables; \(R(H,A)\) denotes the energy system response set based on H and A; \(f_{LC} (r)\) is the total economic loss due to load curtailment; \(f_{C} (r)\) is the costs for optimal operation of energy system.

In this section, we analyze the model from inner minimum problem to outer minimum problem.

2.2.1 Defender response model

In Section 2.1, we have discussed the optimal operation planning for electricity and gas network in normal situation. After an attack, the topology of electricity and gas network may change, which will result in inevitable load curtailment. In defender response model, reactive power is ignored when calculating power flow in electricity network. Therefore, constraints for power system are modified in (17)–(21).

$$P_{ij} x_{ij} = \theta_{i} - \theta_{j}$$
(17)
$$P_{Gi} + \sum\limits_{j \in J} {P_{ij} } + LC_{i}^{elec} = P_{Li}$$
(18)
$$- P_{ij,{{\rm max}} } \le P_{ij} \le P_{ij,{{\rm max}} }$$
(19)
$$P_{Gi,{{\rm min}} } \le P_{Gi} \le P_{Gi,{{\rm max}} } \quad \lambda_{i}^{8} ,\lambda_{i}^{9}$$
(20)
$$- \theta_{i,{{\rm max}} } \le \theta_{i} \le \theta_{i,{{\rm max}} } \quad \lambda_{i}^{10} ,\lambda_{i}^{11}$$
(21)

where \(LC_{i}^{elec}\) is the load curtailments for electricity network. The constraint for natural gas network in (1) is modified as:

$$\begin{aligned} s_{m} & = \sum\limits_{{(m,n) \in N_{g} }} {f_{{H_{mn} }} (H_{mn} )} + d_{m} + e_{m} (P_{Gm} ) \\ & \quad - \sum\limits_{{(m,n) \in N_{CP} }} {f_{mn} } - LC_{m}^{gas} \\ \end{aligned}$$
(22)

where \(LC_{m}^{gas}\) is the load curtailment for natural gas network. The defender response model (DRM) can be simply described in (23).

$$\mathop {\text{min}}\limits_{{\left( {P_{G} ,\theta ,p,s} \right) \in R}} (f_{C} (r) + f_{LC} (r))$$
(23)
$$f_{C} (r) = f_{c}^{elec} (P_{G} ,\theta ) + f_{c}^{gas} (p,s)$$
(24)
$$f_{LC} (r) = \sum\limits_{i \in I} {\tau_{i}^{elec} \cdot LC_{i}^{elec} } + \sum\limits_{m \in M} {\tau_{m}^{gas} \cdot LC_{m}^{gas} }$$
(25)

subject to (2)–(4) and (16)–(21).

where I and M represent the sets of load curtailment nodes for electricity and gas network respectively; \(\tau_{i}^{elec} \;{\text{and}}\;\tau_{m}^{gas}\) denote the coefficients of economic loss for electricity and gas network respectively. Equation (24) is the simplified objective function of (5).

2.2.2 Attacker interdiction model

There are various of attacks that may result in system blackout, which can be generally categorized into two kinds: terrorist attack and natural disaster. However, attackers will destroy network elements, such as generators, compressors or transmission lines, regardless of the type of attacks. As energy transmission failure occurs when network elements are attacked, we assume that attackers only attack transmission lines for simplicity. We introduce \(a_{t,ij}^{elec}\) and \(a_{t,mn}^{gas}\) to represent the attacking states for electricity and gas transmission lines respectively. Then we set:

$$a_{t,ij}^{elec} = \left\{ {\begin{array}{*{20}l} 1 \hfill & {\text{if electricity transmission line }}i{\hbox{-}}j{\text{ is attacked}} \hfill \\ 0 \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$
(26)
$$a_{t,mn}^{gas} = \left\{ {\begin{array}{*{20}l} 1 \hfill & {\text{if gas transmission line }}m{\hbox{-}}n{\text{ is attacked}} \hfill \\ 0 \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$
(27)

The lines constraints (2), (3), (17)–(19) and (22) must be modified due to the introduced variables \(a_{t,ij}^{elec}\) and \(a_{t,mn}^{gas}\).

$$\begin{aligned} s_{m} & = \sum\limits_{{(m,n) \in N_{CP} }} {a_{t,mn}^{gas} f_{{H_{mn} }} (H_{mn} )} + d_{m} + e_{m} (P_{Gm} ) \\ & \quad - \sum\limits_{{(m,n) \in N_{g} }} {a_{t,mn}^{gas} f_{mn} } - LC_{m}^{gas} \\ \end{aligned}$$
(28)
$$f_{mn} = a_{t,mn}^{gas} \cdot {\text{sgn}}(p_{m} - p_{n} ) \cdot C_{mn} \sqrt {\left| {p_{m}^{2} - p_{n}^{2} } \right|}$$
(29)
$$f_{mn}^{cp} = a_{t,mn}^{gas} \frac{{{\text{sgn}}(p_{m} - p_{n} )H_{mn} }}{{k_{mn2} - k_{mn1} \left[ {\frac{{{{\rm max}} (p_{m} ,p_{n} )}}{{{{\rm min}} (p_{m} ,p_{n} )}}} \right]^{\alpha } }}$$
(30)
$$P_{ij} x_{ij} = a_{t,mn}^{gas} (\theta_{i} - \theta_{j} )\quad \lambda_{i}^{15}$$
(31)
$$P_{Gi} + \sum\limits_{j \in J} {P_{ij} } + LC_{i}^{elec} = P_{Li} \quad \lambda_{i}^{16}$$
(32)
$$\left| {P_{ij} } \right| \le a_{t,ij}^{elec} P_{ij,{{\rm max}} } \quad \lambda_{i}^{17} ,\lambda_{i}^{18}$$
(33)

As some of generators in electricity system are customers of natural gas system, attacks on gas transmission lines may cause blackouts of generators, compressors and gas loads, which will result in more damage to entire energy system. Therefore, attackers’ plan for integrated energy system will be quite different from the plan in power system only. Attackers will find the optimal attacking plan which maximizes the economic loss caused by electricity and gas load curtailment. The attacker interdiction problem (AIP) is formulated in (34)–(36).

$$\mathop {\text{max}}\limits_{(I,M) \in A} (f_{C} (r) + f_{LC} (r))$$
(34)

where

$$A = \left\{ {\left. {(LC,I,M)} \right|\mathop {\text{min}}\limits_{{\left( {P_{G} ,\theta ,p,s} \right) \in R}} (f_{C} (r) + f_{LC} (r))} \right\}$$
(35)

subject to (4), (13)–(15), (20), (21), (28)–(33), and

$$a_{t,ij}^{elec} ,\, a_{t,mn}^{gas} \in \left\{ {0,1} \right\}$$
(36)

It is obvious that this model is a two-stage optimization problem.

2.2.3 Defender reinforcement model

Lines reinforcement is an effective measure that defenders use to harden networks and improve resilience of systems. We assume that hardened lines cannot be attacked. Similar to attacker interdiction model discussed in Section 2.2.2, we also introduce \(h_{l,ij}^{elec}\) and \(h_{l,mn}^{gas}\) to represent the hardening states for electricity and gas transmission lines respectively. We set:

$$h_{l,ij}^{elec} = \left\{ {\begin{array}{*{20}l} 1 \hfill & {\text{if electricity transmission line }}i{\hbox{-}}j{\text{ is hardened}} \hfill \\ 0 \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$
(37)
$$h_{l,mn}^{gas} = \left\{ {\begin{array}{*{20}l} 1 \hfill & {\text{if gas transmission line }}m{\hbox{-}}n{\text{ is hardened}} \hfill \\ 0 \hfill & {\text{otherwise}} \hfill \\ \end{array} } \right.$$
(38)

In consideration of economic budget, the number of hardened lines is limited, which is given by (39).

$$\sum\limits_{{i,j \in N_{G} }} {\varsigma_{ij}^{elec} h_{l,ij}^{elec} } { + }\sum\limits_{{m,n \in N_{S} }} {\varsigma_{mn}^{gas} h_{l,mn}^{gas} } \le B_{h}$$
(39)

where \(\varsigma_{ij}^{elec}\) and \(\varsigma_{mn}^{gas}\) are unit cost for lines reinforcement; \(B_{h}\) is the budget for reinforcement.

Defenders will find the optimal reinforcement plan to minimize the economic loss caused by load curtailment. The defender reinforcement problem (DRP) is formulated as follows:

$$\mathop {\text{min}}\limits_{{h = (h_{l,ij}^{elec} ,h_{l,mn}^{gas} ) \in H}} AIP(h)$$
(40)

subject to (30)

$$h_{l,ij}^{elec},\, h_{l,mn}^{gas} \in \left\{ {0,1} \right\}$$
(41)

where AIP is presented in Section 2.2.2. It is obvious that this model is a nested two-stage optimization problem.

3 Solution methodology

The formulated model in Section 2.2 is a nested two-stage robust optimization problem which is shown in Section 2.2.3. We apply a nested Benders decomposition algorithm to solve this problem.

3.1 Master problem

The master problem of the two-stage model is the minimization of economic loss for load curtailment under lines reinforcement.

$${\text{min}}\quad z$$
(42)

subject to

$$z \ge f_{C} (r) + f_{LC} (r)$$
(43)
$$\sum\limits_{{i,j \in N_{G} }} {\varsigma_{ij}^{elec} h_{l,ij}^{elec} } { + }\sum\limits_{{m,n \in N_{S} }} {\varsigma_{mn}^{gas} h_{l,mn}^{gas} } \le B_{h}$$
(44)
$$h_{l,ij}^{elec} ,\,h_{l,mn}^{gas} \in H$$
(45)

Solving the master problem needs a set of network reinforcement plan H and the worst case under attack which can be derived from the subproblem.

3.2 Subproblem

The subproblem is aiming to find the worst case under attack with given reinforcement plan, which is the maximization of economic loss for load curtailment. For a given attacking plan, DRM will formulate an optimal operation plan for integrated system. The subproblem is actually a max-min optimization problem which is formulated in Section 2.2.2 as AIP.

$$\mathop {\text{max}}\limits_{(I,M) \in A} \, \mathop {\text{min}}\limits_{{\left( {P_{G} ,\theta ,p,s} \right) \in R}} (f_{C} (r) + f_{LC} (r))$$
(46)
$${\text{subject}}\,{\text{to}}\quad a_{t,ij}^{elec} ,\,a_{t,mn}^{gas} \in A$$
(47)

and (4), (13)–(15), (20), (21), (28)–(33).

As the objective and constraints in (4), (28)–(30) and (46) are nonlinear, we will find ways to linearize them.

Variables \(a_{mn}\) in (4), \(a_{i}\) in (5) and \(a_{m}\) in (28) are generally quite small, so we take the derivative of the second order terms for linearization. Equations (4), (5), (28) are modified as:

$$f_{{H_{mn} }} (H_{mn} ) = (2a_{mn} { + }b_{mn} )H_{mn} + c_{mn} \quad \lambda_{m}^{1}$$
(48)
$$\text{min} \, \left\{ {\sum\limits_{{i \in N_{G} }} {[(2a_{i} + b_{i} )P_{Gi} + c_{i} ]} + \sum\limits_{{m \in N_{S} }} {\varLambda_{m} s_{m} } } \right\}$$
(49)
$$\begin{aligned} s_{m} & = \sum\limits_{{(m,n) \in N_{CP} }} {a_{t,mn}^{gas} f_{{H_{mn} }} (H_{mn} )} + d_{m} + (2a_{m} + b_{m} )P_{Gm} \\ & \quad + \,c_{m} - \sum\limits_{{(m,n) \in N_{g} }} {a_{t,mn}^{gas} f_{mn} } - LC_{m}^{gas} \quad \lambda_{m}^{12} \\ \end{aligned}$$
(50)

As for the constraints in (29) and (30), we apply Taylor series expansion to linearize them by omitting the higher order terms [6]. Natural gas flow expression (29) is expanded around the neighbor point \((p_{m0} ,p_{n0} )\) as:

$$\begin{aligned} f_{mn} & \approx f_{mn} (p_{m0} ,p_{n0} ) + \frac{{\partial f_{mn} }}{{\partial p_{m} }}(p_{m} - p_{m0} ) + \frac{{\partial f_{mn} }}{{\partial p_{n} }}(p_{n} - p_{n0} ) \\ & = \frac{{a_{t,mn}^{gas} C_{mn} }}{{\sqrt {p_{m0}^{2} - p_{n0}^{2} } }}(p_{m0} p_{m} - p_{n0} p_{n} ) = a_{t,mn}^{gas} (U_{mn} p_{m} - V_{mn} p_{n} )\quad \lambda_{m}^{13} \\ \end{aligned}$$
(51)

where we assume that \(p_{m} > p_{n}\).

Similarly, gas flow through compressors in (30) can be expanded in the neighbor point \((p_{m0} ,\;p_{n0} ,\;H_{mn0} )\) as:

$$\begin{aligned} & f_{mn}^{cp} \approx f_{mn}^{cp} (p_{m0} ,p_{n0} ,H_{mn0} ) + \frac{{\partial f_{mn}^{cp} }}{{\partial p_{m} }}(p_{m} - p_{m0} ) + \frac{{\partial f_{mn}^{cp} }}{{\partial p_{n} }}(p_{n} - p_{n0} ) \\ & \quad + \frac{{\partial f_{mn}^{cp} }}{{\partial H_{mn} }}(H_{mn} - H_{mn0} ) = \frac{{\alpha H_{mn0} k_{1} }}{{\left( {k_{1} \left( {\frac{{p_{n0} }}{{p_{m0} }}} \right)^{\alpha } - k_{2} } \right)^{2} }}\frac{{p_{n0}^{\alpha } }}{{p_{m0}^{\alpha } }}\left( {\frac{{p_{m} }}{{p_{m0} }} - \frac{{p_{n} }}{{p_{n0} }}} \right) \\ & \quad + \frac{{H_{mn} }}{{k_{1} \left( {\frac{{p_{n0} }}{{p_{m0} }}} \right)^{\alpha } - k_{2} }} = a_{t,mn}^{gas} (U_{mn}^{CP} p_{m} - V_{mn}^{CP} p_{n} + W_{mn}^{CP} H_{mn} ) \quad \lambda_{m}^{14} \\ \end{aligned}$$
(52)

As mentioned in [6], the domain of \(f_{mn}\) is divided into (NL)2 grids. NL is the number of segments which is decided by the accuracy requirement. Selection of neighbor points and approximation error analysis are also given in [6], which will not be repeatedly explained here.

To solve this max-min model, the inner minimization problem is transformed into its dual problem. Dual variables for the constraints are given in the end of each formulation, which are \(\lambda^{1} \sim\lambda^{18}\). The dual problem is a maximization program shown in (53)–(67).

$$\begin{aligned} {\text{max }} & \sum\limits_{{m \in N_{g} }} {C_{mn} \lambda_{m}^{1} } + \sum\limits_{{m \in N_{S} }} {(s_{m,{{\rm max}} } \lambda_{m}^{2} - s_{m,{{\rm min}} } \lambda_{m}^{3} )} \\ & + \sum\limits_{{(m,n) \in N_{CP} }} {(H_{mn,{{\rm max}} } \lambda_{m}^{6} - H_{mn,{{\rm min}} } \lambda_{m}^{7} )} \\ & + \sum\limits_{{(i,j) \in N_{el} }} {a_{t,ij}^{elec} P_{ij{{\rm max}} } (\lambda_{i}^{17} - \lambda_{i}^{18} )} + \sum\limits_{{m \in N_{g} }} {(p_{m,{{\rm max}} } \lambda_{m}^{4} - p_{m,{{\rm min}} } \lambda_{m}^{5} )} \\ & + \sum\limits_{{i \in N_{g} }} {(P_{Gi,{{\rm max}} } \lambda_{i}^{8} - P_{Gi,{{\rm min}} } \lambda_{i}^{9} )} + \sum\limits_{{m \in N_{g} }} {d_{m} \lambda_{m}^{12} } + \sum\limits_{{i \in N_{e} }} {P_{Li} \lambda_{i}^{16} } \\ & + \sum\limits_{{i \in N_{e} }} {\theta_{i{{\rm max}} } (\lambda_{i}^{10} - \lambda_{i}^{11} )} \\ \end{aligned}$$
(53)

subject to

$$a_{t,ij}^{elec} ,\,a_{t,mn}^{gas} \in A$$
(54)
$$\lambda_{m}^{2} + \lambda_{m}^{3} + \lambda_{m}^{12} = \varLambda_{m} \quad m \in N_{s}$$
(55)
$$\begin{aligned} & \lambda_{m}^{4} + \lambda_{m}^{5} + a_{t,mn}^{gas} (U_{mn} \lambda_{m}^{13} - V_{mn} \lambda_{n}^{13} )\\ & \quad + a_{t,mn}^{gas} (U_{mn}^{CP} \lambda_{m}^{14} - V_{mn}^{CP} \lambda_{n}^{14} ) = 0\quad (m,n) \in N_{g} \\ \end{aligned}$$
(56)
$$(2a_{mn} + b_{mn} )\lambda_{m}^{1} + \lambda_{m}^{6} + \lambda_{m}^{7} + W_{mn}^{CP} \lambda_{m}^{14} = 0 \, \quad { (}m,n) \in N_{CP}$$
(57)
$$\lambda_{i}^{8} + \lambda_{i}^{9} + (2a_{m} + b_{m} )\lambda_{m}^{12} + \lambda_{i}^{16} = 2a_{i} + b_{i} \quad i \in N_{G} ,\,m \in N_{g}$$
(58)
$$\lambda_{i}^{10} + \lambda_{i}^{11} + a_{t,ij}^{elec} (\lambda_{i}^{15} - \lambda_{j}^{15} ) = 0\quad i,j \in N_{e}$$
(59)
$$\lambda_{i}^{17} + \lambda_{i}^{18} + x_{ij} \lambda_{i}^{15} + \sum {\lambda_{i}^{16} } = 0\quad \, i \in N_{e}$$
(60)
$$\sum {a_{t,mn}^{gas} } \lambda_{m}^{12} + \lambda_{m}^{13} = 0\quad \, (m,n) \in N_{el}$$
(61)
$$- \sum {a_{t,mn}^{gas} } \lambda_{m}^{12} + \lambda_{m}^{14} = 0\quad \, (m,n) \in N_{CP}$$
(62)
$$\lambda_{m}^{12} = \tau_{m}^{gas} \quad \, m \in N_{g}$$
(63)
$$\lambda_{i}^{16} = \tau_{i}^{elec} \quad i \in N_{e}$$
(64)
$$\lambda_{m}^{2} ,\lambda_{m}^{4} ,\lambda_{m}^{6} ,\lambda_{i}^{8} ,\lambda_{i}^{10} ,\lambda_{i}^{17} \le 0$$
(65)
$$\lambda_{m}^{3} ,\lambda_{m}^{5} ,\lambda_{m}^{7} ,\lambda_{i}^{9} ,\lambda_{i}^{11} ,\lambda_{i}^{18} \le 0$$
(66)
$$\lambda_{m}^{1} ,\lambda_{m}^{12} \sim\lambda_{i}^{16} \quad {\text{unlimited}}$$
(67)

The dual maximization problem is linear and can be solved by an MILP solver.

3.3 Solution step for nested two-stage robust optimization algorithm

  1. 1)

    Step 1: Initialization of variables. Set iterations \(K \to 0,\,LB \to - \infty ,\,UB \to \infty .\)

  2. 2)

    Step 2: Solve the subproblem with given reinforcement plan \(H^{ * }\). Get the objective maximum economic loss \(obj_{SP}^{{{\rm max}} }\) and the optimal attacking plan \(A^{ * }\) for the worst case. Set \(UB \to {\text{Min}}(UB,obj_{SP}^{{{\rm max}} } ),\,\xi_{c} \to UB - LB,\,K = K + 1.\)

  3. 3)

    Step 3: Solve the master problem with attacking plan \(A^{ * }\) and optimal system operation variables R derived from Step 2. Get the minimum economic loss \(obj_{MP}^{{{\rm min}} }\). Set \(LB \to {\text{Max}}(LB,obj_{MP}^{{{\rm min}} } ),\,\xi_{c} \to UB - LB\). Update lines reinforcement plan \(H^{ * }\).

  4. 4)

    Step 4: If \(\xi_{c}\) satisfies convergence condition, stop the process. Otherwise, return to Step 2.

4 Case study

To show the performance of our proposed optimization model for integrated electricity and natural gas energy system, we apply two testing systems which are IEEE 30-bus power system with 7-node natural gas system and IEEE 118-bus power system with 14-node natural gas system.

4.1 IEEE 30-bus system

This case is based on the modified IEEE 30-bus power network and 7-node natural gas system.

The modified IEEE-30 bus power system shown in Fig. 2 is consist of 6 thermal units including 3 fossil units and 3 natural gas-fired units, 41 branches, 4 transformers and 20 demand sides. Data for buses, branches and load demands are from [21]. Characteristics for generator locations, costs and limits are taken from [22]. The 7-node natural gas system depicted in Fig. 3 has 2 natural gas suppliers, 5 pipelines, 1 compressor and 5 natural gas loads including 3 demand loads for gas-fired units. Parameters for natural gas network consults [5].

Fig. 2
figure 2

IEEE 30-bus power system

Full size image
Fig. 3
figure 3

7-node natural gas system

Full size image

In normal conditions where no transmission lines are hardened and no elements are damaged, the integrated energy system will operate in a safe and economic manner. Some of the optimal dispatch plan results for the system in a certain time period are given in Tables 1–3. L1, L3 and L5 are variable gas loads for natural gas-fired units G2, G3 and G4, which depends on the outputs and loads of electricity network. In normal operation conditions, no congestions happen and no loads are curtailed.

Table 1 Units outputs for electricity network
Full size table
Table 2 Gas outputs for natural gas network
Full size table
Table 3 Gas loads for natural gas network
Full size table

In resilient operation conditions, lines reinforcement and load curtailment are considered. As both of defenders and attackers have operating budgets, defenders cannot reinforce all transmission lines and attackers cannot damage all elements in the system. Therefore, in this case, we assume that attackers have budgets for damaging up to 5 lines in electricity network and 1 line in natural gas network. Then we analyze the reinforcement plan for defenders of which budgets can support from 0 to 5 lines reinforcement for electricity network and 0 to 1 line reinforcement for natural gas network. The lines reinforcement plans for integrated system under different budgets are given in Table 4. In Table 4, HLelec and HLgas are hardened lines for electricity and natural gas networks respectively.

Table 4 Reinforcement plans for integrated energy system
Full size table

When no line is hardened, attackers will cause damage up to 115.1 MW electricity load curtailments and 4764 kcf natural gas load curtailments. By hardening 5 lines in electricity network and 1 line in gas network, the worst case only causes 57.5 MW electricity load curtailments and 2351 kcf natural gas load curtailments, which is a reduction of nearly 50%. Congestions caused by attacking for 6 different defending plans are shown in Table 5. Bcg is the branch that congestion happens. Sij and Slim are apparent power and limits for transmission lines. It is obvious that more hardened lines will alleviate congestions of transmission lines. Therefore, resilience of integrated energy system is improved greatly.

Table 5 Congestions caused by attacks for 6 reinforcemet plans
Full size table

Load curtailments for power system and natural gas system are shown in Figs. 4 and 5 respectively. Curtailed loads decrease with the increase of hardened transmission lines. However, the slope for the reduction of load curtailments tend to be gentle. It means that more money will be cost to reduce the same quantities of load curtailments with the increase of hardened lines. Defenders must decide the optimal number of transmission lines to be hardened to gain more profit despite unpredictable attacks.

Fig. 4
figure 4

Load curtailment for IEEE 30-bus power system

Full size image
Fig. 5
figure 5

Load curtailments for 7-node natural gas system

Full size image

4.2 IEEE 118-bus system

This case is studied to test the validity of proposed method in large-scaled systems. It is based on a modified IEEE 118-bus power network and 14-node natural gas system.

The modified IEEE 118-bus power system depicted in Fig. 6 is consists of 54 thermal units including 42 fossil units and 12 natural gas-fired units, 186 branches, 9 transformers and 99 loads. Data for generators, buses, branches and load demands are from [23]. The 14-node natural gas system depicted in Fig. 7 has 3 natural gas suppliers, 12 pipelines, 2 compressors and 16 natural gas loads including 8 demand loads for gas-fired units. Parameters for natural gas network consults [5].

Fig. 6
figure 6

IEEE 118-bus power system

Full size image
Fig. 7
figure 7

14-node natural gas system

Full size image

We assume that attackers can damage up to 10 transmission lines in electricity network and 3 transmission lines in natural gas network. Defender has budget for hardening 0–10 power transmission lines and 0–3 gas transmission lines. Figures 8 and 9 depicts load curtailments in power system and natural gas network respectively. As 12 of the generators in power system are driven by gas fuel, natural gas network is quite important to the safe and stable operating of electricity system. It is obvious from the following figures that gas lines hardening plays a significant role in improving resilience of integrated energy system.

Fig. 8
figure 8

Load curtailments for IEEE 118-bus power system

Full size image
Fig. 9
figure 9

Load curtailments for 14-node natural gas system

Full size image

From the results shown in Figs. 8 and 9, we can see that two gas pipelines hardened and three gas pipelines hardened has nearly the same gas and power load curtailments. However, compared with none gas line hardened and one gas line hardened, load curtailments obviously decrease. Similarly, 7–10 power transmission lines hardened also have nearly the same results. As the costs will increase if more lines are hardened, 7 power transmission lines and 2 gas pipelines hardened is the most economic and effective plan to protect integrated energy system from attacks.

From the results for this large-scale system, we can also draw that our proposed model and method is feasible to improve and analyze resilience of integrated energy systems. Optimal plan for hardening power and gas lines can be derived from simulation results, which breaks the traditional view that resilience of system will improve as long as the number of hardened lines increases.

When applied for large-scale system or more defense and attack budgets, our proposed method is more time saving. Compared with column and constraint generation algorithm used in [16], computation time is much shorter as computational complexity increases. As stated in [24], C&CG algorithm has difficulties in operating in reasonable time when dealing with large-scale real problems. Therefore, algorithms based on benders decomposition or L-shape method should be developed [24]. The computation time is shown in Tables 6 and 7.

Table 6 Computation time for C&CG in [16] and the proposed method in the case of IEEE 30-bus system
Full size table
Table 7 Computation time for C&CG in [16] and our proposed method in the case of IEEE 118-bus system
Full size table

4.3 Impact of energy storage system

The concern on energy storage technologies is rapidly increasing due to high penetration of renewable and intermittent energy plug-in. As an effective manner to improve energy utilization efficiency, energy storage systems can solve the problem of mismatch between energy supply and demand side in time and space [25].

The system discussed in this paper integrates electricity and natural gas. Therefore, energy storage systems may include both electricity and gas storage systems. In power system, storage system is an effective manner to adjust peak. Excess power is stored into storage devices at off-peak period and stored energy will be transported back to grid when load demand is at peak. Moreover, as the unpredictable and intermittent nature of wind, solar and other renewable energy generation have great influence on resilience of power system, storage devices are usually installed nearby for tracking load changes [25]. However, electricity has the feature of easy to transport but difficult to store. Therefore, large-scale electricity storage technology is still a challenge. Compared with electricity, gas has the feature of easy to store but difficult to transport. Natural gas transportation costs mainly depend on the volume of gas supply and transport distance [26]. In fact, natural gas consumers are usually far away from gas sources and the cost for transporting is quite expansive, as a result of which the storage facilities for natural gas is of vital importance.

Either in power system or natural gas system, energy storage is an effective way to improve the resilience for the integrated system. We place 3 electricity storage devices at buses 21, 50, 96 in IEEE 118-bus power system and 2 gas storage devices at nodes 3, 5 in natural gas system. The capacity of each electricity storage devices is 50 MW and the capacity of each gas storage devices is 1000 kcf. The results considering energy storage systems is shown in Figs. 10 and 11.

Fig. 10
figure 10

Load curtailments for IEEE 118-bus power system with energy storage devices

Full size image
Fig. 11
figure 11

Load curtailments for 14-node gas system with energy storage devices

Full size image

Compared with Figs. 8 and 9, load curtailments in both electricity and natural system decrease apparently. With 7–10 electricity transmission lines and 2–3 gas transmission lines hardened, gas load curtailments are prevented thoroughly and electricity load curtailments reduced by nearly 80%. The energy storage system plays a significant role in improving resilience of integrated energy system.

In Section 2.1, we have stated that we assume natural gas system operates in steady state and line pack is ignored. Actually, large amounts of gas stored in pipelines, such as line-pack, can provide additional gas supply after the occurrence of contingencies. In this regard, the load shedding results obtained from steady-state model are conservative. Certain amounts of load shedding can be compensated by line-pack and hence can be avoided.

5 Conclusion

This paper proposes a robust optimization model for resilient operation of integrated energy system with electricity and natural gas infrastructures. The proposed model is formulated as a three-stage optimization problem which considers network reinforcement, damage caused by attackers and defenders response of both electricity and natural gas systems. Linearization technologies and decomposition algorithms are applied to reformulate and solve this defender-attacker-defender problem. Numerical results validate the effectiveness of our proposed model. Studies also point to the importance of energy storage systems in improving resilience of integrated energy system against contingencies.

In future works, we will focus on the resilience of energy systems which integrate other forms of energies and demand response management. And distribution networks or micogrids will also be concerned in future works.