*Article* **Demand Management for Resilience Enhancement of Integrated Energy Distribution System against Natural Disasters**

**Yuting Xu \*, Songsong Chen, Shiming Tian and Feixiang Gong**

China Electric Power Research Institute Co., Ltd., Beijing 100192, China; 15901168062@163.com (S.C.); oldtian@sina.com.cn (S.T.); gongfeixiangouc@126.com (F.G.)

**\*** Correspondence: yutingxu163@163.com

**Abstract:** For energy sustainability, the integrated energy distribution system (IEDS) is an efficient and clean energy system, which is based on the coordinated operation of a power distribution network, a gas distribution network and a district heating system. In this paper, considering the damage of natural disasters to IEDS, a demand management strategy is proposed to improve resilience of IEDS and ensure stable operation, which is divided into three stages. In the first stage, the electricity, natural gas and thermal energy are co-optimized in the simulating fault state to develop the importance ranking of transmission lines and gas pipelines. In the second stage, the natural disasters are classified as surface natural disasters and geological natural disasters. According to the types of natural disasters, the demand management strategy includes semi-emergency demand management scheme and full-emergency demand management scheme in the electrical resilience mode and the integrated resilience mode, respectively. In the third stage, the non-sequential Monte-Carlo simulation and scenario reduction algorithm are applied to describe potential natural disaster scenarios. According to the importance ranking of transmission lines and gas pipelines, a demand management strategy is formulated. Finally, the proposed strategy is applied on an IEEE 33-bus power system and a 19-node natural gas system. Its effectiveness is verified by numerical case studies.

**Keywords:** demand management; integrated energy distribution system; resilience; co-optimization; non-sequential Monte-Carlo simulation; scenario reduction algorithm

### **1. Introduction**

With the development of society and economy, predatory energy consumption has caused environmental pollution [1] and energy crisis [2]. The integrated energy distribution system (IEDS) [3,4] takes full account of electricity, natural gas, heat and other forms of energy coupling. It can achieve the effect of energy mutual benefit according to the energy consumption characteristics of electricity, natural gas and heat. Hence, IEDS is an efficient and clean energy system.

On the other hand, frequent natural disasters have severely affected the energy system. In 2011, the Great East Japan Earthquake caused power outages in 8.71 million homes in the affected area [5]. In 2012, approximately 7.5 million customers suffered power outages in the Hurricane Sandy in New York and the disaster caused an economic loss of 65 billion US dollars [6]. Many typhoons from the Pacific will land in China every year. For example, Jiangsu Province was hit by a typhoon in 2016, which caused two 500-kV transmission lines, four 220-kV transmission lines, and eight 110-kV transmission lines to trip and left many customers without power [7]. In this regard, it is necessary and exigent to enhance the IEDS resilience. The IEDS resilience is derived from the extension of power system resilience, which can be defined as the ability to anticipate, resist, absorb and recover from disruptions caused by extreme natural disasters such as earthquakes and hurricanes [8].

**Citation:** Xu, Y.; Chen, S.; Tian, S.; Gong, F. Demand Management for Resilience Enhancement of Integrated Energy Distribution System against Natural Disasters. *Sustainability* **2022**, *14*, 5. https://doi.org/10.3390/ su14010005

Academic Editors: Luis Hernández-Callejo, Sergio Nesmachnow and Sara Gallardo Saavedra

Received: 24 October 2021 Accepted: 30 November 2021 Published: 21 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

At present, the previous studies mainly focused on power system resilience [9–15]. Two types of strategies can be adopted to enhance resilience: operational measures [9–11] and hardening measures [12–15]. The operational measures include scheduling flexible back-up resources [9], using decentralized control strategies [10] and altering network topology [11]. The study in [9] proposes dispatching mobile emergency generators to restore critical loads and improve power distribution system resilience. In [10], the networked micro-grids (MGs) are scheduled by a decentralized control strategy, which can improve power quality by supporting and interchanging electricity among the networked MGs. In [11], power grid reconfiguration is adopted in an active distribution system to reduce load demand response and improve power grid resilience. On the other hand, the hardening measures mean physically enhancing the infrastructure to reduce its susceptibility to disruptions [12], In [13], a tri-level defender-attacker-defender (DAD) model is proposed and its result can provide the system hardening decisions.

Due to the increasing interaction among the electricity distribution system, natural gas distribution system and district heating system, it is not suitable to consider power system resilience solely and it is necessary to consider integrated energy system resilience. However, the literature contains little research on IEDS resilience [14,15]. The study in [14] considers that the overhead power grid can be hardened by replacing fragile overhead transmission lines with underground natural gas pipelines and proposes a two-stage robust model to formulate this issue. Combined with the natural gas system, the DAD model is further expanded in [15] to accommodate electricity and gas storage facilities. The existing enhanced resilience methods in the aforementioned research are listed in Table 1.


**Table 1.** The existing enhanced resilience methods.

Although the IEDS resilience has not been fully studied yet, there is much research focusing on the operation of integrated energy system [16–19]. In [16], an integrated framework based on the Newton-Raphson technique is proposed to solve the steady-state energy flow among electrical, natural gas, and district heating networks. Furthermore, the study in [17] proposes a coordinated optimal operation method of the regional energy internet, considering the combined cooling, heating and power (CCHP) units. The study in [18] proposes a robust day-ahead scheduling model for electricity and natural gas system, which minimizes the total cost including fuel cost, spinning reserve cost and cost of operational risk while ensuring the feasibility for all scenarios within the uncertainty set. In particular, the study in [19] proposes a novel linear method for Weymouth equation of natural gas system, which develops a more robust, flexible and tractable formulation of the integrated power and gas network.

The above research on the integrated energy system operation lay the foundation for IEDS resilience. Therefore, motivated by the aforementioned facts, this paper proposes a novel demand management strategy for improving the resilience of the power distribution network and gas distribution network in different types of natural disasters. The demand management strategy is divided into three stages. In the first stage, the critical transmission lines and gas pipelines are identified by co-optimizing electricity, natural gas and thermal energy under the assumptive fault. In other word, the importance ranking of transmission lines and gas pipelines are developed by simulating fault state. In the second stage, the natural disasters are classified as surface natural disasters (SNDs) and geological natural disasters (GNDs). The demand management strategy includes two modes. In the case of SNDs, the demand management strategy adopts the electrical resilience mode, in which the overhead transmission lines in power distribution system will be damaged, while the underground gas pipelines in gas distribution system still work without fault. In the case of GNDs, the demand management strategy adopts the integrated resilience mode, in which both overhead transmission lines and underground gas pipelines will be damaged. In the third stage, the non-sequential Monte-Carlo simulation and scenario reduction algorithm are applied to describe potential natural disaster scenarios. According to the importance ranking of transmission lines and gas pipelines, the demand management strategy is formulated. The advantages of the proposed strategy are demonstrated by numerical case studies.

The objective of the study is to propose a demand management strategy for resilience enhancement of IEDS against natural disasters and prove the effectiveness and practicality of this strategy. The major contributions of this paper are summarized as follows:

(1) Aiming at improving the IEDS resilience, this paper proposes a demand management strategy, which innovatively expands the traditional demand management in the power system.

