Enhancing Integrated Energy Distribution System Resilience through a Hierarchical Management Strategy in District Multi-Energy Systems

Abstract

The requirement for energy sustainability drives the development of integrated energy distribution systems (IEDSs). In this paper, considering the coordination of district multi-energy systems (DMESs), a hierarchical management strategy is proposed to enhance IEDS resilience. The proposed strategy is divided into three modes: the normal operation mode, the preventive operation mode and the resilient operation mode. In the normal operation mode, the objective of DEMSs is to minimize the operation costs. In the preventive operation mode, the objective of DEMSs is to maximize the stored energy for mitigating outage. The resilient operation mode consists of two stages. DMESs schedule their available resources and broadcast excess generation capacities or unserved loads to neighboring DMESs through the cyber communication network in the first stage. In the second stage, DMESs interchange electricity and natural gas with each other through the physical common bus for global optimization. A consensus algorithm was applied to determine the allocated proportions of exported or imported electricity and natural gas for each DMES in a distributed way. An IEDS including five DMESs was used as a test system. The results of the case studies demonstrate the effectiveness of the proposed hierarchical management strategy and algorithm.

1. Introduction

With the development of integrated energy distribution system (IEDSs) [1,2], the interaction among electricity distribution systems, natural gas distribution systems and district heating systems is increasingly intensified by the continuing penetration of gas-fired combined heat and power (CHP) units [3]. Moreover, based on the energy hub (EH) concept [4], district multi-energy systems (DMESs) consisting of micro-grids (MGs) and district gas-heating networks can coordinately schedule the interrelated electricity, natural gas and heat. However, frequent natural disasters have been causing severe threats to the security of IEDSs. For example, in 2012, after Hurricane Sandy hit the eastern shore of the United States, approximately 7.5 million customers were reportedly without power [5] and the cost of outages estimated ranged from $27 billion to $52 billion [6]. In this regard, it is necessary to develop methods to enhance IEDS resilience. Furthermore, IEDS resilience and sustainability are closely related [7].

However, compared with IEDS resilience, many studies have been made in the literature on power system resilience [8,9,10,11,12,13,14,15]. Power system resilience can be defined as the ability to anticipate, withstand, absorb and recover rapidly from disruptions caused by extreme natural catastrophes such as earthquakes, hurricanes and flooding [16].

To enhance power system resilience, the strategies in previous research could be roughly categorized into the following aspects: placing flexible back-up resources [8,9], improving response ability [10,11], altering network topology [12,13,14] and hardening distribution lines [15]. The study in [8] proposes a two-stage dispatch framework consisting of the pre-positioning and real-time allocation of mobile emergency generators to restore critical loads in electricity distribution systems. In addition, the study in [9] proposes a stochastic mixed-integer second-order conic program for the optimal investments and transportation routes of mobile energy storages to avoid load shedding by disasters. For another, stationary devices can also improve resilience after extreme events. In [10], a method is presented for the optimal switch placement in distribution systems, taking into account resilience ability against major faults caused by natural disasters. In [11], a Markov model is proposed to construct a sequentially proactive operation strategy with generation re-dispatch to reduce potential power loss. Besides, an integrated framework is proposed in [12] to schedule topology switching and adjustable load shedding for resilience response. To utilize existing resources more efficiently, the study in [13] develops a self-healing strategy by sectionalizing a distribution system into self-adequate MGs with distributed generators (DGs). The study in [14] proposes an operational approach for forming multiple MGs energized by DG via controlling the automatic switch devices so as to provide power for interrupted loads. To withstand disruption from natural disasters, a tri-level robust optimization is proposed to harden distribution lines in [15]. Furthermore, MGs can cover small geographical areas and maintain critical loads in the inland mode, so they can mitigate load reductions during transmission system outages. For example, after the Great East Japan Earthquake in March 2011, the Sendai MG continuously supplied power to some critical loads within small area until grid restoration [17]. The study in [18] presents the mesh grid approach for calculating proposed resilience indices, which demonstrate that island MDs can improve the resiliency of power grids. Additionally, the networked MGs can further improve power grid resilience by supporting and interchanging electricity with each other [19]. In particular, a hierarchical outage management framework is proposed in [20] for the networked MGs to coordinate available resources.

Actually, with the enhancement of coupling among electricity, natural gas and heat, it is not suitable that interactive support and exchange only consider electricity in several DMESs. Hence, it is significant and urgent to research IEDS resilience. IEDS resilience, which is an extension of the author’s review in [21], represents the ability of IEDS to withstand low probability high impact events in efficient means while ensuring the least outages of electricity, natural gas and heat, and further enabling quick restoration to the normal operation mode. However, the literature contains little research on IEDS resilience [22,23]. The study in [22] proposes a two-stage robust model for the integrated planning of electricity and natural gas system to improve power grid resilience via the replacement of some power lines by a natural gas transportation system. Moreover, an integrated electricity–heat–gas network model is presented for a holistic evaluation of energy system resilience in [23].

From the above, IEDS resilience has not been sufficiently studied yet. In particular, the restoration approach, which includes multiple energies interchange and supporting each other in DMESs to avoid interrupting power and gas supply during electricity distribution system and natural gas distribution system disruption by disaster, has not yet been investigated. Therefore, motivated by the aforementioned facts, this paper proposes a novel hierarchical management strategy for the coordinated operation and management of DMESs to enhance IEDS resilience. This paper divides the operation of DMESs into three modes: the normal operation mode, the preventive operation mode and the resilient operation mode. In the normal operation mode, the objective of DMESs is to minimize the operation costs. After receiving a disaster alert, DMESs enter the preventive operation mode in which the objective of DEMSs is to maximize the stored energy, so that adequate energy in batteries, gasholders and thermal storages of DEMSs can support electrical loads, gas loads and heat loads after IEDS outages. When natural disaster disrupts IEDS, the power and gas purchased from IEDS will be limited or even interrupted severely. In the resilient operation mode, the hierarchical management strategy is a two-stage optimization model, which is formulated as a mixed-integer linear programming (MILP), for scheduling electricity and natural gas between the networked DMESs to enhance IEDS resilience. A consensus algorithm [24] is applied to determine the allocated proportions of exported or imported electricity and natural gas for each DMES in a distributed way and the MILP can be solved by commercial solvers. The prominent features of the proposed strategy are demonstrated by case studies.

The major contributions of this paper are summarized as follows:

(1)

This paper co-optimizes the electricity, natural gas and thermal energy for three operation modes in DMES, which promotes the comprehensive and efficient utilization of different energy in normal and destroyed situations.

(2)

The hierarchical management strategy in DMES is proposed, which has enhancement for IEDS resilience to withstand the extreme events, such as natural disaster.

(3)

The consensus algorithm is applied to obtain the most suitable energy exchange with the least visible information on DMESs, which can not only calculate the allocated proportions accurately in a distributed way but also protect the privacy of DMESs effectively.

The rest of this paper is organized as follows. The concept of the networked DMESs is introduced and the general framework of the hierarchical management strategy is discussed in Section 2. Section 3 proposes mathematical formulations of the optimal normal, preventive and resilient operation and develops a consensus algorithm to solve this problem in a distributed way. In Section 4, numerical results are provided to validate the proposed strategy. Finally, Section 5 concludes the paper.

2. Hierarchical Management Strategy

2.1. Networked District Multi-Energy Systems

The DMESs, consisting of MGs and district gas-heating networks, are small-scale integrated energy systems covering small geographical areas with DGs, energy storages (ESs), energy coupling devices and various energy loads, which are shown in Figure 1.