(2) The demand management strategy includes semi-emergency demand management scheme and full-emergency demand management scheme in the electrical resilience mode and the integrated resilience mode, respectively. In the electrical resilience mode, gas distribution system can operate normally and play a role of energy storage to help power distribution system reduce demand response load. In the integrated resilience mode, the whole integrated energy system suffers damage and demand management is required.

(3) The demand management strategy is formulated in accordance with the importance of transmission lines and gas pipelines, taking into account the impact of the transmission grid and gas transmission network structure on demand management, which can increase the resilience of IEDS and reduce the economic subsidy cost of demand management.

The rest of this paper is organized as follows. The type of natural disasters, the stochastic approaches, and the framework of demand management strategy are discussed in Section 2. Then, Section 3 presents the scheduling model for the IEDS. In Section 4, numerical results are provided to validate the proposed strategy. Finally, Section 5 draws main conclusions in this paper.

### **2. The Demand Management Strategy for the IEDS under Natural Disasters**

### *2.1. The Type of Natural Disasters*

In this paper, natural disasters are divided into SNDs and GNDs according to whether the underground gas pipeline is destroyed or not. The SNDs include hurricane, typhoon and blizzard, while the GNDs include earthquake, landslide and mud-rock flow. The typical surface natural disaster and geological natural disaster are shown in Figure 1.

Because the overhead transmission lines are exposed to the environment, they are vulnerable to SNDs, such as typhoons and hurricanes. However, the underground gas pipelines have a certain resistibility to SNDs, so the gas distribution system can operate normally in the case of SNDs. Hence, the gas distribution system plays a role of energy storage to helps power distribution system improve resilience. In summary, the mode in which the overhead transmission lines are destroyed and the underground gas pipelines are not destroyed is called the electrical resilience mode.

On the other hand, both overhead transmission lines and underground gas pipelines may be damaged by the GNDs. In this case, the underground gas pipelines will be destroyed, and the gas distribution system also needs to be restored. Therefore, the mode in which the overhead transmission lines and the underground gas pipelines are destroyed is called the integrated resilience mode.

In this paper, the proposed demand management strategy can be applied to two types of natural disasters. The main difference between the two modes is whether the gas loads need to be reduced in the recovery process.

**Figure 1.** The typical surface natural disaster and geological natural disaster. (**a**) Satellite cloud picture of "Sangmei" typhoon; (**b**) Satellite telemetry picture of Wenchuan Earthquake.

### *2.2. The Fragility Curves and Stochastic Approaches*

The non-sequential Monte-Carlo simulation method generates many disaster scenarios with the random parameters of typhoon or earthquake. In order to get a typical disaster scenario, the scenario reduction algorithm is applied to reduce the number of samples. Then, as a calculation result of the scenario reduction algorithm, several main scenarios and their weights are derived. A detailed description of this algorithm can be found in [20].

The failure probabilities of the overhead transmission lines and underground gas pipelines in each simulated scenario are provided through their fragility curves, which express their failure probability as a function of the disaster parameter. The generic fragility curves of the overhead transmission lines and underground gas pipelines are shown in Figure 2 [8].

### *2.3. The Framework of the Demand Management Strategy*

As noted earlier, the demand management strategy is divided into three stages and its framework is depicted in Figure 3.

The first stage is the pre-disaster stage, in which the power distribution system, gas distribution system and district heating system co-optimize under the assumptive fault to identify critical transmission lines and gas pipelines, which is formulated in Section 3. Then, the second stage is mainly to distinguish the modes of demand management strategy

corresponding to different types of natural disasters. On this basis, the third stage generates the simulated scenarios and develops the demand management scheme.

**Figure 3.** The framework of the demand management strategy.

### **3. Model Formulation**

*3.1. Objective Function*

The goal of scheduling is to minimize the total operation cost of IEDS as shown in (1). Please refer to the Abbreviations for the explanation of the variables, indices, sets and parameters.

$$\begin{aligned} \min & \underbrace{\sum\_{s=1}^{N\_{\text{f}}} \lambda\_{s} \left( \sum\_{t=1}^{T} \underbrace{c\_{t}^{p} p\_{i^{\text{pr}},t}^{p^{\text{pr}}}}\_{i} + \underbrace{\sum\_{j^{\text{r}} \in \delta\_{\text{pr}}(b)} \epsilon\_{i}^{q} q\_{i^{\text{pr}},t,s}^{p^{\text{pr}}}}\_{\overline{\text{ii}}} + \underbrace{\sum\_{j^{\text{r}} \in \delta\_{\text{pr}}(b)} \epsilon\_{i}^{\text{r}} f\_{i^{\text{pr}},t,s}^{\text{dir}}}\_{\overline{\text{iii}}} + \underbrace{\sum\_{j^{\text{r}} \in \delta\_{\text{pr}}(b)} \epsilon\_{i}^{\text{f}} f\_{i^{\text{pr}},t,s}^{\text{dir}}}\_{\overline{\text{iii}}} + \underbrace{\sum\_{i^{\text{r}} \in \delta\_{\text{fr}}(b)} \epsilon\_{i}^{q} q\_{i^{\text{r}},t,s}^{\text{CH}}}\_{i^{\text{r}}} \right) \tag{1} \\ \underbrace{\sum\_{j^{\text{r}} \in \delta\_{\text{fr}}(b)} \epsilon\_{i}^{q} p\_{i^{\text{r}},t,s}^{\text{dir}}}\_{\overline{\text{i}}} + \underbrace{\sum\_{i^{\text{r}} \in \delta\_{\text{fr}}(b)} \epsilon\_{i}^{\text{f}} f\_{i^{\text{r}},t,s}^{\text{dir}}}\_{\overline{\text{iii}}} + \underbrace{\sum\_{i^{\text{r}} \in \delta\_{\text{fr}}(b)} \epsilon\_{i}^{q} p\_{i^{\text{r}},t,s}^{\text{LH}}}\_{i^{\text{r}}} \end{aligned}}\_{i^{\text{r}}} \tag{1}$$

The goal of the optimization is to minimize operation cost. The operation cost includes 9 parts, which can be divided into 3 categories. The first category is purchasing cost, including active power purchased (i) from electricity transmission system, reactive power purchased (ii) from electricity transmission system, natural gas purchased (iii) from gas transmission system, reactive power purchased (iv) from combined heating and power (CHP) units. The second category is demand response cost, including active power load demand response (v), reactive power load demand response (vi), gas load demand response (vii), and heat load demand response (viii). The third category is penalty cost, including penalty cost of wind power curtailment (ix).

### *3.2. The Power Distribution System Constraints*

The branch flow model is widely applied to compute the alternative current power flow of a power distribution system, but it is not suitable for optimization problems, because this model includes the power flow constraints, which is non-convex. For this difficulty, the second-order cone (SOC) relaxation is employed to relax the non-convex constraints in the branch flow model [21]. In detail, following the SOC relaxation method, the power flow equality constraints are replaced by inequality constraints and this modification still ensures the exactness, so the non-convex branch flow model is converted to the DistFlow model. Furthermore, the DistFlow model can be linearized to formulate microgrid operation constraints. In fact, the linear DistFlow model has been extensively used and justified in power distribution systems.

The linear DistFlow model gives constraints (2)–(13).