As shown in Figure 1, ESs include battery, gasholder and thermal storage and energy coupling devices include CHP units, power to gas, electrical heat pumps (EHPs) [25] and boilers. In addition to the mentioned advantages, when the power and gas in the main grid are not available or the IEDS is faulty, customers in a DMES could benefit from the power and gas supplied from stored resources. Hence, DMESs could be employed to enhance IEDS resilience by reducing the possibility of load shedding.

However, each DMES can only support unserved loads by scheduling their own resources. Since each DMES is controlled independently, it is not an optimal strategy for multiple DMESs. Therefore, the networked DMESs, which is further development and application of the concept of DMES, need to be studied, which can coordinate energy management among connective multiple DMESs. Figure 2 illustrates such the networked DMESs model.

As shown in Figure 2, the networked DMESs consist of multiple DMESs connected through the physical common electricity bus and gas bus. Moreover, for information exchange, the DMESs are also connected through the cyber communication network. The networked DMESs can support each other with the remainder power and gas through the common buses to supply more loads after a natural disaster.

In IEDS scheduling management, each DMES is managed by an integrated energy serving company, which is the lower-layer controller (LLC). However, since there is no central controller in the upper-layer networked DMES, each DMES can only exchange information between neighboring DMESs. In the model, two conditions need to be met. First, each LLC has ownership of a DMES. This means that each DMES is autonomously managed by its own LLC, so the hierarchical management strategy should fulfil the autonomy of DMESs as much as possible. Second, the privacy in each DMES needs to be protected, so the energy and transaction information on DMESs are visible minimally.

2.2. General Framework of the Hierarchical Management Strategy

As noted earlier, the operation of DMESs can be divided into three modes, as depicted in Figure 3. In the normal operation mode, the LLC manages DGs, ESs and energy coupling devices, purchases energy from IEDSs and sells energy to consumers without energy exchange with other DMESs. Hence, the network lines and common bus are idle. After a disaster alert is announced, the proactive scheduling objective of DMESs is to maximize the stored energy, so the LLCs purchase energy to fill ESs as soon as possible in the preventive operation mode. When an IEDS is damaged by disaster, the DMESs cannot purchase adequate electricity and natural gas from IEDSs, so they have to enter the resilient operation mode.

In the resilient operation mode, the hierarchical management strategy proposed in this paper is comprised of two stages. In the first stage, the LLC develops a scheduling plan for the next timeslot, which includes DGs and ESs output and load consumptions in its own DMES. Then, the LLC broadcasts surplus generation capacities or unserved loads to neighboring DMESs through the cyber communication network. In the second stage, the scheduling plans are updated by exchange energy. According to the announced amount of excess and deficit electricity and natural gas, the allocated proportions of exported or imported electricity and natural gas, which can be used to determine the amount of exchange energy between DMESs, are calculated for each DMES by the proposed consensus algorithm. With the updated scheduling plan implemented, the excess DMESs provide the surplus energy to the deficit DMESs through the common bus so as to reduce load shedding when purchased energy is limited. To put it more clearly, the networked DMESs share the available generation resources and storage capacities to supply more loads as a resilient response, which is a self-healing process to enhance IEDS resilience.

It is worth highlighting that the scheduling plan is updated in two steps. In detail, after disaster disruption, the purchased electricity and natural gas from IEDSs are limited. Actually, besides the stored energy in battery and thermal storage, the supply source of electrical loads, gas loads and heat loads is the stored natural gas in gasholders. Therefore, in this paper, according to load shedding cost in a short time, the electrical loads should be supplied as much as possible, followed by the gas loads, and if necessary, the heat loads can be curtailed first. In the first step, the gas consumption and output of CHP units and the output of the battery and gasholder are updated, so that the electrical load curtailment is mitigated preferentially. In the second step, the output of the gasholder is updated again for the mitigation of gas load curtailment. All steps of the proposed hierarchical management strategy are briefly depicted in Figure 4.

Moreover, this paper has two assumptions. First, it is assumed that the common buses and the cyber network are strongly connected and no island exists in the networked DMESs after extreme disaster. Second, in order to cope with potential natural disasters, all the LLCs in the networked DMESs have signed the control protocols beforehand, which is based on the strategy proposed in this paper. What is more, the scheduling horizon is one hour. The purchased energy will be available after IEDS eliminates the faults caused by the disaster. In the case of an uncertain IEDS repair time, the networked DMESs should provide the existing loads to customers reliably. After a scheduling horizon, the proposed strategy will be implemented in the next hour, and the procedure will continue until IEDS restoration. The model formulation of this procedure is discussed in detail in the following section.

3. Methodology

3.1. Operation Model Formulation

The operation of DMESs has three modes, including the normal, preventive and resilient operation mode, and different operation modes have different objectives.

Objective function:

$$\begin{array}{l}\min \{{\displaystyle \sum _{t=1}^{T_1}({\displaystyle \sum _{n=1}^{N^{DMES}}(\underset{\mathrm{i}}{\underbrace{{\lambda }_t^e{p}_{n,t}^{pur}}}+\underset{\mathrm{ii}}{\underbrace{{\lambda }_t^g{f}_{n,t}^{pur}}}+\underset{\mathrm{iii}}{\underbrace{{\pi }_n^e{p}_{n,t}^{shed}}}+\underset{\mathrm{iv}}{\underbrace{{\pi }_n^g{f}_{n,t}^{shed}}}+\underset{\mathrm{v}}{\underbrace{{\pi }_n^h{q}_{n,t}^{shed}}}}}\\ +\underset{\mathrm{vi}}{\underbrace{{\displaystyle \sum _{i^{bat}=1}^{N_n^{bat}}({\theta }_{i^{bat},n}^{bat,ch}{p}_{i^{bat},n,t}^{ch}+{\theta }_{i^{bat},n}^{bat,dch}{p}_{i^{bat},n,t}^{dch})}}}+\underset{\mathrm{vii}}{\underbrace{{\displaystyle \sum _{i^{gas}=1}^{N_n^{gas}}({\theta }_{i^{gas},n}^{gas,ch}{f}_{i^{gas},n,t}^{ch}+{\theta }_{i^{gas},n}^{gas,dch}{f}_{i^{gas},n,t}^{dch})}}}\\ +\underset{\mathrm{viii}}{\underbrace{{\displaystyle \sum _{i^{heat}=1}^{N_n^{heat}}({\theta }_{i^{heat},n}^{heat,ch}{q}_{i^{heat},n,t}^{ch}+{\theta }_{i^{heat},n}^{heat,dch}{q}_{i^{heat},n,t}^{dch})}}}+\underset{\mathrm{ix}}{\underbrace{{\displaystyle \sum _{i^{wind}=1}^{N_n^{wind}}{\kappa }_{i^{wind}}^{wind}{p}_{i^{wind},n,t}^{wind,cut}}}}+\underset{\mathrm{x}}{\underbrace{{\displaystyle \sum _{i^{solar}=1}^{N_n^{solar}}{\kappa }_{i^{solar}}^{solar}{p}_{i^{solar},n,t}^{solar,cut}}}}))\}\end{array}$$

(1)