$$\begin{array}{cccc}\sum\limits\_{\mathbf{i}^{p}\in\delta\_{P}(\mathbf{b})}p^{pur}\_{\mathbf{i}^{p}\mathbf{i},t,s}+\sum\limits\_{\mathbf{i}^{l\varepsilon}\in\delta\_{I\varepsilon}(\mathbf{b})}p^{l}\_{\mathbf{i}^{l}\mathbf{i}^{l},t,s}+\sum\limits\_{\mathbf{i}^{\mathrm{CHP}}\in\delta\_{\mathrm{CHP}}(\mathbf{b})}p^{\mathrm{CHP}}\_{\mathbf{i}^{\mathrm{CHP}},t,s}+\sum\limits\_{\mathbf{i}^{\mathrm{w}}\in\delta\_{\mathrm{W}}(\mathbf{b})}p^{\mathrm{w}}\_{\mathbf{i}^{\mathrm{w}},t,s}\\\vdots&\vdots&=I(\mathbf{i}^{1\mathrm{c}})\\\mathbf{i}^{1\mathrm{s}}\in\delta\_{ls}(\mathbf{b})\\\mathbf{i}^{1\mathrm{c}}\in I(\mathbf{i}^{1\mathrm{s}})\end{array}\tag{2}$$

$$\begin{array}{ccccc}\displaystyle\sum\_{\mathbf{i}^{\mathsf{lrs}}\in\mathcal{S}\_{\mathsf{pr}}(\mathbf{b})}q^{\mathrm{pure}}\_{\mathbf{i}^{\mathsf{lrs}},\mathbf{t},\mathbf{s}}+&\sum\_{\mathbf{i}^{\mathsf{lrs}}\in\mathcal{S}\_{\mathsf{lrs}}(\mathbf{b})}q^{\mathrm{I}}\_{\mathbf{i}^{\mathsf{lrs}}\mathbf{j}^{\mathsf{lrs}},\mathbf{t},\mathbf{s}}+&\sum\_{\mathbf{i}^{\mathsf{lrs}}\in\mathcal{S}\_{\mathsf{l}^{\mathsf{l}}\mathbf{H}^{\mathsf{l}}}(\mathbf{b})}q^{\mathrm{C}\mathbf{H}\mathbf{P}}\_{\mathbf{i}^{\mathsf{lrs}}\mathbf{H}^{\mathsf{l}},\mathbf{t},\mathbf{s}}+&\sum\_{\mathbf{i}^{\mathsf{lrs}}\in\mathcal{S}\_{\mathsf{l}^{\mathsf{l}}\mathbf{s}}(\mathbf{b})}q^{\mathrm{W}}\_{\mathbf{i}^{\mathsf{lrs}},\mathbf{t},\mathbf{s}}\\&\displaystyle&\mathbf{j}^{\mathsf{l}\mathbf{s}}\in\mathcal{I}\left(\mathbf{i}^{\mathsf{l}\mathbf{s}}\right)\\&\displaystyle&\mathbf{k}^{\mathsf{l}\mathbf{s}}\in\mathcal{I}\left(\mathbf{i}^{\mathsf{l}\mathbf{s}}\right)\\&\displaystyle&\mathbf{k}^{\mathsf{l}\mathbf{s}}\text{ and }&\mathbf{k}^{\mathsf{l}\mathbf{s}}\in\mathcal{I}\left(\mathbf{k}^{\mathsf{l}}\right)\end{array}$$

$$\mathbf{x} = \sum\_{\substack{\mathbf{i}^{ls} \in \mathcal{S}\_{ls}(\mathbf{b}) \\ \mathbf{j}^{lc} \in I(\mathbf{i}^{ls})}} q\_{l^{ls}}^{l'} \mathbf{i}\_{\mathbf{i}^{ls},\mathbf{j}^{lc},\mathbf{s}} + q\_{b,t,s}^{l'oad} - q\_{b,t,s}^{demand} \qquad \forall b, \forall t, \forall s \tag{3}$$

$$V\_{i^{l\*},t,s} - V\_{j^{l\*},t,s} = \frac{p\_{i^{l\*} \circ t, t, s}^{l} r\_{i^{l\*} \circ t} + q\_{i^{l\*} \circ t, s}^{l} \mathbf{x}\_{i^{l\*} \circ t}}{V\_{\text{base}}} \qquad \dot{\mathbf{r}}^{\text{l}s} \in \delta\_{ls}(b), \mathbf{y}^{\text{l}t} \in l(\mathbf{i}^{\text{l\*}}), \forall t, \forall s \tag{4}$$

$$1 - \overline{p}\_{i^{l s} j^{l s}}^l \le \overline{p}\_{i^{l s} j^{l s}, t, s}^l \le \overline{p}\_{i^{l s} j^{l s}}^l \qquad i^{l s} \in \delta\_{ls}(b), j^{l a} \in l(i^{l s}), \forall t, \forall s \tag{5}$$

$$\mathbf{I} - \overline{\eta}^{l}\_{l^{ls}\dot{f}^{lc}} \le \eta^{l}\_{l^{ls}\dot{f}^{lc}, l^{ls}s} \le \overline{\eta}^{l}\_{l^{ls}\dot{f}^{lc}} \qquad \mathrm{i}^{ls} \in \delta\_{ls}(b), \mathbf{j}^{l\varepsilon} \in \mathbf{l}(\mathbf{i}^{ls}), \forall t, \forall s \tag{6}$$

$$0 \le p\_{i^{\rm pr}, t, s}^{\rm pur} \le \overline{p}^{\rm pr} \qquad i^{\rm pr} \in \delta\_{\rm pc}(b), \forall t, \forall s \tag{7}$$

$$0 \le \overline{q}\_{i^{\rm pr},t,\mathbf{s}}^{\rm prr} \le \overline{q}^{\rm prr} \qquad i^{\rm pr} \in \mathcal{S}\_{\rm pr}(\mathbf{b})\_{\prime} \forall \mathbf{t}, \forall \mathbf{s} \tag{8}$$

$$
\underline{V}\_b \le V\_{b,t,s} \le \nabla\_b \qquad \forall b \; \forall t \; \forall s \tag{9}
$$

$$0 \le p\_{b,t,s}^{demand} \le p\_{b,t,s}^{load} \qquad \forall b \text{ } \forall t \text{ } \forall s \tag{10}$$

$$0 \le q\_{b,t,s}^{demand} \le q\_{b,t,s}^{load} \qquad \forall b \text{ } \forall t \text{ } \forall s \tag{11}$$

$$p\_{i^w, t, s}^{w} = \overline{p}\_{i^w, t, s}^{w} - p\_{i^w, t, s}^{w, cut} \qquad i^w \in \delta\_{w}(b), \forall t, \forall s \tag{12}$$

$$q\_{i^w, t, s}^w = \tan(\arccos(\phi\_{i^w})) \, p\_{i^w, t, s}^w \qquad \text{i^w} \in \delta\_w(b), \forall t, \forall s \tag{13}$$

Specifically, (2) and (3) imply active power balance and reactive power balance, which means that the active (or reactive) power injection amount is equivalent to outflow amount at each bus. (4) is the DistFlow equation, which relates the active power and reactive power flows of a transmission line to the voltage magnitude of its two terminal buses. (5) and (6) denote active power limit and reactive power limit on overhead transmission lines, respectively. (7) and (8) refer to the active (or reactive) power purchasing limits. (9) gives voltage limit at each bus. (10) and (11) imply the load demand response of active power and reactive power, respectively. (12) and (13) describe power outputs of wind generators.