$$\begin{array}{l}\min \{{\displaystyle \sum _{t=T_1+1}^{T_1+T_2}({\displaystyle \sum _{n=1}^{N^{DMES}}(\underset{\mathrm{i}}{\underbrace{{\lambda }_t^e{p}_{n,t}^{pur}}}+\underset{\mathrm{ii}}{\underbrace{{\lambda }_t^g{f}_{n,t}^{pur}}}+\underset{\mathrm{iii}}{\underbrace{{\pi }_n^e{p}_{n,t}^{shed}}}+\underset{\mathrm{iv}}{\underbrace{{\pi }_n^g{f}_{n,t}^{shed}}}+\underset{\mathrm{v}}{\underbrace{{\pi }_n^h{q}_{n,t}^{shed}}}}}\\ +\underset{\mathrm{vi}}{\underbrace{{\displaystyle \sum _{i^{bat}=1}^{N_n^{bat}}({\theta }_{i^{bat},n}^{bat,ch}{p}_{i^{bat},n,t}^{ch}+{\theta }_{i^{bat},n}^{bat,dch}{p}_{i^{bat},n,t}^{dch})}}}+\underset{\mathrm{vii}}{\underbrace{{\displaystyle \sum _{i^{gas}=1}^{N_n^{gas}}({\theta }_{i^{gas},n}^{gas,ch}{f}_{i^{gas},n,t}^{ch}+{\theta }_{i^{gas},n}^{gas,dch}{f}_{i^{gas},n,t}^{dch})}}}\\ +\underset{\mathrm{viii}}{\underbrace{{\displaystyle \sum _{i^{heat}=1}^{N_n^{heat}}({\theta }_{i^{heat},n}^{heat,ch}{q}_{i^{heat},n,t}^{ch}+{\theta }_{i^{heat},n}^{heat,dch}{q}_{i^{heat},n,t}^{dch})}}}+\underset{\mathrm{ix}}{\underbrace{{\displaystyle \sum _{i^{wind}=1}^{N_n^{wind}}{\kappa }_{i^{wind}}^{wind}{p}_{i^{wind},n,t}^{wind,cut}}}}+\underset{\mathrm{x}}{\underbrace{{\displaystyle \sum _{i^{solar}=1}^{N_n^{solar}}{\kappa }_{i^{solar}}^{solar}{p}_{i^{solar},n,t}^{solar,cut}}}}\\ +\underset{\mathrm{xi}}{\underbrace{{\displaystyle \sum _{i^{bat}=1}^{N_n^{bat}}{\omega }_{i^{bat}}^{bat}({\overline{E}}_{i^{bat},n}^{bat}-E_{i^{bat},n,t}^{bat})}}}+\underset{\mathrm{xii}}{\underbrace{{\displaystyle \sum _{i^{gas}=1}^{N_n^{gas}}{\omega }_{i^{gas}}^{gas}({\overline{E}}_{i^{gas},n}^{gas}-E_{i^{gas},n,t}^{gas})}}}+\underset{\mathrm{xiii}}{\underbrace{{\displaystyle \sum _{i^{heat}=1}^{N_n^{heat}}{\omega }_{i^{heat}}^{heat}({\overline{E}}_{i^{heat},n}^{heat}-E_{i^{heat},n,t}^{heat})}}}))\}\end{array}$$

(2)

$$\begin{array}{l}\min \{{\displaystyle \sum _{t=T_1+T_2+1}^{T_1+T_2+T_3}({\displaystyle \sum _{n=1}^{N^{DMES}}(\underset{\mathrm{iii}}{\underbrace{{\pi }_n^e{p}_{n,t}^{shed}}}+\underset{\mathrm{iv}}{\underbrace{{\pi }_n^g{f}_{n,t}^{shed}}}+\underset{\mathrm{v}}{\underbrace{{\pi }_n^h{q}_{n,t}^{shed}}}}}+\underset{\mathrm{vi}}{\underbrace{{\displaystyle \sum _{i^{bat}=1}^{N_n^{bat}}({\theta }_{i^{bat},n}^{bat,ch}{p}_{i^{bat},n,t}^{ch}+{\theta }_{i^{bat},n}^{bat,dch}{p}_{i^{bat},n,t}^{dch})}}}\\ +\underset{\mathrm{vii}}{\underbrace{{\displaystyle \sum _{i^{gas}=1}^{N_n^{gas}}({\theta }_{i^{gas},n}^{gas,ch}{f}_{i^{gas},n,t}^{ch}+{\theta }_{i^{gas},n}^{gas,dch}{f}_{i^{gas},n,t}^{dch})}}}+\underset{\mathrm{viii}}{\underbrace{{\displaystyle \sum _{i^{heat}=1}^{N_n^{heat}}({\theta }_{i^{heat},n}^{heat,ch}{q}_{i^{heat},n,t}^{ch}+{\theta }_{i^{heat},n}^{heat,dch}{q}_{i^{heat},n,t}^{dch})}}}\\ +\underset{\mathrm{ix}}{\underbrace{{\displaystyle \sum _{i^{wind}=1}^{N_n^{wind}}{\kappa }_{i^{wind}}^{wind}{p}_{i^{wind},n,t}^{wind,cut}}}}+\underset{\mathrm{x}}{\underbrace{{\displaystyle \sum _{i^{solar}=1}^{N_n^{solar}}{\kappa }_{i^{solar}}^{solar}{p}_{i^{solar},n,t}^{solar,cut}}}}))\}\end{array}$$

(3)

In the normal operation mode, the main objectives are to minimize the operation costs under the supply–demand balance of each energy in DMESs. The goal of the LLCs is to minimize the total cost of DMESs, as shown in (1). The total cost consists of the following terms: (i) electricity purchasing cost, (ii) natural gas purchasing cost, (iii) electrical load shedding cost, (iv) gas load shedding cost, (v) heat load shedding cost, (vi) charging and discharging cost of batteries, (vii) charging and discharging cost of gasholders, (viii) charging and discharging cost of thermal storages, (ix) penalty cost of wind power curtailment, and (x) penalty cost of solar power curtailment.

In the preventive operation mode, the objectives are not only to minimize the operation costs, but also to maximize the energy stored, so that DMESs can supply the load as much as possible after the energy purchased from IEDSs is limited. To facilitate modeling, the idle capacity of energy storages is calculated as the penalty cost: (xi) penalty cost of idle batteries, (xii) penalty cost of idle gasholders, and (xii) penalty cost of idle thermal storages.

In the resilient operation mode, after disaster disruption, the IEDS is destroyed and the energy purchased is limited. Hence, the main objectives are to minimize the load shedding. Each DMES supplies its own loads preferentially and broadcasts excess and deficit energy to adjacent DMESs. And then, the amount of energy exchanged is determined by the consensus algorithm, which will be introduced in the Section 3.2.

Furthermore, the constraints of the operation include energy balance, the operation status of equipment, the energy purchased, renewable energy output and load shedding.

Energy balance constraints:

$$p_{n,t}^{pur}+{\displaystyle \sum _{i^{wind}=1}^{N_n^{wind}}{p}_{i^{wind},n,t}^{wind}}+{\displaystyle \sum _{i^{solar}=1}^{N_n^{solar}}{p}_{i^{solar},n,t}^{solar}}+{\displaystyle \sum _{i^{CHP}=1}^{N_n^{CHP}}{p}_{i^{CHP},n,t}^{CHP}}+{\displaystyle \sum _{i^{bat}=1}^{N_n^{bat}}(p_{i^{bat},n,t}^{dch}-p_{i^{bat},n,t}^{ch})}+p_{n,t}^{shed}=p_{n,t}^{load}+{\displaystyle \sum _{i^{EHP}=1}^{N_n^{EHP}}{p}_{i^{EHP},n,t}^{EHP}} \forall n, \forall t$$

(4)

$$f_{n,t}^{pur}+{\displaystyle \sum _{i^{gas}=1}^{N_n^{gas}}(f_{i^{gas},n,t}^{dch}-}{f}_{i^{gas},n,t}^{ch})+f_{n,t}^{shed}=f_{n,t}^{load}+{\displaystyle \sum _{i^{CHP}=1}^{N_n^{CHP}}{f}_{i^{CHP},n,t}^{CHP}}+{\displaystyle \sum _{i^{boi}=1}^{N_n^{boi}}{f}_{i^{boi},n,t}^{boi}} \forall n, \forall t$$

(5)

$${\displaystyle \sum _{i^{CHP}=1}^{N_n^{CHP}}{q}_{i^{CHP},n,t}^{CHP}}+{\displaystyle \sum _{i^{boi}=1}^{N_n^{boi}}{q}_{i^{boi},n,t}^{boi}}+{\displaystyle \sum _{i^{EHP}=1}^{N_n^{EHP}}{q}_{i^{EHP},n,t}^{EHP}}+{\displaystyle \sum _{i^{heat}=1}^{N_n^{heat}}(q_{i^{heat},n,t}^{dch}-q_{i^{heat},n,t}^{ch})}+q_{n,t}^{shed}=q_{n,t}^{load} \forall n, \forall t$$

(6)

Constraints (4)–(6) are power, natural gas and heat balance constraints, respectively. For each DMES and at each hour, (4) the sum of the total generated power by electrical power purchased, wind generators, photovoltaic plants, CHP units, charging/discharging power of batteries and curtailed electrical load must be equal to the sum of the total electrical load by conventional electrical load and power consumption of EHPs; (5) the sum of the total generated gas by gas purchased, charging/discharging gas of gasholders and curtailed gas load must be equal to the sum of the total gas load by conventional gas load and the gas consumption of CHP units and boilers; (6) the sum of the total generated heat by CHP units, boilers, EHPs, charging/discharging heat of thermal storages and curtailed heat load must be equal to the sum of the total heat load by conventional gas load.

The constraints of the operation status of the equipment include not only coupling equipment between different energy systems, such as CHP unit constraints, EHP constraints and boiler constraints, but also energy storage equipment, such as battery constraints, gasholder constraints and thermal storage constraints.

CHP unit constraints:

$$p_{i^{CHP},n,t}^{CHP}={\tau }_{i^{CHP},n}{\eta }_{i^{CHP},n}^{CHP,ele}{f}_{i^{CHP},n,t}^{CHP} \forall i^{CHP}, \forall n, \forall t$$

(7)

$$q_{i^{CHP},n,t}^{CHP}=(1-{\tau }_{i^{CHP},n}){\eta }_{i^{CHP},n}^{CHP,heat}{f}_{i^{CHP},n,t}^{CHP} \forall i^{CHP}, \forall n, \forall t$$

(8)

$$0\le p_{i^{CHP},n,t}^{CHP}\le {\overline{p}}_{i^{CHP},n}^{CHP} \forall i^{CHP}, \forall n, \forall t$$

(9)

$$\Delta {\underset{¯}{p}}_{i^{CHP},n}^{CHP}\le p_{i^{CHP},n,t+1}^{CHP}-p_{i^{CHP},n,t}^{CHP}\le \Delta {\overline{p}}_{i^{CHP},n}^{CHP} \forall i^{CHP}, \forall n, \forall t$$

(10)

Constraints (7) and (8) show the energy conversion among the electricity generation, heat generation and natural gas consumption of CHP units. Constraints (9) limit the power output of CHP units. Constraints (10) limit the ramp-up and ramp-down of CHP units.

EHP constraints:

$$q_{i^{EHP},n,t}^{EHP}={\eta }_{i^{EHP},n}^{EHP}{p}_{i^{EHP},n,t}^{EHP} \forall i^{EHP}, \forall n, \forall t$$

(11)

$$0\le q_{i^{EHP},n,t}^{EHP}\le {\overline{q}}_{i^{EHP},n}^{EHP} \forall i^{EHP}, \forall n, \forall t$$

(12)

$$\Delta {\underset{¯}{q}}_{i^{EHP},n}^{EHP}\le q_{i^{EHP},n,t+1}^{EHP}-q_{i^{EHP},n,t}^{EHP}\le \Delta {\overline{q}}_{i^{EHP},n}^{EHP} \forall i^{EHP}, \forall n, \forall t$$

(13)

Constraints (11) show the energy conversion between the heat generation and electrical power consumption of EHP. Constraints (12) limit the heat output of EHP. Constraints (13) limit the ramp-up and ramp-down of EHP.

Boiler constraints:

$$q_{i^{boi},n,t}^{boi}={\eta }_{i^{boi},n}^{boi}{f}_{i^{boi},n,t}^{boi} \forall i^{boi}, \forall n, \forall t$$

(14)

$$0\le q_{i^{boi},n,t}^{boi}\le {\overline{q}}_{i^{boi},n}^{boi} \forall i^{boi}, \forall n, \forall t$$

(15)

$$\Delta {\underset{¯}{q}}_{i^{boi},n}^{boi}\le q_{i^{boi},n,t+1}^{boi}-q_{i^{boi},n,t}^{boi}\le \Delta {\overline{q}}_{i^{boi},n}^{boi} \forall i^{boi}, \forall n, \forall t$$

(16)

Constraints (14) show the energy conversion between the heat generation and natural gas consumption of the boiler. Constraints (12) limit the heat output of the boiler. Constraints (13) limit the ramp-up and ramp-down of the boiler.

Battery constraints:

$$E_{i^{bat},n,t+1}^{bat}=E_{i^{bat},n,t}^{bat}+p_{i^{bat},n,t}^{ch}{\eta }_{i^{bat},n}^{bat,ch}\Delta t-\frac{p_{i^{bat},n,t}^{dch}\Delta t}{{\eta }_{i^{bat},n}^{bat,dch}} \forall i^{bat}, \forall n, \forall t$$

(17)

$$0\le E_{i^{bat},n,t}^{bat}\le {\overline{E}}_{i^{bat},n}^{bat} \forall i^{bat}, \forall n, \forall t$$

(18)

$$0\le p_{i^{bat},n,t}^{ch}\le {\chi }_{i^{bat},n,t}^{bat}{\overline{p}}_{i^{bat},n}^{ch} \forall i^{bat}, \forall n, \forall t$$

(19)

$$0\le p_{i^{bat},n,t}^{dch}\le (1-{\chi }_{i^{bat},n,t}^{bat}){\overline{p}}_{i^{bat},n}^{dch} \forall i^{bat}, \forall n, \forall t$$

(20)

Constraints (17) show the relationship of electricity stored and charging/discharging power. Constraints (18) show the allowable limit of battery capacity. Allowable charging/discharging power limits are imposed through constraints (19) and (20).

Gasholder constraints:

$$E_{i^{gas},n,t+1}^{gas}=E_{i^{gsa},n,t}^{gas}+f_{i^{gas},n,t}^{ch}{\eta }_{i^{gas},n}^{gas,ch}\Delta t-\frac{f_{i^{gas},n,t}^{dch}\Delta t}{{\eta }_{i^{gas},n}^{gas,dch}} \forall i^{gas}, \forall n, \forall t$$

(21)

$$0\le E_{i^{gas},n,t}^{gas}\le {\overline{E}}_{i^{gas},n}^{gas} \forall i^{gas}, \forall n, \forall t$$

(22)

$$0\le f_{i^{gas},n,t}^{ch}\le {\chi }_{i^{gas},n,t}^{gas}{\overline{f}}_{i^{gas},n}^{ch} \forall i^{gas}, \forall n, \forall t$$

(23)

$$0\le f_{i^{gas},n,t}^{dch}\le (1-{\chi }_{i^{gas},n,t}^{gas}){\overline{f}}_{i^{gas},n}^{dch} \forall i^{gas}, \forall n, \forall t$$

(24)