### *3.3. The Gas Distribution System Constraints*

The operation constraints of the gas distribution system are shown in (14)–(19).

$$\begin{array}{cccc}\sum\_{\boldsymbol{i}^{\rm FI}\in\mathcal{S}\_{\rm FI}(\boldsymbol{n})}f^{\rm pur}\_{\boldsymbol{i}^{\rm FI},\boldsymbol{t},\boldsymbol{s}}+&\sum\_{\boldsymbol{i}^{\rm SC}\in\mathcal{S}\_{\rm SC}(\boldsymbol{n})}f^{\rm S}\_{\boldsymbol{i}\boldsymbol{i}^{\rm S},\boldsymbol{t}^{\rm s},\boldsymbol{t},\boldsymbol{s}}&=&\sum\_{\boldsymbol{i}^{\rm S}\in\mathcal{S}\_{\rm SC}(\boldsymbol{n})}f^{\rm S}\_{\boldsymbol{i}^{\rm S},\boldsymbol{i}^{\rm res},\boldsymbol{t},\boldsymbol{s}}\\&\boldsymbol{j}^{\rm S^{\rm s}}\in\mathcal{G}(\boldsymbol{i}^{\rm c^{\rm c}})&\boldsymbol{j}^{\rm c^{\rm c}}\in\mathcal{G}(\boldsymbol{i}^{\rm c^{\rm c}})\\&+\sum\_{\boldsymbol{i}^{\rm CHP}\in\mathcal{S}\_{\rm CHP}(\boldsymbol{n})}f^{\rm CHP}\_{\boldsymbol{i}^{\rm s},\boldsymbol{t},\boldsymbol{s}}+&\sum\_{\boldsymbol{i}^{\rm iso}\in\mathcal{S}\_{\rm Ni}(\boldsymbol{n})}f^{\rm iso}\_{\boldsymbol{i}^{\rm iso},\boldsymbol{t},\boldsymbol{s}}+f^{\rm load}\_{\boldsymbol{n},\boldsymbol{t},\boldsymbol{s}}-f^{\rm demand}\_{\boldsymbol{n},\boldsymbol{t},\boldsymbol{s}}&\quad\forall\boldsymbol{n},\,\forall\boldsymbol{t},\,\forall\boldsymbol{s}\end{array}\tag{14}$$

$$-\overline{f}^{\mathbb{R}}\_{i\mathbb{R}^{\mathfrak{s}}/\mathbb{R}^{\mathfrak{s}}} \le f^{\mathbb{R}}\_{i\mathbb{R}^{\mathfrak{s}}/\mathbb{R}^{\mathfrak{s}},\mathfrak{t},\mathfrak{s}} \le \overline{f}^{\mathbb{R}}\_{i\mathbb{R}^{\mathfrak{s}}/\mathbb{R}^{\mathfrak{s}}} \qquad i\mathbb{S}^{\mathfrak{s}} \in \delta\_{\mathbb{S}^{\mathfrak{s}}}(\mathfrak{n}), j^{\mathbb{S}^{\mathfrak{s}}} \in l(i^{\mathbb{S}^{\mathfrak{s}}}), \forall t, \forall \mathbf{s} \tag{15}$$

$$0 \le f\_{\overline{\mathcal{P}}^{\mathfrak{P}\mathfrak{F}},t,s}^{\text{pure}} \le \overline{f}^{\mathfrak{P}^{\mathfrak{p}\mathfrak{r}}} \qquad i^{\mathfrak{p}\mathfrak{g}} \in \delta\_{\mathfrak{p}\mathfrak{g}}(\mathfrak{n})\_{\mathfrak{r}} \,\forall t\_{\prime} \,\forall s \tag{16}$$

$$
\Xi\_{n} \leq \Xi\_{n,t,s} \leq \Xi\_{n} \qquad \forall n \, \forall t \, \forall s \tag{17}
$$

$$0 \le f\_{n,t,s}^{dcmd} \le f\_{n,t,s}^{load} \quad \forall n, \forall t, \forall s \tag{18}$$

$$f^{\mathcal{S}}\_{\text{i}\text{i}^{p}\text{j}^{\text{e}},\text{t},\text{s}} = K\_{\text{i}\text{i}^{p}\text{j}^{\text{e}}} \sqrt{\beta (z\_{\text{i}\text{v}}, z\_{\text{j}\text{v}}, z\_{\text{j}\text{e}^{\text{e}},\text{t},\text{s}}) (z\_{\text{i}\text{v}}^{2}, -z\_{\text{j}\text{v}}^{2})} \qquad i\text{\textquotedblleft} \in \delta\_{\text{\mathcal{S}}^{\text{e}}}(\text{n}), j\text{\textquotedblright} \in l(\text{i}\text{\textquotedblleft}), \forall \text{t}, \forall \text{s} \quad \text{(19)} \right)$$

(14) denotes natural gas balance, which means that the natural gas flow injection amount is equivalent to amount at each node. (15) refers to gas flow limit in gas pipelines. (16) refers to the active natural gas purchasing limits. In order to avoid damage or malfunction of gas pipeline caused by too low or too high gas pressure, (17) gives gas pressure limit at each node. (18) implies load demand response of natural gas. (19) is the Weymouth equation, which relates gas flow in gas pipeline and pressure of its two terminal nodes. Specifically, the direction coefficient *β*(*zi gs*,*t*,*s*, *zj ge*,*t*,*s*) describes the relationship between gas flow direction and gas pressure value at both ends of gas pipeline, which is denoted as follows:

$$\beta(z\_{i^{\Re},t,s}, z\_{j^{\Re},t,s}) = \begin{cases} 1\_{i^{\Re},t,s} \ge z\_{j^{\Re},t,s} \\ -1\_{i^{\Re},t,s} < z\_{j^{\Re},t,s} \end{cases} \quad \bar{r}^{\circledast} \in \delta\_{\mathbb{S}^3}(n), j^{\Re} \in l(i^{\Re}), \forall t, \forall s \tag{20}$$

(20) denotes the fact that gas flow runs from higher pressure to lower pressure. Moreover, the Weymouth equation (19) can be rewritten as (21).