Constraints (21) show the relationship of natural gas stored and charging/discharging gas. Constraints (22) show the allowable limit of gasholder capacity. Constraints (23) and (24) limit charging/discharging gas.

Thermal storage constraints:

$$E_{i^{heat},n,t+1}^{heat}=E_{i^{heat},n,t}^{heat}+q_{i^{heat},n,t}^{ch}{\eta }_{i^{heat},n}^{heat,ch}\Delta t-\frac{q_{i^{heat},n,t}^{dch}\Delta t}{{\eta }_{i^{heat},n}^{heat,dch}} \forall i^{heat}, \forall n, \forall t$$

(25)

$$0\le E_{i^{heat},n,t}^{heat}\le {\overline{E}}_{i^{heat},n}^{heat} \forall i^{heat}, \forall n, \forall t$$

(26)

$$0\le q_{i^{heat},n,t}^{ch}\le {\chi }_{i^{heat},n,t}^{heat}{\overline{q}}_{i^{heat},n}^{ch} \forall i^{heat}, \forall n, \forall t$$

(27)

$$0\le q_{i^{heat},n,t}^{dch}\le (1-{\chi }_{i^{heat},n,t}^{heat}){\overline{q}}_{i^{heat},n}^{dch} \forall i^{heat}, \forall n, \forall t$$

(28)

Constraints (25) show the relationship of heat stored and charging/discharging heat. Constraints (26) show the allowable limit of thermal storage capacity. Constraints (27) and (28) limit charging/discharging heat. What is more, it is worth emphasizing that the binary variables ${\chi }_{i^{bat},n,t}^{bat}$, ${\chi }_{i^{gas},n,t}^{gas}$ and ${\chi }_{i^{heat},n,t}^{heat}$ are guaranteed to avoid the simultaneous charging and discharging of energy storage equipment.

Since heat is not suitable for long-distance transmission, the energy purchased off DMESs is mainly electrical power and natural gas. Hence, heat is supplied by CHP units, EHPs and boilers in DMESs.

Energy purchased constraints:

$$0\le p_{n,t}^{pur}\le {\overline{p}}_n^{pur} \forall n, \forall t$$

(29)

$$0\le f_{n,t}^{pur}\le {\overline{f}}_n^{pur} \forall n, \forall t$$

(30)

Constraints (29) and (30) limit electrical power and natural gas purchased from IEDS.

On the other hand, renewable energy should be consumed as much as possible for low carbon and environmental protection.

Renewable energy output constraints:

$$0\le p_{i^{wind},n,t}^{wind}\le {\overline{p}}_{i^{wind},n,t}^{wind} \forall i^{wind},\forall n,\forall t$$

(31)

$$p_{i^{wind},n,t}^{wind,cut}={\overline{p}}_{i^{wind},n,t}^{wind}-p_{i^{wind},n,t}^{wind} \forall i^{wind},\forall n,\forall t$$

(32)

$$0\le p_{i^{solar},n,t}^{solar}\le {\overline{p}}_{i^{solar},n,t}^{solar} \forall i^{solar},\forall n,\forall t$$

(33)

$$p_{i^{solar},n,t}^{solar,cut}={\overline{p}}_{i^{solar},n,t}^{solar}-p_{i^{solar},n,t}^{solar} \forall i^{solar},\forall n,\forall t$$

(34)

Constraints (31) and (32) describe the power outputs of wind generators. Constraints (33) and (34) describe the power outputs of photovoltaic plants.

Finally, DMESs ensure that the load shedding should be less than the total load.

Load shedding constraints:

$$0\le p_{n,t}^{shed}\le p_{n,t}^{load} \forall n,\forall t$$

(35)

$$0\le f_{n,t}^{shed}\le f_{n,t}^{load} \forall n,\forall t$$

(36)

$$0\le q_{n,t}^{shed}\le q_{n,t}^{load} \forall n,\forall t$$

(37)

Constraints (35)–(37) limit the electricity, gas and heat load shedding amount, respectively.

The proposed scheduling model of DMES is a mixed-integer linear programming (MILP), which can be solved with the solvers, such as cplex, gurobi.

The proposed scheduling model is solved in each DMES simultaneously. Next, each DMES calculates the excess and deficit electrical power and natural gas. The excess electrical power includes unused CHP unit and battery capacities. Besides, the deficit electrical power refers to unsupplied electrical demands. The calculation formulation is presented as follows:

$$\begin{array}{ll}{C}_{n,t}^{max}=& {\displaystyle \sum _{i^{CHP}=1}^{N_n^{CHP}}\min \{[{\overline{p}}_{i^{CHP},n}^{CHP}-{\widehat{p}}_{i^{CHP},n,t}^{CHP}], [{\widehat{p}}_{i^{CHP},n,t-1}^{CHP}+\Delta {\overline{p}}_{i^{CHP},n}^{CHP}-{\widehat{p}}_{i^{CHP},n,t}^{CHP}]},\\ & [{\tau }_{i^{CHP},n}{\eta }_{i^{CHP},n}^{CHP,ele}({\overline{f}}_{i^{gas},n}^{dch}-{\widehat{f}}_{i^{gas},n,t}^{dch})], [{\tau }_{i^{CHP},n}{\eta }_{i^{CHP},n}^{CHP,ele}({\widehat{E}}_{i^{gas},n,t}^{gas}{\eta }_{i^{gas},n}^{gas,dch}-{\widehat{f}}_{i^{gas},n,t}^{dch})]\}\end{array}$$

(38)

$$B_{n,t}^{max}={\displaystyle \sum _{i^{bat}=1}^{N_n^{bat}}\min \{[{\overline{p}}_{i^{bat},n}^{dch}-{\widehat{p}}_{i^{bat},n,t}^{dch}],[{\widehat{E}}_{i^{bat},n,t}^{bat}{\eta }_{i^{bat},n}^{bat,dch}-{\widehat{p}}_{i^{bat},n,t}^{dch}]\}}$$

(39)

$$L_{n,t}^e={\widehat{p}}_{n,t}^{shed}$$

(40)

where symbol ^ denotes that the corresponding variables are optimal values from the proposed scheduling model. The further generating capacity $C_{n,t}^{max}$ is the sum of the unused capacities of each CHP unit in DMES n, which is constrained by capacity, the ramp-up of CHP units and discharging gas, and gas stored in gasholders. The further discharging capacity $B_{n,t}^{max}$ is the sum of the unused capacities of each battery in DMES n, which is constrained by discharging power, and electricity stored in batteries. $L_{n,t}^e$ is the total electrical load shedding in DMES n. The excess natural gas includes unused gasholder capacities, and the deficit natural gas refers to unsupplied gas demands. They are calculated as follows:

$$G_{n,t}^{max}={\displaystyle \sum _{i^{gas}=1}^{N_n^{gas}}\min \{[{\overline{f}}_{i^{gas},n}^{dch}-{\widehat{f}}_{i^{gas},n,t}^{dch}],[{\widehat{E}}_{i^{gas},n,t}^{gas}{\eta }_{i^{gas},n}^{gas,dch}-{\widehat{f}}_{i^{gas},n,t}^{dch}]\}}$$

(41)

$$L_{n,t}^g={\widehat{f}}_{n,t}^{shed}$$

(42)

The further discharging capacity $G_{n,t}^{max}$ is the sum of the unused capacities of each gasholder in DMES n, which is constrained by discharging gas, and gas stored in gasholders. $L_{n,t}^g$ is the total natural gas load shedding in DMES n.

After the excess and deficit electrical power and natural gas is calculated, each DMES broadcasts this information to neighboring DMESs at the end of stage I.

3.2. The Consensus Algorithm in Networked DMESs

In stage II, the amount of energy exchange supported by each other will be determined between the DMESs. This paper defines an excess DMES as a DMES that has excess energy after supplying its own load, and a deficit DMES as a DMES that has load shedding and needs energy support from other DMESs. If a centralized approach is used for calculating energy exchange, it will suffer high communication infrastructure costs and low reliability. Alternatively, a consensus algorithm can be applied in which the allocated proportions of energy exchange are computed locally at each DMES by information exchange only with neighboring DMESs. The essence of the consensus algorithm is to update the variable of local DMESs according to the variable of adjacent DMESs, iteratively. The iteration variable of DMESs can be linearly updated as follows:

$${\Phi }_{n,t}(k+1)={\Phi }_{n,t}(k)+\xi \cdot {\displaystyle \sum _{m\in {\Omega }_n}({\Phi }_{m,t}(k)-{\Phi }_{n,t}(k))} \forall n,\forall t$$

(43)

where k is the iteration times. ξ is the step size, which satisfies $\xi \in (0,1/\Theta )$. $\Theta$ is the maximum node out-degree of the communication network [26,27]. ${\Omega }_n$ is set of neighbors of DMES n. ${\Phi }_{n,t}(k)$ is the kth iteration variable, which represents $C_{n,t}^{max}$, $B_{n,t}^{max}$, $G_{n,t}^{max}$, $L_{n,t}^e$ and $L_{n,t}^g$. Based on research in [27], the convergence of the consensus algorithm is related to the eigenvalues of the Laplacian matrix representing the communication network. After enough runs of iteration, variables ${\Phi }_{n,t}(k)$ will converge to the average value ${\tilde{\Phi }}_t$:

$${\tilde{\Phi }}_t=\frac{1}{N^{DMES}}{\displaystyle \sum _{n=1}^{N^{DMES}}{\Phi }_{n,t}(0)} \forall n,\forall t$$

(44)

where ${\Phi }_{n,t}(0)$ denotes the initial value of variable. Actually, the consensus algorithm is to calculate the average value under the condition that only the local and adjacent variables are known. It protects the privacy of DMESs, since only excess and deficit energy are visible. Furthermore, a distributed calculation method is applied for developing scheduling plan of each DMESs in parallel and only local communication between neighboring DMESs is involved, so it has low communication complexity and cost.

Next, the allocated proportions can be calculated as follows:

$${\mu }_t^e=\{\begin{array}{ll}\frac{{\tilde{L}}_t^e}{{\tilde{C}}_t+{\tilde{B}}_t}& {\tilde{L}}_t^e<{\tilde{C}}_t+{\tilde{B}}_t\\ \frac{{\tilde{C}}_t+{\tilde{B}}_t}{{\tilde{L}}_t^e}& {\tilde{L}}_t^e\ge {\tilde{C}}_t+{\tilde{B}}_t\end{array}$$

(45)

$${\mu }_t^g=\{\begin{array}{ll}\frac{{\tilde{L}}_t^g}{{\tilde{G}}_t}& {\tilde{L}}_t^g<{\tilde{G}}_t\\ \frac{{\tilde{G}}_t}{{\tilde{L}}_t^g}& {\tilde{L}}_t^g\ge {\tilde{G}}_t\end{array}$$

(46)

where ${\mu }_t^e$ and ${\mu }_t^g$ are allocated proportions of electricity and natural gas, respectively. ${\tilde{C}}_t$, ${\tilde{B}}_t$, ${\tilde{L}}_t^e$, ${\tilde{G}}_t$ and ${\tilde{L}}_t^g$ are convergent values calculated by (43). According to the allocated proportions, the exported or imported electricity and natural gas for each DMES can be calculated.

If excess energy is more than deficit energy in networked DMESs at hour t (${\tilde{L}}_t^e<{\tilde{C}}_t+{\tilde{B}}_t$, ${\tilde{L}}_t^g<{\tilde{G}}_t$), allocated proportions are applied to calculated the exported energy of excess DMESs as follows:

$$p_{n,t}^{out}={\mu }_t^e(C_{n,t}^{\max }+B_{n,t}^{max}) \forall n,\forall t$$

(47)

$$f_{n,t}^{out}={\mu }_t^g{G}_{n,t}^{max} \forall n,\forall t$$

(48)

where $p_{n,t}^{out}$ and $f_{n,t}^{out}$ are the exported electricity and natural gas of excess DMESs, respectively. The electricity and natural gas are enough to alleviate load shedding in networked DMESs after disaster and improve the resilience of IEDS.

On the other hand, if the excess energy is less than the deficit energy in networked DMESs at hour t (${\tilde{L}}_t^e\ge {\tilde{C}}_t+{\tilde{B}}_t$, ${\tilde{L}}_t^g\ge {\tilde{G}}_t$), allocated proportions are applied to the calculated imported energy of deficit DMESs as follows:

$$p_{n,t}^{in}={\mu }_t^e{L}_{n,t}^e \forall n,\forall t$$

(49)

$$f_{n,t}^{in}={\mu }_t^g{L}_{n,t}^g \forall n,\forall t$$

(50)

where $p_{n,t}^{in}$ and $f_{n,t}^{in}$ are the imported electricity and natural gas of deficit DMESs, respectively. When excess energy is not enough, the energy support for deficit DMESs is based on their amount of load shedding.

In general, it is fair to all DMESs, since the allocated exported or imported energy for each DMES is set according to its excess capacity or amount of load shedding. After DMESs determine the amount of energy exchange by the consensus algorithm, the scheduling plan developed in stage I will be updated in two steps as shown in Figure 4. And then, the updated scheduling plan will be implemented at the end of stage II.

4. Numerical Results

4.1. Test System

The proposed hierarchical management strategy is examined on a test system with five DMESs as shown in Figure 2. The information on CHP units, EHPs, boilers, batteries, gasholders, heat storages and price is summarized in

Table 1. Information on combined heat and power (CHP) units.

n	${\mathit{\tau }}_{{\mathit{i}}^{\mathit{C}\mathit{H}\mathit{P}},\mathit{n}}$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{C}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{C}\mathit{H}\mathit{P},\mathit{e}\mathit{l}\mathit{e}}$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{C}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{C}\mathit{H}\mathit{P},\mathit{h}\mathit{e}\mathit{a}\mathit{t}}$	${\overline{\mathit{p}}}_{{\mathit{i}}^{\mathit{C}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{C}\mathit{H}\mathit{P}} (\mathbf{MW})$	$\mathbf{\Delta }{\overline{\mathit{p}}}_{{\mathit{i}}^{\mathit{C}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{C}\mathit{H}\mathit{P}} (\mathbf{MW})$	$\mathbf{\Delta }{\underset{¯}{\mathit{p}}}_{{\mathit{i}}^{\mathit{C}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{C}\mathit{H}\mathit{P}} (\mathbf{MW})$
1, 2, 3, 5	0.8	0.75	0.75	4	2	−2
4	0.8	0.75	0.75	5	2	−2

,

Table 2. Information on electrical heat pumps (EHPs).

n	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{E}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{E}\mathit{H}\mathit{P}}$	${\overline{\mathit{q}}}_{{\mathit{i}}^{\mathit{E}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{E}\mathit{H}\mathit{P}} (\mathbf{MBtu}/\mathbf{h})$	$\mathbf{\Delta }{\overline{\mathit{q}}}_{{\mathit{i}}^{\mathit{E}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{E}\mathit{H}\mathit{P}} (\mathbf{MBtu}/\mathbf{h})$	$\mathbf{\Delta }{\underset{¯}{\mathit{q}}}_{{\mathit{i}}^{\mathit{E}\mathit{H}\mathit{P}},\mathit{n}}^{\mathit{E}\mathit{H}\mathit{P}} (\mathbf{MBtu}/\mathbf{h})$
1, 2, 3, 4, 5	0.75	2	0.8	−0.8

,

Table 3. Information on boilers.

n	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{b}\mathit{o}\mathit{i}},\mathit{n}}^{\mathit{b}\mathit{o}\mathit{i}}$	${\overline{\mathit{q}}}_{{\mathit{i}}^{\mathit{b}\mathit{o}\mathit{i}},\mathit{n}}^{\mathit{b}\mathit{o}\mathit{i}} (\mathbf{MBtu}/\mathbf{h})$	$\mathbf{\Delta }{\overline{\mathit{q}}}_{{\mathit{i}}^{\mathit{b}\mathit{o}\mathit{i}},\mathit{n}}^{\mathit{b}\mathit{o}\mathit{i}} (\mathbf{MBtu}/\mathbf{h})$	$\mathbf{\Delta }{\underset{¯}{\mathit{q}}}_{{\mathit{i}}^{\mathit{b}\mathit{o}\mathit{i}},\mathit{n}}^{\mathit{b}\mathit{o}\mathit{i}} (\mathbf{MBtu}/\mathbf{h})$
1, 2, 3, 4, 5	0.8	2	1	−1

,

Table 4. Information on batteries.

n	${\overline{\mathit{E}}}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{b}\mathit{a}\mathit{t}} (\mathbf{MWh})$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{b}\mathit{a}\mathit{t},\mathit{c}\mathit{h}}$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{b}\mathit{a}\mathit{t},\mathit{d}\mathit{c}\mathit{h}}$	${\overline{\mathit{p}}}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{c}\mathit{h}} (\mathbf{MW})$	${\overline{\mathit{p}}}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{d}\mathit{c}\mathit{h}} (\mathbf{MW})$	${\mathit{\theta }}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{b}\mathit{a}\mathit{t},\mathit{c}\mathit{h}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\theta }}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{b}\mathit{a}\mathit{t},\mathit{d}\mathit{c}\mathit{h}} (\mathbf{€}/\mathbf{MWh})$
1	10	0.9	0.9	1.5	1.5	1	1
2	6	0.9	0.9	1	1	1	1
3	4	0.9	0.9	0.8	1	1	1
4	12	0.9	0.9	1.5	2	1	1
5	8	0.9	0.9	1.5	1.5	1	1