$$\left(f\_{\left[\delta^{\mathcal{S}}\right]^{\mathcal{B}^{\mathcal{S}}},t,\mathsf{s}}^{\mathcal{S}}\right)^{2} = \left\{\mathfrak{f}(z\_{\delta^{\mathcal{S}},\mathsf{s},\mathsf{s}},z\_{\beta^{\mathcal{S}},\mathsf{t},\mathsf{s}})\left(\mathcal{K}\_{\delta^{\mathcal{S}}\beta^{\mathcal{S}}}\right)^{2}\left(z\_{\delta^{\mathcal{S}},\mathsf{t},\mathsf{s}}^{2} - z\_{\beta^{\mathcal{S}},\mathsf{t},\mathsf{s}}^{2}\right) \quad \text{ if } \mathsf{s}^{\mathcal{S}} \in \delta\_{\mathcal{S}^{\mathcal{S}}}(\mathsf{n}),\,\mathsf{j}\beta^{\mathcal{s}} \in I(i\mathbb{\delta}^{\mathcal{S}}),\,\forall t,\forall\mathsf{s} \tag{21}$$

Apparently, since gas flow is squared, its direction cannot connect to gas pressure value. In other words, the direction coefficient is useless in (21). For the flow direction problem, (21) is equivalently reformulated as (22)–(25) by the Big-M theory with additional binary variables [22].

$$1 - M(1 - \delta\_{\vec{\nu}^y|\vec{\nu}^x, t, s}) \le \left(f\_{\vec{\nu}^y|\vec{\nu}^x, t, s}^{\vec{\nu}}\right)^2 - \left(K\_{\vec{\nu}^y|\vec{\nu}^x}\right)^2 \left(z\_{\vec{\nu}^y, t, s}^2 - z\_{\vec{\nu}^x, t, s}^2\right) \le M(1 - \delta\_{\vec{\nu}^y|\vec{\nu}^x, t, s}) i^{y\varepsilon} \in \delta\_{\vec{\nu}^y}(\mathfrak{n}), j^{y\varepsilon} \in l(i^{\infty}), \forall t, \forall s \tag{22}$$

$$1 - M\delta\_{\mathcal{V}^{\mathcal{B}}|\mathcal{F},t,\mathbf{s}} \le \left(f^{\mathcal{S}}\_{\mathcal{V}^{\mathcal{B}}|\mathcal{F},t,\mathbf{s}}\right)^{2} - \left(K\_{\mathcal{V}^{\mathcal{B}}|\mathcal{F}}\right)^{2} \left(z^{2}\_{\mathcal{V}^{\mathcal{B}},t,\mathbf{s}} - z^{2}\_{\mathcal{V}^{\mathcal{B}},t,\mathbf{s}}\right) \le M\delta\_{\mathcal{V}^{\mathcal{B}}|\mathcal{F},t,\mathbf{s}} \qquad \vec{\mathbf{r}}^{\mathcal{S}} \in \delta\_{\mathcal{S}^{\mathcal{S}}}(\mathbf{n}), \vec{\mathbf{y}}^{\mathcal{B}} \in l(\vec{\mathbf{r}}^{\mathcal{B}}), \forall t,\forall \mathbf{s} \tag{23}$$

$$-M(1 - \delta\_{\mathbb{N}^s|\mathbb{N}^s, t, \mathbf{s}}) \le z\_{\mathbb{N}^s, t, \mathbf{s}} - z\_{\mathbb{j}^{\text{tr}}, t, \mathbf{s}} \le M \delta\_{\mathbb{N}^s|\mathbb{N}^{\text{tr}}, t, \mathbf{s}} \qquad i^{\text{gs}} \in \delta\_{\mathbb{S}^\text{s}}(n), j^{\text{gc}} \in l(i^{\text{gs}}), \forall t, \forall \mathbf{s} \tag{24}$$

$$-M(1 - \delta\_{i^{\rm g}j^{\rm gr}, t, s}) \le f^{\rm g}\_{i^{\rm g}j^{\rm gr}, t, s} \le M\delta\_{i^{\rm g}j^{\rm gr}, t, s} \qquad i^{\rm g} \in \delta\_{\rm g^{\rm g}}(n), j^{\rm g} \in I(i^{\rm gs}), \forall t, \forall s \tag{25}$$

It is worth mentioning that the quadratic terms in (22) and (23) are nonlinear, which makes the optimization model difficult to solve, so the quadratic terms are linearized by the piecewise linear approximation method. The gas pressure square magnitude (i.e.,*z*<sup>2</sup> *i gs*,*t*,*s*) in (19) can be directly replaced by a positive continuous variable *Zi gs*,*t*,*s*. However, the same approach cannot be applied to tackle gas flow square magnitude, because the linear terms of gas flow are expressed in natural gas balance constraint (14) and gas pipelines limit constraint (15). Therefore, as an efficient piecewise linear approximation method, the incremental linearization method [19] is employed to deal with the nonlinearity of gas flow square magnitude, which is appropriate for mixed-integer linear programming.

The domain of gas flow in pipeline is *f g i gsj ge*,*t*,*s* ∈ −*f g i gsj ge*, *f g i gsj ge* , which is denoted by (15). The domain −*f g i gsj ge*, *f g i gsj ge* is divided into *KP* intervals by breakpoints *xi gsj ge*,*<sup>k</sup>* as shown in (26):

$$\cdots \overline{f}\_{j\mathbb{P}^{\mathbb{P}}\mathbb{P}}^{\mathbb{R}} = \mathbf{x}\_{i\mathbb{P}^{\mathbb{P}}}, 1 \le \mathbf{x}\_{i\mathbb{P}^{\mathbb{P}}} \mathbf{z} \le \cdots \le \mathbf{x}\_{i\mathbb{P}^{\mathbb{P}}} \mathbf{z}^{\mathbf{r}} \le \cdots \le \mathbf{x}\_{i\mathbb{P}^{\mathbb{P}}} \mathbf{z}^{\mathbf{r}} \\ \mathbf{x}^{\mathbf{r}} \le \mathbf{x}\_{i\mathbb{P}^{\mathbb{P}}\mathbb{X}} \mathbf{z}^{\mathbf{r}} + 1 = \overline{f}\_{i\mathbb{P}^{\mathbb{P}}\mathbb{P}}^{\mathbb{R}} \qquad i\mathbb{P}^{\mathbb{P}} \in \delta\_{\mathbb{R}^{r}}(\mathbf{n}), \\ i^{\otimes \mathbf{r}} \in \mathbb{I}(i^{\otimes \mathbf{r}}) \tag{26}$$

$$\varepsilon\_{i^{j\_{\mathcal{F}}}j^{\varepsilon}} = \mathbf{x}\_{i^{j\_{\mathcal{F}}}j^{\varepsilon},k+1} - \mathbf{x}\_{i^{j\_{\mathcal{F}}}j^{\varepsilon},k} \qquad i^{\mathcal{G}} \in \delta\_{\mathcal{G}}(n), j^{\mathcal{G}} \in l(i^{\mathcal{G}}), k = 1, 2, \dots, KP-1 \tag{27}$$

$$KP = \left[2\overline{f}^{\mathbb{g}^{\varepsilon}}\_{i^{\mathbb{g}^{\varepsilon}}j^{\mathbb{g}^{\varepsilon}}} / \varepsilon\_{i\mathbb{f}^{\varepsilon}\bar{j}^{\mathbb{g}^{\varepsilon}}}\right] \qquad i^{\mathbb{g}^{\varepsilon}} \in \delta\_{\mathbb{g}^{\varepsilon}}(n), j^{\mathbb{g}^{\varepsilon}} \in l(i^{\mathbb{g}^{\varepsilon}}) \tag{28}$$

The size of interval is calculated by (27) and the number of intervals is calculated by (28). Corresponding to each *xi gsj ge*,*k*, the ordinate value is calculated by *yi gsj ge*,*<sup>k</sup>* = *x*<sup>2</sup> *i gsj ge*,*k*. Figure 4 shows the piecewise linearization of nonlinear function. The quadratic function (solid red line) is approximated by piecewise linear function (solid blue line).

**Figure 4.** The piecewise linearization of nonlinear function.

Then, the linear terms of (22) and (23) are formulated by

$$\begin{split} \left( \left( \mathbf{K}\_{i\mathbb{P}^{\mathcal{I}}\boldsymbol{j}\mathbb{R}^{\mathsf{c}}} \right)^{2} \left( \mathbf{Z}\_{i\mathbb{P}^{\mathcal{I}}\boldsymbol{j},\mathbb{R}} - \mathbf{Z}\_{j\mathbb{P}^{\mathcal{I}},\mathbb{I},\mathsf{s}} \right) = y\_{i\mathbb{P}^{\mathcal{I}}\boldsymbol{j}\mathbb{R}^{\mathsf{c}},\mathbbm{1}} + \sum\_{k=1}^{\mathsf{KP}} \left( y\_{i\mathbb{P}^{\mathcal{I}}\boldsymbol{j}\mathbb{R}^{\mathsf{c}},k+1} - y\_{i\mathbb{P}^{\mathcal{I}}\boldsymbol{j}\mathbb{R}^{\mathsf{c}}} \right) \mu\_{i\mathbb{P}^{\mathcal{I}}\boldsymbol{j}\mathbb{R}^{\mathsf{c}},\mathsf{t},\mathsf{s},\mathsf{k}} \\ \mathbf{i\mathbb{P}^{\mathsf{c}}} \in \delta\_{\mathcal{S}^{\mathsf{c}}}(\mathsf{n})\_{\boldsymbol{\tau}} \boldsymbol{j}\mathbb{R}^{\mathsf{c}} \in \mathrm{I} \left( \mathbf{i}\mathbb{P}^{\mathcal{R}} \right), \forall \boldsymbol{t}, \forall \boldsymbol{s},\boldsymbol{k} = \mathbf{1}, \mathsf{2}, \mathsf{c} \cdot \mathsf{N} - 1 \end{split} \tag{29}$$

$$\begin{split} f^{\mathcal{S}}\_{i\bar{i}^{\mathcal{P}}j^{\mathcal{S}},t,s} &= \mathbf{x}\_{i\bar{i}^{\mathcal{P}}j^{\mathcal{P}},1} + \sum\_{k=1}^{\mathrm{KP}} \left( \mathbf{x}\_{i\bar{i}^{\mathcal{P}}j^{\mathcal{S}},k+1} - \mathbf{x}\_{i\bar{i}^{\mathcal{P}}j^{\mathcal{P}},k} \right) \mu\_{i\bar{i}^{\mathcal{S}}j^{\mathcal{P}},t,s,k} \\ &i^{\mathcal{S}^{\mathrm{s}}} \in \delta\_{\mathcal{S}^{\mathrm{s}}}(n), j^{\mathcal{S}^{\mathrm{e}}} \in l(i^{\mathcal{S}^{\mathrm{s}}}), \forall t, \forall s, k = 1, 2, \dots, KP - 1 \end{split} \tag{30}$$

$$0 \le \mu\_{i^{\otimes t} j^{\otimes t}, t, s, k} \le 1 \qquad i^{\otimes s} \in \delta\_{\mathbb{S}^\otimes}(n), j^{\otimes t} \in I(i^{\otimes s}), \forall t, \forall s, \forall k \tag{31}$$

$$\mu\_{\mathbb{R}^6} \mu\_{\mathbb{R}^6, t, \succ k+1} \le \mathbb{Z}\_{\mathbb{R}^6} \mu\_{\mathbb{R}^4, t, \succ k} \le \mu\_{\mathbb{R}^6} \mu\_{\mathbb{R}^5, t, \succ k} \qquad \mathbb{R}^{\mathbb{R}} \in \delta\_{\mathbb{R}^5}(n), j^{\otimes t} \in l(\mathbb{R}^6), \forall t, \forall s, k = 1, 2, \dots, KP-1 \tag{32}$$

where *μ<sup>i</sup> gsj ge*,*t*,*s*,*<sup>k</sup>* and *ξ<sup>i</sup> gsj ge*,*t*,*s*,*<sup>k</sup>* are auxiliary variables at *k*th interval of piecewise linear function. Constraints (32) ensures that there is at most one index *k* of an interval with 0 < *μ<sup>i</sup> gsj ge*,*t*,*s*,*<sup>k</sup>* < 1 and other indexes are equal to 0 or 1. Specifically, *μ<sup>i</sup> gsj ge*,*t*,*s*,*k*+<sup>1</sup> = 0, if *ξi gsj ge*,*t*,*s*,*<sup>k</sup>* = 0; *μ<sup>i</sup> gsj ge*,*t*,*s*,*k*+<sup>1</sup> = 1, if *ξ<sup>i</sup> gsj ge*,*t*,*s*,*<sup>k</sup>* = 1.

From the above, the operation constraints of the gas distribution system consists of (14)–(19) and the nonlinear constraint (19) is reformulated by (22)–(25) and (29)–(32).

### *3.4. The District Heating System Constraints*

The district heating system constraints consist of heat balance and operation of heat generating equipment, which includes CHP units, electrical heat pumps (EHPs) and boilers.

The heat balance constraint:

$$\sum\_{i^\* \in \mathcal{H}^p \in \delta\_{\mathcal{C}HP}(n)} h^{\mathcal{C}HP}\_{i^\* \mathcal{C}HP} + \sum\_{i^\* \in \mathcal{H}^p \in \delta\_{\mathcal{E}HP}(b)} h^{\mathcal{E}HP}\_{i^\* \mathcal{C}HP} + \sum\_{i^{\text{pair}} \in \delta\_{\mathcal{S}wi}(n)} h^{\text{pair}}\_{i^{\text{pair}} \mathcal{J}, \mathfrak{s}} = h^{\text{load}}\_{t, \mathfrak{s}} - h^{\text{demand}}\_{t, \mathfrak{s}} \quad \forall t, \forall \mathbf{s} \tag{33}$$

$$0 \le h\_{t,s}^{demand} \le h\_{t,s}^{load} \qquad \forall t, \forall s \tag{34}$$

The constraints (33) means that the heat generating amount is equivalent to heat demand amount in the district heating system. (34) implies load demand response of heat.

The CHP units constraints are as follows:

$$\underline{p}\_{i^{\rm CHP}}^{\rm CHP} \le p\_{i^{\rm CHP}, t\_{\rm r^\*}}^{\rm CHP} \le \overline{p}\_{i^{\rm CHP}}^{\rm CHP} \qquad i^{\rm CHP} \in \delta\_{\rm CHP}(b), \forall t\_{\rm r} \forall s \tag{35}$$

$$
\underline{q}\_{i^{\rm CHP}}^{\rm CHP} \le \eta\_{i^{\rm CHP}, t, \rm s}^{\rm CHP} \le \overline{q}\_{i^{\rm CHP}}^{\rm CHP} \qquad i^{\rm CHP} \in \delta\_{\rm CHP}(b), \forall t, \forall s \tag{36}
$$

$$p\_{i^{\rm CHP}\_{\rm r},t,s}^{\rm CHP} = \eta\_{i^{\rm CHP}}^{\rm CHP,clle} f\_{i^{\rm CHP}\_{\rm r},t,s}^{\rm CHP} \qquad i^{\rm CHP} \in \delta\_{\rm CHP}(b)\_{\prime} \,\forall t, \,\forall s \tag{37}$$

$$h\_{i^{\rm CHP},t,s}^{\rm CHP} = \eta\_{i^{\rm CHP}}^{\rm CHP,h\rm heat} f\_{i^{\rm CHP},t,s}^{\rm CHP} \qquad \text{i}^{\rm CHP} \in \delta\_{\rm CHP}(b)\_{\prime} \forall t, \forall s \tag{38}$$

(35) and (36) denote output limit of CHP units, which means that active power and reactive power must be between the minimum output and maximum output. (37) and (38) imply power generation efficiency and heat production efficiency of CHP units, respectively.

The EHPs constraints are as follows:

$$0 \le p\_{i^{EHP}, t, s}^{EHP} \le \overline{p}\_{i^{EHP}}^{EHP} \qquad \text{i}^{EHP} \in \delta\_{EHP}(b), \forall t, \forall s \tag{39}$$

$$h\_{i^{\rm EHP},t,\mathbf{s}}^{\rm EHP} = \eta\_{i^{\rm EHP}}^{\rm EHP} p\_{i^{\rm EHP},t,\mathbf{s}}^{\rm EHP} \qquad i^{\rm EHP} \in \mathcal{S}\_{\rm EHP}(\mathbf{b})\_{\star} \,\forall t, \forall \mathbf{s} \tag{40}$$

(39) denotes power consumption limit of EHP. (40) refers to heat production efficiency of EHP.

The boilers constraints are as follows:

$$0 \le f\_{i^{\text{bvi}}, t, s}^{\text{bvi}} \le \overline{f}^{\text{bvi}} \quad {}\_{i^{\text{bvi}}} \text{ $i^{\text{bvi}}$ } \in \mathcal{S}\_{\text{bvi}}(n) \text{ } \forall \text{t} \text{ } \forall \text{s} \tag{41}$$

$$h\_{i^{bwi},t\_{\prime}s}^{bwi} = \eta\_{i^{bwi}}^{bwi} f\_{i^{bwi},t\_{\prime}s}^{bwi} \qquad i^{bwi} \in \delta\_{bwi}(n)\_{\prime} \forall t\_{\prime} \forall s \tag{42}$$

(41) denotes natural gas consumption limit of boiler. (42) refers to heat production efficiency of boiler.

The proposed scheduling model of IEDS is a mixed integer linear programming (MILP), which can be solved with the solvers, such as CPLEX, GUROBI.

### **4. Numerical Results**

### *4.1. The Operation of IEDS*

The proposed hardening strategy is examined on a test system, which consists of IEEE 33-bus power distribution system, 19-node gas distribution system and district heating system, shown in Figure 5. The numbers without brackets in Figure 5 represent the node number, and the numbers with brackets represent the line number.

**Figure 5.** The integrated energy distribution system.

As shown in Figure 5, there are three CHP units, two EHPs, two boilers and three wind generators. In fact, the three energy systems are geographically overlapping. For convenience, Figure 5 is a schematic diagram of the wiring of three systems.

First of all, in order to facilitate the explanation of the test system, the operation of IEDS is introduced in normal state, which is the basis of fault operation. The supply of active power load, gas load and heat load are shown in Figure 6.

As shown in Figure 6, the supply and demand of IEDS are balanced at each time slot in normal operation. However, if the overhead transmission lines or underground gas pipelines are broken, the energy supply will be restricted and load demand will not be met, so an effective demand management strategy is urgently needed to alleviate the tight supply of IEDS.

### *4.2. The Importance Ranking of Transmission Lines and Pipelines*

After breaking transmission lines and pipelines, the amount of load demand response is calculated at each simulating fault state. This section calculates the economic subsidy caused by load demand response. The corresponding economic subsidy caused by the damaged transmission lines and pipelines is shown in Tables 2 and 3.

Then, according to the economic subsidy in each fault state, the importance of transmission lines and pipelines is ranked, which is divided into three levels. In the A-level and B-level, economic subsidy is greater than 120 K and 60 K, respectively. Economic subsidy less than 60 K is the C-level.

**Figure 6.** The supply of load in normal state. (**a**) The supply of active power load; (**b**) The supply of gas load; (**c**) The supply of heat load.

### *4.3. The Electrical Resilience Mode and Integrated Resilience Mode*

In the electrical resilience mode, the surface natural disasters are simulated. In total, 100 fault scenarios are generated by non-sequential Monte-Carlo simulation, and then 5 main scenarios are obtained by scenario reduction algorithm, which are used to calculate the fault operation conditions caused by the surface natural disasters. The main scenario with the largest probability is 0.34. Transmission lines No. 6, 14, 17, 19, 29, and 30 have failed. After the failure, the power load cannot be fully supplied. Comparing the normal operation state, a demand response is required, as shown in Figure 7.


**Table 2.** The economic subsidy of the damaged transmission lines.

**Table 3.** The economic subsidy of the damaged pipelines.


**Figure 7.** The power load supply situation. (**a**) The power load supply situation in normal operation state; (**b**) The power load supply situation in failure operation state.

As shown in Figure 7, destruction of transmission lines results in an amount of power loads that cannot obtain electricity through outsourcing. The wind turbines are damaged due to natural disasters and the output decreases. At this time, the distribution network forms many small island grids, and CHP provides power supply to the island grids.

However, affected by the supply capacity of CHP, a certain amount of power load is still needed as a demand response resource to participate in demand management. This demand management scheme that considers the support of the natural gas system to the distribution network is called a semi-emergency scheme in this paper, and it is not necessary to use all demand response capabilities.

In the integrated resilience mode, the geological natural disasters are simulated. Similarly, 100 failure scenarios are generated and reduced to 5 main scenarios to analyze the operation of IEDS after geological natural disaster. The main scenario with the largest probability is 0.29. Transmission lines No. 7, 9, 14, 16, 20, 26, and 29 have failed. Gas pipelines No. 4, 8, 11, and 15 have failed. After the failure, not only can the power load not be fully supplied, but the natural gas is also unable to supply all the gas load due to the failure of some pipelines. The power load and gas load supply conditions are shown in Figure 8.

**Figure 8.** The power load and gas load supply situation. (**a**) The power load supply situation in failure operation state; (**b**) The gas load supply situation in failure operation state.

As shown in Figure 8, after the natural gas distribution network was destroyed by geological natural disasters, some gas pipeline failures caused a gap in natural gas supply. The red shaded part in Figure 8b is the amount of natural gas involved in gas load demand response, which is equal to the total gas load minus the largest purchased natural gas in each period.

In addition, due to the difficulty of gas supply, CHP's ability to support power supply is extremely weak, so they all operate at a very low level of power generation. At this time, since the power load cannot get assistance from the natural gas system, almost all loads have to participate in demand management. This demand management scheme that does not consider natural gas support is called full-emergency scheme in this paper, and it is necessary to use all demand response capabilities.

It is worth mentioning that regarding the supply of heat load, in the electrical resilience mode, although pumps cannot provide heat due to insufficient power supply, gas can be used to provide sufficient heat through the boilers without causing heat load loss. In the integrated resilience mode, due to the failure of No. 8, 11, and 15 pipelines, No. 1 boiler cannot obtain gas for heating, and No. 2 boiler is used for heating, supplementally. Although part of the heat load is involved in demand response, it is not the focus of this paper and will not be analyzed in detail.

### *4.4. The Demand Management Strategy*

From the above, the demand management strategy includes a semi-emergency demand management scheme and a full-emergency demand management scheme, according to the type of natural disaster. Moreover, the demand management strategy is based on the importance of lines and pipelines. As shown in Tables 1 and 2, the failure of each transmission line and gas pipeline corresponds to a certain economic subsidy to the load participating in the demand response. According to the semi-emergency scheme and the full-emergency scheme, the upper and lower cost limits of the corresponding demand management for each failure transmission line and gas pipeline can be obtained, as shown in Figures 9 and 10.

As shown in Figures 9 and 10, according to the faults caused by natural disasters generated by non-sequential Monte-Carlo simulation, the economic subsidies for demand management caused by the failure of various transmission lines and gas pipelines fluctuate relatively widely. Therefore, the demand management strategy proposed in this paper can obtain the optimal demand management cost according to the failure situation, which can ensure that IEDS can not only improve the resilience, but also have a lower operating cost.

**Figure 9.** The economic subsidy of power demand response.

**Figure 10.** The economic subsidy of gas demand response.

### **5. Conclusions**

IEDS is an efficient, clean and sustainable energy system. However, increasingly frequent natural disasters pose a severe threat to the security of IDES. Hence, this paper studies demand management for IEDS to improve resilience against extreme events. The demand management strategy includes a semi-emergency demand management scheme and a full-emergency demand management scheme, according to the resilience mode. The strategy generation process is divided into three stages. In short, combined with component generic fragility curves, through the non-sequential Monte-Carlo simulation method and the scene reduction algorithm, the main scene of the failure caused by the natural disaster is simulated, and then according to the importance ranking of transmission line and gas pipeline the demand management economic subsidy costs under different failure situations are obtained. The main conclusions are summarized as follows:

(1) The proposed scheduling model of IEDS co-optimizes the electricity, natural gas and thermal energy in normal and failure operation state, which promotes efficient and sustainable energy consumption.

(2) This paper proposes a demand management strategy, which innovatively expands the traditional demand management in a power system. The demand management of the IEDS includes two parts: power demand management and gas demand management. This paper analyzes the economic subsidies for demand response in these two parts, respectively.

(3) The demand management strategy proposed in this paper can increase the resilience of IEDS and reduce the economic subsidy cost of demand management.

**Author Contributions:** Conceptualization, Y.X. and S.C.; methodology, S.T. and F.G.; software, Y.X.; validation, Y.X., S.C., S.T. and F.G.; formal analysis, Y.X.; investigation, Y.X.; resources, S.C.; data curation, S.T.; writing—original draft preparation, Y.X.; writing—review and editing, Y.X. and F.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by Science and Technology Project of State Grid "Research on Largescale Interactive Response Technology of Demand-Side Flexible Resources with High Penetration of Renewable Energy" (No. 5100-202114296A-0-0-00).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Abbreviations**


*f load n*,*t*,*s* Natural gas load at node *n* in Scenarios *s* at hour *t*. *hload <sup>t</sup>*,*<sup>s</sup>* Heat load in Scenarios *s* at hour *t*. *rilsjle*/*xilsjle* Resistance/reactance of overhead transmission line *ij*. *Vbase* Base value of voltage. *pl ilsjle*/*q<sup>l</sup> ilsjle* Maximum active/reactive power of overhead transmission line *ij*. *<sup>p</sup>pur*/*qpur* Maximum active/reactive power purchased from IEDS. *Vb*/*Vb* Maximum/minimum voltage at bus *b*. *pw iw*,*t*,*s* Forecasted output of wind generator *i <sup>w</sup>* in Scenarios *s* at hour *t*. *φiw* Power factor of wind generator *i w*. *Kigsjge* Gas flow constant of gas pipeline *ij*. *<sup>β</sup>*(*zigs*,*t*,*s*, *zjge*,*t*,*s*) Direction coefficient of gas flow in gas pipeline *ij*. *f g <sup>i</sup>gsjge* Maximum natural gas flow of gas pipeline *ij*. *f pur* Maximum natural gas purchased from IEDS. *zn*/*zn* Maximum/minimum pressure at node *n*. *M* Sufficiently big number *xigsjge*,*<sup>k</sup>* Abscissa breakpoints of domain *yigsjge*,*<sup>k</sup>* Ordinate value corresponding to *xigsjge*,*<sup>k</sup> εigsjge*/*KP* Size/number of interval. *pCHP iCHP* /*pCHP iCHP* Maximum/minimum active power output of CHP units *i CHP*. *qCHP iCHP* /*qCHP iCHP* Maximum/minimum reactive power output of CHP units *i CHP*. *ηCHP*,*ele iCHP* /*ηCHP*,*heat iCHP* Electrical/heat efficiency of CHP units *<sup>i</sup> CHP*. *pEHP iEHP* / *f boi iboi* Maximum power/gas consumption of EHP *i EHP*/boiler *i boi*. *ηEHP iEHP* /*ηboi iboi* Heat efficiency of EHP*i EHP*/boiler *i boi*. **Variables** *ppur <sup>i</sup>pe*,*t*,*s*/*q pur <sup>i</sup>pe*,*t*,*s*/ *f pur ipg*,*t*,*s* Active power/reactive power/natural gas purchased from IEDS in scenarios *s* at hour *t*. *pdemand <sup>b</sup>*,*t*,*<sup>s</sup>* /*qdemand <sup>b</sup>*,*t*,*<sup>s</sup>* / *<sup>f</sup> demand <sup>n</sup>*,*t*,*<sup>s</sup>* /*hdemand t*,*s* Active power/reactive power/natural gas/heat load demand response at bus *b*/at node *n* in Scenarios *s* at hour *t*. *pw*,*cut iw*,*t*,*s* Power curtailment of wind generator *i <sup>w</sup>* in Scenarios *s* at hour *t*. *pw iw*,*t*,*s*/ *<sup>q</sup><sup>w</sup> iw*,*t*,*s* Scheduled active/reactive power output of wind generator*i <sup>w</sup>* in Scenarios *s* at hour *t*. *qw iw*,*t*,*s* Scheduled reactive power output of wind generator *i <sup>w</sup>* in Scenarios *s* at hour *t*. *pCHP iCHP*,*t*,*s* /*qCHP iCHP*,*t*,*s* Active power/reactive power of CHP unit *i CHP* in Scenarios *s* at hour *t*. *pEHP iEHP*,*t*,*s* Power consumption of electrical heat pump *i EHP* in Scenarios *s* at hour *t*. *pl ilsjle*,*t*,*s* /*q<sup>l</sup> ilsjle*,*t*,*s* Active power/reactive power of overhead transmission line *ij* in Scenarios *s* at hour *t*. *Vils*,*t*,*s*/*Vjle*,*t*,*<sup>s</sup>* Voltage at overhead transmission line start/end point in Scenarios *s* at hour *t*. *Vb*,*t*,*<sup>s</sup>* Voltage at bus *b* in Scenarios *s* at hour *t*. *f CHP iCHP*,*t*,*s* Gas consumption of CHP unit *i CHP* in Scenarios *s* at hour *t*.


### **References**