,

Table 5. Information on gasholders.

n	${\overline{\mathit{E}}}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}},\mathit{n}}^{\mathit{g}\mathit{a}\mathit{s}} (\mathbf{kcf})$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}},\mathit{n}}^{\mathit{g}\mathit{a}\mathit{s},\mathit{c}\mathit{h}}$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}},\mathit{n}}^{\mathit{g}\mathit{a}\mathit{s},\mathit{d}\mathit{c}\mathit{h}}$	${\overline{\mathit{f}}}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}},\mathit{n}}^{\mathit{c}\mathit{h}} (\mathbf{kcf}/\mathbf{h})$	${\overline{\mathit{f}}}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}},\mathit{n}}^{\mathit{d}\mathit{c}\mathit{h}} (\mathbf{kcf}/\mathbf{h})$	${\mathit{\theta }}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}},\mathit{n}}^{\mathit{g}\mathit{a}\mathit{s},\mathit{c}\mathit{h}} (\mathbf{€}/\mathbf{kcf})$	${\mathit{\theta }}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}},\mathit{n}}^{\mathit{g}\mathit{a}\mathit{s},\mathit{d}\mathit{c}\mathit{h}} (\mathbf{€}/\mathbf{kcf})$
1	120	0.95	0.95	15	25	1	1
2	60	0.95	0.95	10	20	1	1
3	40	0.95	0.95	10	20	1	1
4	160	0.95	0.95	30	30	1	1
5	100	0.95	0.95	15	20	1	1

,

Table 6. Information on heat storages.

n	${\overline{\mathit{E}}}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}} (\mathbf{MBtu})$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t},\mathit{c}\mathit{h}}$	${\mathit{\eta }}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t},\mathit{d}\mathit{c}\mathit{h}}$	${\overline{\mathit{q}}}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{c}\mathit{h}} (\mathbf{MBtu}/\mathbf{h})$	${\overline{\mathit{q}}}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{d}\mathit{c}\mathit{h}} (\mathbf{MBtu}/\mathbf{h})$	${\mathit{\theta }}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t},\mathit{c}\mathit{h}} (\mathbf{€}/\mathbf{MBtu})$	${\mathit{\theta }}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}},\mathit{n}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t},\mathit{d}\mathit{c}\mathit{h}} (\mathbf{€}/\mathbf{MBtu})$
1	120	0.95	0.95	15	25	1	1
2	60	0.95	0.95	10	20	1	1
3	40	0.95	0.95	10	20	1	1
4	160	0.95	0.95	30	30	1	1
5	100	0.95	0.95	15	20	1	1

and

Table 7. Information on price.

t	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{e}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{g}} (\mathbf{€}/\mathbf{kcf})$	t	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{e}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{g}} (\mathbf{€}/\mathbf{kcf})$	t	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{e}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{g}} (\mathbf{€}/\mathbf{kcf})$	t	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{e}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\lambda }}_{\mathit{t}}^{\mathit{g}} (\mathbf{€}/\mathbf{kcf})$
1	67.09	16.93	7	71.34	22.15	13	70.88	20.25	19	72.56	24.60
2	66.89	16.60	8	71.04	22.60	14	71.10	20.25	20	71.19	24.43
3	67.64	15.75	9	71.41	23.77	15	71.52	20.48	21	70.42	20.53
4	68.64	16.85	10	71.61	23.73	16	72.27	24.50	22	70.16	17.98
5	69.54	20.54	11	71.23	23.93	17	72.82	24.65	23	69.36	17.88
6	70.16	21.30	12	71.17	23.75	18	73.15	25.00	24	66.92	17.18

.

The forecasted output of DGs in each DMES, which is the sum of wind output and photovoltaic output, is shown in Figure 5.

And then, the electrical load, gas load and heat load of DMESs are shown in Figure 6a–c, respectively.

It is assumed that DMESs receive a disaster alert at the 7th hour, and disaster attacks the IEDS at the 17th hour. This means that the normal operation lasts 6 h, the preventive operation lasts 10 h and the resilient operation lasts 8 h. In other words, T₁, T₂ and T₃ are set to 6, 10 and 8 h. The remaining penalty cost parameters are summarized in

Table 8. Information on penalty cost parameters.

${\mathit{\pi }}_{\mathit{n}}^{\mathit{e}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\pi }}_{\mathit{n}}^{\mathit{g}} (\mathbf{€}/\mathbf{kcf})$	${\mathit{\pi }}_{\mathit{n}}^{\mathit{h}} (\mathbf{€}/\mathbf{MBtu})$	${\mathit{\kappa }}_{{\mathit{i}}^{\mathit{w}\mathit{i}\mathit{n}\mathit{d}}}^{\mathit{w}\mathit{i}\mathit{n}\mathit{d}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\kappa }}_{{\mathit{i}}^{\mathit{s}\mathit{o}\mathit{l}\mathit{a}\mathit{r}}}^{\mathit{s}\mathit{o}\mathit{l}\mathit{a}\mathit{r}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\omega }}_{{\mathit{i}}^{\mathit{b}\mathit{a}\mathit{t}}}^{\mathit{b}\mathit{a}\mathit{t}} (\mathbf{€}/\mathbf{MWh})$	${\mathit{\omega }}_{{\mathit{i}}^{\mathit{g}\mathit{a}\mathit{s}}}^{\mathit{g}\mathit{a}\mathit{s}} (\mathbf{€}/\mathbf{kcf})$	${\mathit{\omega }}_{{\mathit{i}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}}}^{\mathit{h}\mathit{e}\mathit{a}\mathit{t}} (\mathbf{€}/\mathbf{MBtu})$
5000	2000	1000	80	80	100	50	50

.

4.2. Results and Analysis

First, all DMESs have sufficient energy purchased from IEDSs in the normal and preventive operation mode, the electrical power and natural gas purchased curves are shown in Figure 7.

As shown in Figure 7, the energy purchased increases and maintains a high level in the preventive operation mode, because it not only supplies the current load, but also charges the energy storages. In addition, in order to discuss the proposed hierarchical management strategy more clearly, the energy storages should be filled during the preventive operation mode, so that there is sufficient energy stored to schedule between networked DMESs. The decline of energy purchased curves at the end of the preventive operation mode indicates that energy storages were full.

The resilient operation mode covers 8 h, from the 17th hour to the 24th hour. The excess electricity and natural gas of DMESs are summarized in

Table 9. The excess electricity and natural gas.

Time Slot t (h)		1	2	3	4	5	6	7	8
$C_{n,t}^{max}$ + $B_{n,t}^{max}$ (MW)	DMES 1	1.34	1.72	2.01	2.51	2.56	2.76	0	0
	DMES 2	0.97	0.68	1.43	0	0	0	0	0
	DMES 3	0.48	0.84	0	0	0	0	0	0
	DMES 4	1.16	2.05	1.80	1.91	3.14	1.98	0.77	0
	DMES 5	1.17	0.97	1.52	1.64	1.16	0.08	0	0
$G_{n,t}^{max}$ (kcf/h)	DMES 1	13.52	11.86	8.89	8.48	6.82	2.23	0	0
	DMES 2	6.84	3.67	5.68	0	0	0	0	0
	DMES 3	3.89	4.52	0	0	0	0	0	0
	DMES 4	16.61	16.61	13.16	9.14	8.48	1.60	0	0
	DMES 5	7.40	5.23	6.06	5.53	3.15	0.06	0	0

. The deficit electricity and natural gas (load shedding) of DMESs are summarized in

Table 10. The deficit electricity and natural gas.

Time Slot t (h)		1	2	3	4	5	6	7	8
$L_{n,t}^e$ (MW)	DMES 1	0	0	0	0	0	0	1.70	1.18
	DMES 2	0	0	0	0.20	1.18	1.82	2.00	1.40
	DMES 3	0	0	1.75	2.05	2.23	2.27	2.05	1.09
	DMES 4	0	0	0	0	0	0	0	2.18
	DMES 5	0	0	0	0	0	0	2.26	1.94
$L_{n,t}^g$ (kcf/h)	DMES 1	0	0	0	0	0	0	2.15	5.84
	DMES 2	0	0	0	0	7.84	7.39	7.80	6.90
	DMES 3	0	0	2.13	8.01	7.86	7.46	7.81	7.32
	DMES 4	0	0	0	0	0	0	0	8.04
	DMES 5	0	0	0	0	0	0	7.16	7.53

.

Comparing Table 9 and Table 10, we can see that excess energy and deficit energy are almost complementary. Approximately, each DMES can supply its own loads at time slot 1~2, and DMES 1, DMES 4 and DMES 5 have more excess energy to support DMES 2 and DMES 3 at time slot 3~6. However, with the consumption of energy, the DMESs have to interrupt their loads.

As shown in Figure 2, the DMESs can exchange information through the cyber communication network according to the route 1 ↔ 2 ↔ 3 ↔ 4 ↔ 5 ↔ 1. Since there is no energy exchange at time slot 1~2 and 8, it mainly shows the iterative process of the consensus algorithm at time slot 3~7 in Figure 8, especially the excess energy at time slot 3~6 and deficit energy at time slot 7, which are representative.

As shown in Figure 8, after five runs of iteration, the variables nearly converge to the convergent value, which demonstrates that the consensus algorithm has efficient convergence. In summary, the convergent values and allocated proportions calculated by (45) and (46) are summarized in

Table 11. The convergent values and allocated proportions.

Time Slot t (h)	1	2	3	4	5	6	7	8
${\tilde{C}}_t$ + ${\tilde{B}}_t$ (MW)	1.02	1.25	1.35	1.21	1.37	0.96	0.15	0
${\tilde{L}}_t^e$ (MW)	0	0	0.35	0.45	0.68	0.82	1.50	1.56
${\mu }_t^e$	0	0	0.2588	0.3715	0.4968	0.8496	0.0963	0
${\tilde{G}}_t$ (kcf/h)	9.65	8.38	6.76	4.63	3.69	0.78	0	0
${\tilde{L}}_t^g$ (kcf/h)	0	0	0.43	1.60	3.14	2.97	4.98	7.13
${\mu }_t^g$	0	0	0.063	0.346	0.8461	0.2621	0	0

.

According to the allocated proportions, exported or imported electricity and natural gas for each DMES at time slot 3~6 are shown in Figure 9.

Figure 9 shows that the amount of exchange energy between DMESs is calculated for each DMES by the consensus algorithm. And then, the scheduling plan is updated. As noted earlier, the updated process is divided into two steps: the electrical load curtailment is mitigated in the first step and the gas load curtailment is mitigated in the second step. To clarify, the first step is called electricity update and the second step is called gas update in this paper. The load shedding of networked DMESs in the initial scheduling plan and the two updated scheduling plans are shown in Figure 10.

As shown in Figure 10, compared with the initial scheduling plan, the updated scheduling plan has less load shedding. In other words, after implementing the hierarchical management strategy, the networked DMESs can serve more consumers, avoiding severe loading shedding. If the IEDS has been repaired in time slot 1~5, the networked DMESs will not interrupt the load supply at all, which shows that IEDSs have the ability to absorb and withstand extreme events, so the hierarchical management strategy can enhance the resilience of IEDS.

5. Conclusions

IEDSs are efficient, clean and sustainable energy systems. However, increasingly frequent natural disasters pose a severe threat to the security of IEDSs. Hence, this paper studied the resilience of IEDSs against extreme events. A hierarchical management strategy was proposed, which includes the normal, preventive and resilient operation mode of DMESs. In the normal operation mode, the objective is to minimize costs. In the preventive operation mode, the objective is to maximize the stored energy. It is noted that the resilient operation mode is divided into two stages. In the first stage, each DMES develops its own scheduling plan with the objective of minimizing load shedding. And then, the excess and deficit energy is broadcasted to neighboring DMESs by the LLCs. In the second stage, a consensus algorithm is applied to calculate energy exchange between the networked DMESs. Next, the scheduling plan is updated in two steps and implemented. Finally, the numerical results demonstrate the effectiveness of the proposed hierarchical management strategy and algorithm. The main conclusions are summarized as follows:

(1)

This paper co-optimizes the electricity, natural gas and thermal energy for three operation modes at the distribution and utilization level, which promotes efficient and sustainable energy consumption.

(2)

The hierarchical management strategy can dispatch energy exchange fairly and reduce load shedding efficiently after extreme events, which demonstrates its enhancement for IEDS resilience.

(3)

The privacy of each DMES is respected since only information on excess and deficit energy is visible to neighboring DMESs through the consensus algorithm, which has efficient convergence.

Future research directions include the topology formation of IEDS, flexible resources for sectionalized restoration and the self-healing mode of IEDS.