*Article* **Distributionally Robust Joint Chance-Constrained Dispatch for Electricity–Gas–Heat Integrated Energy System Considering Wind Uncertainty**

**Hui Liu 1,2, Zhenggang Fan 1,2, Haimin Xie 1,2,\* and Ni Wang 1,2**


**Abstract:** With the increasing penetration of wind power, the uncertainty associated with it brings more challenges to the operation of the integrated energy system (IES), especially the power subsystem. However, the typical strategies to deal with wind power uncertainty have poor performance in balancing economy and robustness. Therefore, this paper proposes a distributionally robust joint chance-constrained dispatch (DR-JCCD) model to coordinate the economy and robustness of the IES with uncertain wind power. The optimization dispatch model is formulated as a two-stage problem to minimize both the day-ahead and the real-time operation costs. Moreover, the ambiguity set is generated using Wasserstein distance, and the joint chance constraints are used to ensure that the safety constraints (e.g., ramping limit and transmission limit) can be satisfied jointly under the worst-case probability distribution of wind power. The model is remodeled as a mixed-integer tractable programming issue, which can be solved efficiently by ready-made solvers using linear decision rules and linearization methods. Case studies on an electricity–gas–heat regional integrated system, which includes a modified IEEE 24-bus system, 20 natural gas-nodes, and 6 heat-node system, are investigated for verification. Numerical simulation results demonstrate that the proposed DR-JCCD approach effectively coordinates the economy and robustness of IES and can offer operators a reasonable energy management scheme with an acceptable risk level.

**Keywords:** distributionally robust optimization (DRO); integrated energy system (IES); joint chance constraints; linear decision rules (LDRs); Wasserstein distance

**1. Introduction**

In order to achieve the 1.5 ◦C temperature control target set by the Paris Climate Agreement [1], the proportion of global power generation via renewable sources will continue to rise. By 2050, renewable energy is expected to account for 86% of the power generation source. Wind power, in particular, will meet more than 35% of power demand and become the main source of power generation at that time [2]. However, as the penetration of renewable energy sources (RESs) increases, the power network will be exposed to greater risks due to the uncertainty of RESs. Therefore, there is an urgent need to improve the flexibility of power systems or mitigate the variability. Constructing regional integrated energy systems (IESs) has been proved as an effective way to provide more flexibility to accommodate renewable sources and reduce the impact of uncertainty on the power system [3].

Many researches have concentrated on the optimal dispatch of IESs to cope with the uncertainties associated with renewable energy. Two common strategies include stochastic programming (SP) [4–9] and robust optimization (RO) [10,11]. To handle uncertainties of load demand and renewable energy, Yong et al. [4] propose a low-carbon optimal

**Citation:** Liu, H.; Fan, Z.; Xie, H.; Wang, N. Distributionally Robust Joint Chance-Constrained Dispatch for Electricity–Gas–Heat Integrated Energy System Considering Wind Uncertainty. *Energies* **2022**, *15*, 1796. https://doi.org/10.3390/en15051796

Academic Editor: Luigi Fortuna

Received: 23 January 2022 Accepted: 10 February 2022 Published: 28 February 2022

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

**Copyright:** © 2022 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/).

stochastic operation model using power-to-gas technology. In a building energy system, a multistage-based scenario-driven approach is proposed to deal with solar power uncertainty and nonschedulable load uncertainty [5]. However, stochastic programming either relies on scene samples to approximate deterministic distributions [6] or assumes a predefined probability distribution that random variables follow [7]. As a result, it imposes a substantial computational burden on optimization [8] and adds difficulty to scenario selection [9]. Compared to stochastic programming, the robust optimization approach does not require any assumptions about wind power probability distribution, because it can ensure system-operational robustness by using uncertainty sets to make optimal decisions under the worst renewable fluctuation cases [10]. Robust optimization can provide a more reliable scheme while considering wind power uncertainty [11], but it compromises system cost-effectiveness and may result in over-conservative solutions.

Distributionally robust optimization (DRO), an effective strategy to overcome weaknesses of stochastic programming and robust optimization, loads all possible wind power probability distributions information into an ambiguous set to incorporate uncertain wind power distributions. In addition, DRO can ensure that all the possible wind power probability distributions in the ambiguous set are met by making the best decision under the worst-case probability distribution [12]. There are several studies on distributionally robust energy models based on moment-based ambiguity sets such as mean vector, covariance matrix, and higher moment information [13–15]. Reference [13] develops a two-stage voltage and natural gas pipe pressure management model for photovoltaic power in IES, where the photovoltaic power uncertainty is modeled by an ambiguity set containing the first-order moment and second-order moment information. Using the same moment information, a distributionally robust optimal power flow problem is formulated to solve renewable energy and load uncertainties [14]. Due to the distributions of renewable forecast errors practically containing higher moment information, the first two moment-based ambiguous sets may cause unnecessary conservatism. Therefore, Reference [15] proposes a DRO model for an energy hub system with an energy storage function, where the ambiguity set contains the first two moments and multimodal information of photovoltaic power forecasting errors. Despite all this, there may be the same moment information among different distributions, which makes it difficult to determine the worst-case probability distribution.

The other approach to characterizing ambiguity sets is on account of the statistical distance between the true probability distribution and possible probability distributions. One type of discrepancy-based ambiguity set is established by Kullback–Leibler divergence. To guarantee the safe operation of the natural gas system under hydrogen injection when utilizing power-to-gas technology, Ref. [16] develops a natural gas security DRO programming for IES using a Kullback–Leibler divergence-based ambiguity set to capture wind uncertainty. However, only in the circumstances that potential distributions are supported on a set of limited values, the Kullback–Leibler divergence-based ambiguity set can be observed through historical data [17]. In contrast, the ambiguity sets based on Wasserstein distance, which include all possible probability distributions that have a narrow gap with the discrete empirical distribution, are introduced and have been increasingly used. With a Wasserstein-distance-based ambiguity set to deal with renewable energy uncertainties, a power-flow DRO problem with multi-stage feedback policies is formulated in [18]. In [19], considering dynamic line rating and operational risk, a power flow DRO approach is established, which constructs the ambiguity set via combining the moment information and Wasserstein distance. To avoid the calculation issue arising from a large number of historical data sets, Ref. [20] proposes a distance-based aggregation method, and the Wassersteindistance ambiguity set is introduced to a distributionally robust unit commitment problem. Regarding wind power uncertainty, the Wasserstein-distance-based ambiguity set has shown good performance in both finite-sample guarantees and confidence sets.

Although the DRO method with chance-constrained problems has been extensively addressed in optimal power flow [21–23], its research on energy optimization and management is still in the early stage. A distributionally robust individual chance-constrained

energy dispatch model is put forward for an islanded heat and electricity system in [24] while considering the uncertain renewable generation. However, due to the confidence levels considered separately, the individual chance constraints will result in high-risk costs and may even result in confidence levels as low as 0% for any individual constraint [25]. Therefore, a joint chance-constrained DRO model is proposed for the combined electricity and natural gas system to address renewable energy uncertainty while using the ambiguity set with the confidence bands of the true density function [26]. Joint chance constraints can improve the simultaneous satisfaction of multiple safety conditions with a high probability, but the ambiguity set only includes marginal distribution information, which will result in a conservative solution.

Therefore, a two-stage distributionally robust joint chance-constrained dispatch (DR-JCCD) model is proposed for the electricity–gas–heat IES with the Wasserstein distancebased ambiguity set, considering the wind power uncertainty. The main contributions of this paper are as follows:


The remainder of this paper is organized as follows. Section 2 presents the detailed mathematical formulation for both the day-ahead and the real-time operation, where the nomenclature could be found in Appendix A. Section 3 develops effective approximate and re-modeled schemes for the distributionally robust joint chance-constrained model under the Wasserstein ambiguity set. Numerical results are shown in Section 4. Finally, conclusions and future work are given in Section 5.

#### **2. DR-JCCD Modeling of IES**

## *2.1. Framework for Electricity–Gas–Heat IES*

As illustrated in Figure 1, a typical framework of the electricity, gas, and heat IES uses a constrained transmission infrastructure to coordinate power generation and natural gas resources. The configuration of multiple energy carriers employed in the energy hub includes Combined Heat and Power (CHP), gas-fired generations, and an electric boiler. CHP transforms natural gas into electricity and heat concurrently. Natural gas-fired units transform natural gas into electricity, allowing them to respond swiftly to power fluctuations. An electric boiler is introduced to supply enough heat power flexibly, and three different energy sources are used to meet local electricity, gas, and heat demands.

**Figure 1.** Multi‐energy system framework for the IES dispatch. **Figure 1.** Multi-energy system framework for the IES dispatch.

#### *2.2. Formulation of Day-Ahead Operation*

*2.2. Formulation of Day‐Ahead Operation* The day‐ahead operation schedules the output and reserve capacity of thermal gen‐ erations, natural gas‐fired generations, and CHP. In addition, it determines the natural gas output of gas sources. The objective of the day‐ahead operation is to minimize the expected operation cost, including generation cost, reserve capacity costs for traditional units, natural gas‐fired units, and CHP, as well as the cost of consuming natural gas from The day-ahead operation schedules the output and reserve capacity of thermal generations, natural gas-fired generations, and CHP. In addition, it determines the natural gas output of gas sources. The objective of the day-ahead operation is to minimize the expected operation cost, including generation cost, reserve capacity costs for traditional units, natural gas-fired units, and CHP, as well as the cost of consuming natural gas from the source, as shown in Equation (1)

$$\min\_{\begin{array}{c} P^{\text{DA}}, P^{\text{RT}}(\boldsymbol{\zeta}),\\ r^{+}, r^{-} \end{array}} \sum\_{\begin{subarray}{c} \mathcal{I}\_{\mathcal{I}} \in I\_{\mathcal{I}}, i\_{\mathcal{S}} \in I\_{\mathcal{S}} \\ t \in T \end{subarray}} \lambda\_{\{\cdot\}} P^{DA}\_{\{\cdot\}, t} + \lambda\_{i\_{\mathcal{S}}} \mathsf{G}\_{\text{gas}, t} + \lambda\_{\{\cdot\}}^{+} r^{+}\_{\{\cdot\}, t} + \lambda\_{\{\cdot\}}^{-} r^{-}\_{\{\cdot\}, t} + \mathsf{E}^{\text{P}}[\gamma\_{\mathcal{E}} P^{\text{P}}\_{\{\cdot\}, t}(\boldsymbol{\zeta})] \tag{1}$$

The detail constraints for day‐ahead operation include constraints of the power sys‐ tem, natural‐gas system, heating system, and multiple energy converters, as listed in Equations (2)–(29). The detail constraints for day-ahead operation include constraints of the power system, natural-gas system, heating system, and multiple energy converters, as listed in Equations (2)–(29).

#### 2.2.1. Constraints for the Power System

2.2.1. Constraints for the Power System The reserve capacity from thermal generations, natural gas‐fired generations, and CHP is shown in Equations (2) and (3), followed by their power output limits in Equations (4) and (5). The adjustments of multiple energy devices for the uncertain wind power fore‐ casting errors are limited by Equation (6), which must be within the reserve capacity range. The ramping up/down constraints are restricted by Equations (7) and (8) [17]. En‐ ergy balance is presented in Equation (9). Constrain Equation (10) is used to ensure that the power flows are well within the capacity limits of the transmission lines. Note that The reserve capacity from thermal generations, natural gas-fired generations, and CHP is shown in Equations (2) and (3), followed by their power output limits in Equations (4) and (5). The adjustments of multiple energy devices for the uncertain wind power forecasting errors are limited by Equation (6), which must be within the reserve capacity range. The ramping up/down constraints are restricted by Equations (7) and (8) [17]. Energy balance is presented in Equation (9). Constrain Equation (10) is used to ensure that the power flows are well within the capacity limits of the transmission lines. Note that {·} is the index and set of thermal units, gas-fired generations, and CHP.

$$0 \le r^+\_{\{\cdot\},t} \le \mathcal{R}^+\_{\{\cdot\}'} \qquad \{\cdot\} = i\_{\text{e}}, \text{gg}, \text{ch}p \tag{2}$$

(1)

{ }

$$0 \le r\_{\{\cdot\},t}^{-} \le R\_{\{\cdot\}'}^{-} \qquad \{\cdot\} = i\_{\text{e}} \text{ g\,g\text{,}} \\ \text{d}p \tag{3}$$

$$P\_{\{\cdot\},t}^{\text{DA}} + r\_{\{\cdot\},t}^{+} \le P\_{\{\cdot\},t}^{\text{max}}, \quad \{\cdot\} = i\_{\text{e}} \text{ g\{g,chp\}} \tag{4}$$

$$P\_{\{\cdot\}}^{\min} \le P\_{\{\cdot\},t}^{\text{DA}} - r\_{\{\cdot\},t'}^{-} \quad \{\cdot\} = i\_{\text{e}} \text{ g\{g,chp\}} \tag{5}$$

$$-r\_{\{\cdot\},t}^{-} \le P\_{\{\cdot\},t}^{w}(\zeta) \le r\_{\{\cdot\},t'}^{+} \quad \{\cdot\} = i\_{\text{e}} \ggcd \tag{6}$$

*P r P r P i gg chp t t tt* (8)

$$(P\_{\{\cdot\},t}^{\text{DA}} + r\_{\{\cdot\},t}^{+}) - (P\_{\{\cdot\},t-1}^{\text{DA}} - r\_{\{\cdot\},t-1}^{-}) \le P\_{\{\cdot\},\prime}^{\text{RU}} \quad \{\cdot\} = i\_{\text{ev}} \text{gg}, \text{chp} \tag{7}$$

$$(P\_{\{\cdot\},t-1}^{\text{DA}} + r\_{\{\cdot\},t-1}^{+}) - (P\_{\{\cdot\},t}^{\text{DA}} - r\_{\{\cdot\},t}^{-}) \le P\_{\{\cdot\},\prime}^{\text{RD}} \cdot \{\cdot\} = i\_{\text{e}} \text{gg}\_{\prime} \text{ch}p \tag{8}$$

DA DA RD

{ }, 1 { }, 1 { }, { }, { } e ( ) ( ) , {} , ,

$$\begin{split} \sum\_{i\_{\ell} \in I\_{\ell}} & \left( P\_{i\_{\ell},t}^{DA} + P\_{i\_{\ell},t}^{RT}(\zeta) \right) + \sum\_{\mathcal{S} \mathcal{S}} \left( P\_{\mathcal{S} \mathcal{S},t}^{DA} + P\_{\mathcal{S} \mathcal{S},t}^{RT}(\zeta) \right) + \left( P\_{\mathcal{S} \mathcal{I},t}^{DA} + P\_{\mathcal{S} \mathcal{I},t}^{RT}(\zeta) \right) \\ & + \sum\_{j \in I\_{j}} \left( \omega\_{j,t}^{DA} + \zeta\_{j,t} \right) + \sum\_{l\_{\ell} \in L\_{\ell}} f\_{l\_{\epsilon,t}}^{DA, jini} = \sum\_{k\_{\ell} \in K\_{\ell}} P\_{k\epsilon,t}^{d} + \sum\_{b\_{\ell} \in B\_{\ell}} P\_{b\_{\ell},t}^{EB} + \sum\_{l\_{\ell} \in L\_{\ell}} f\_{l\_{\ell,t}}^{DA, jini} \\ & \left\{ \begin{array}{c} -f\_{l\_{\ell}}^{\max} \le Q \big( P\_{l\_{\ell},t}^{DA} + P\_{l\_{\ell},t}^{w}(\zeta) \big) + Q^{w} (\omega\_{j,t}^{DA} + \zeta\_{j,t}) - Q^{d} (P\_{k\epsilon,t}^{d} + P\_{b\_{\ell},t}^{EB}) \\ \text{bmatrix} & -\zeta\_{j,t} \end{array} \right\} \end{split} \tag{9}$$

$$\begin{cases} -f\_{l\_{\varepsilon}}^{\max} \le \mathcal{Q}^{\rm g}(P\_{\dot{i}\_{\varepsilon},t}^{\rm DFA} + P\_{\dot{i}\_{\varepsilon},t}^{w}(\zeta)) + \mathcal{Q}^{w}(\omega\_{\dot{j},t}^{\rm DFA} + \zeta\_{\dot{j},t}) - \mathcal{Q}^{u}(P\_{\dot{k}\varepsilon,t}^{\rm e} + P\_{\dot{b}\_{\varepsilon},t}^{\rm e}) \\\ \mathcal{Q}^{\rm g}(P\_{\dot{i}\_{\varepsilon},t}^{\rm DA} + P\_{\dot{i}\_{\varepsilon},t}^{w}(\zeta)) + \mathcal{Q}^{w}(\omega\_{\dot{j},t}^{\rm DA} + \zeta\_{\dot{j},t}) - \mathcal{Q}^{d}(P\_{\dot{k}\varepsilon,t}^{\rm d} + P\_{\dot{b}\_{\varepsilon},t}^{\rm EB}) \le f\_{l\_{\varepsilon}}^{\rm max} \end{cases} \tag{10}$$

However, constraints Equations (6) and (10) may not always be satisfied [27] due to the uncertain wind power forecasting errors, or the strict restrictions may result in high operational costs. Additionally, individual chance constraints may not be satisfied simultaneously with a certain confidence level. To deal with this problem, the two constraints are converted into joint chance constraints in Equations (11) and (12).

$$\mathbb{P}\{-r^{-}\_{\{\cdot\},t} \le P^{\mathbb{RT}}\_{\{\cdot\},t}(\mathbb{Q}) \le r^{+}\_{\{\cdot\},t}\} \ge 1 - \varepsilon^{\text{gen}}, \{\cdot\} = i\_{\text{e}}, \text{gg}, \text{clp} \quad \forall i\_{\text{e}}, \text{gg}, \text{clp} \tag{11}$$

$$\mathbb{P}\{-f\_{l\_{\varepsilon}}^{\max}\leq Q\xi(P\_{i\_{\mathrm{e}},l}^{\mathrm{DA}}+P\_{i\_{\mathrm{e}},l}^{\mathrm{RT}}(\zeta))+Q^{\mathrm{w}}(\omega\_{j,l}^{\mathrm{DA}}+\zeta\_{j,l})-Q^{\mathrm{d}}(P\_{k\_{\mathrm{e}},l}^{\mathrm{d}}+P\_{b\_{\mathrm{e}},l}^{\mathrm{EB}})\leq f\_{l\_{\mathrm{e}}}^{\max}\}\geq 1-\varepsilon^{\mathrm{grid}}\,\forall l\_{\mathrm{e}}\tag{12}$$

#### 2.2.2. Constraints for Natural Gas System

Equation (13) explains the relation between the natural gas pressure of gas compressors' headend nodes and terminals. The nodal natural gas balance is given in Equation (14), whose total demand includes gas load, gas consumed by CHP, and gas-fired units. Constraint Equation (15) implies that gas flows only in the positive/negative direction, which cannot exist simultaneously in the gas pipelines. Equation (16) is the Weymouth gas flow equation [28], which describes the relationship between natural gas pressure and natural gas flow in a steady-state condition. The constraints of natural gas pressure at each node and the limits of the natural gas flow in each gas pipeline are represented by Equations (17) and (18), respectively.

$$
\psi\_{j\_{\mathbf{g}'}t} = \rho\_{\mathbf{c}} \psi\_{i\_{\mathbf{g}'}t} \tag{13}
$$

$$w\_{\mathbf{j}\_{\mathbf{g},t}}^{\text{source}} + \sum\_{\mathbf{i}\_{\mathbf{g}} \mathbf{j}\_{\mathbf{g}} \in \mathbf{Z}(j\_{\mathbf{g}})} w\_{\mathbf{i}\_{\mathbf{g}} \mathbf{j}\_{\mathbf{g}},t} - w\_{\mathbf{j}\_{\mathbf{g}},t}^{\text{ch}p} - w\_{\mathbf{j}\_{\mathbf{g}},t}^{\text{gg}} - w\_{\mathbf{j}\_{\mathbf{g}},t}^{\text{load}} = \sum\_{\mathbf{j}\_{\mathbf{g}} \mathbf{k}\_{\mathbf{g}} \in \mathbf{Z}(k)} w\_{\mathbf{j}\_{\mathbf{g}} \mathbf{k}\_{\mathbf{g}},t} \tag{14}$$

$$w\_{i\_\mathsf{g}j\_\mathsf{g},t} + w\_{j\_\mathsf{g}i\_\mathsf{g},t} = 0 \tag{15}$$

$$w\_{i\_\text{g}j\_\text{g},t} = \mathbb{C}\_{i\_\text{g}j\_\text{g}} \sqrt{\left| \psi\_{i\_\text{g},t}^2 - \psi\_{j\_\text{g},t}^2 \right|} \tag{16}$$

$$
\psi\_{\rm min} \le \psi\_{i\_{\rm \beta},t} \le \psi\_{\rm max} \tag{17}
$$

$$w\_{i\_\mathsf{g}j\_\mathsf{g},\min} \le w\_{i\_\mathsf{g}j\_\mathsf{g},t} \le w\_{i\_\mathsf{g}j\_\mathsf{g},\max} \tag{18}$$

#### 2.2.3. Constraints for Heating System

A closed-cycle heating system is employed, which comprises supply and return pipelines. Hot water is selected as the heat medium for transmission, and quality adjustment is adopted by adjusting the temperature of the heat medium. When the water from different pipelines flows into the same node, there will be a mixture of water described by Equations (19) and (20) [29]. Then, Equations (21) and (22) show that the temperature of the output at each node is equal to the mixed water temperature. Equation (23) presents the relationship between heat demand and the heat flow. Furthermore, Equation (24) takes the loss of heat transmission into consideration. The nodal heat balance is given in Equation (25) [30].

$$\sum\_{b \in S\_j^-} \left( T\_{b,t}^{\text{ps,outb}} \cdot m\_{\text{s}\_b} \right) = T\_{i\_{\text{b}},t}^{\text{ms}} \cdot \sum\_{b \in S\_b^-} m\_{\text{s}\_b} \tag{19}$$

$$\sum\_{b \in S\_{\mathfrak{h}}^{+}} (T\_{b,t}^{\text{pr,outb}} \cdot m\_{r\_{\mathfrak{h}}}) = T\_{\dot{\mathfrak{h}},t}^{\text{mr}} \cdot \sum\_{b \in S\_{\mathfrak{h}}^{+}} m\_{r\_{\mathfrak{h}}} \tag{20}$$

$$T\_{b,t}^{\rm ps,inb} = T\_{i\_{\rm h},t'}^{\rm ms} \, b \in \mathcal{S}\_{i\_{\rm h}}^{+} \tag{21}$$

$$T\_{b,t}^{\text{pr,inb}} = T\_{\dot{\imath}\_{\text{h}}}^{\text{mr}} \; b \in \mathbb{S}\_{\dot{\imath}\_{\text{h}}}^{-} \tag{22}$$

$$d\_{i\_{\rm h},t}^{\rm heat} = C\_{\rm \bf \bf \bf} m\_{s\_{\rm b}} (T\_{i\_{\rm h},t}^{\rm ms} - T\_{i\_{\rm h},t}^{\rm mr}) \tag{23}$$

$$T\_{b,t}^{\text{outb}} = (T\_{b,t}^{\text{inb}} - T^{\text{a}}) \exp(-\frac{\xi^b L^b}{\mathcal{C}\_{\mathbb{P}} m^b}) + T^{\text{a}} \tag{24}$$

$$H\_{\rm chp,t}^{\rm DA} + H\_{\rm EB,t}^{\rm DA} = \sum\_{\dot{\mathbf{h}} \in I\_{\rm h}} d\_{\dot{\mathbf{h}},t}^{\rm heat} + \sum\_{\dot{\mathbf{h}}\_{\rm h} \dot{\mathbf{h}} \in I \{\dot{\mathbf{h}}\_{\rm h}\}} H\_{\rm h,\dot{\mathbf{h}},t}^{\rm loss} \tag{25}$$

#### 2.2.4. Constraints for Multiple Energy Converters

Various energy converters are applied to make use of natural gas and electricity resources conjointly to realize the policy to adopt a balanced energy mix. The constraints Equations (26)–(29) detail the energy input–output relationship energy converters, i.e., CHP, gas-fired generations, and electric boiler.

$$H\_{\rm clp,t}^{\rm DA} = \eta\_{\rm he}^{\rm clp} P\_{\rm clp,t}^{\rm DA} \tag{26}$$

$$P\_{\rm chp,t}^{\rm DA} = \eta\_{\rm ge}^{\rm chp} w\_{i\_{\rm g,t}}^{\rm chp,DA} \tag{27}$$

$$P\_{\text{gg},t}^{\text{DA}} = \eta\_{\text{ge}}^{\text{gg}} w\_{i\_{\text{ge}},t}^{\text{gg,DA}} \tag{28}$$

$$H\_{\rm EB,t}^{\rm DA} = \eta\_{\rm he}^{\rm EB} P\_{b\_{\rm e},t}^{\rm EB,DA} \tag{29}$$

#### *2.3. Formulation of Real-Time Operation*

The real-time stage considers the re-dispatch and adjustive actions to address wind power uncertainty. The second-stage objective function contains two parts: (I) fines for the overrated or underrated schedule in the day-ahead stage and (II) fines for wind power curtailment or load shedding, as shown in Equation (30).

$$\min\_{\begin{subarray}{c}P^{DA}, P^{RT}(\boldsymbol{\zeta}),\\ \tau^{+}, \tau^{-}\end{subarray}} \sum\_{\begin{subarray}{c}i\_{\{\cdot\}}, \{P^{DA}\_{\{\cdot\}}\},\\ i\_{\mathcal{E}} \in I\_{\mathcal{E}}, i\_{\mathcal{E}} \in I\_{\mathcal{E}}\end{subarray}} \lambda\_{\{\cdot\}} |P^{DA}\_{\{\cdot\}} - P^{RT}\_{\{\cdot\}, i}| + \lambda\_{i\_{\mathcal{E}}} \mathcal{G}\_{\text{gs}, i} + \lambda\_{i\_{\mathcal{E}}}^{\text{sled}} \sum\_{\begin{subarray}{c}k\_{\mathcal{E}} \in \mathcal{E}\_{\text{i}} \end{subarray}} l^{\text{sled}}\_{\mathbf{k}\boldsymbol{\varepsilon}} + \lambda\_{i\_{\mathcal{E}}}^{\text{sied}} |\sum\_{\begin{subarray}{c}j \in \mathcal{J} \\ j \in \mathcal{J}\end{subarray}} \omega^{RT}\_{\mathbf{j}, \mathcal{I}} - (\omega^{DA}\_{\boldsymbol{j}, \mathcal{I}} + \boldsymbol{\zeta}\_{\boldsymbol{j}\boldsymbol{\beta}})|\_{\{\cdot\}} \tag{30}$$

Constraints for the real-time operation are the same as the constraints in the day-ahead stage, Equations (7), (8), (10) and (13)–(29), and at the same time include Equations (31)–(34). The constraint Equation (31) requires the load shedding quantity to be no more than actual energy demands. Wind power curtailment quantity is restricted by Equation (32). Constraint Equation (33) limits the adjusted power outputs of traditional generations, gas-fired generations, and CHPs. Equation (34) is the real-time power balance constraint considering wind spills and load shedding.

$$0 \le \sum\_{k\_{\mathbf{c}} \in K\_{\mathbf{c}}} l\_{k\_{\mathbf{c}}, t}^{\text{shell}} \le \sum\_{k\_{\mathbf{c}} \in K\_{\mathbf{c}}} P\_{k\_{\mathbf{c}}, t}^{\text{d}} + \sum\_{b\_{\mathbf{c}} \in B\_{\mathbf{c}}} P\_{b\_{\mathbf{c}}, t}^{\text{EB}} \tag{31}$$

$$0 \le P\_{\rm spi,t}^{\rm RT} \le \sum\_{j \in I} \left( \omega\_{j,t}^{\rm DA} + \zeta\_{j,t} \right) \tag{32}$$

$$P\_{\{\cdot\},t}^{\rm DA} - r\_{\{\cdot\},t}^{-} \le P\_{\{\cdot\},t}^{\rm RT} \le P\_{\{\cdot\},t}^{\rm DA} + r\_{\{\cdot\},t'}^{+} \{\cdot\} = i\_{\rm e} \ggcd \tag{33}$$

$$\begin{aligned} \sum\_{\mathbf{i}\_{\mathbf{e}} \in I\_{\mathbf{e}}} & \left( P\_{\mathbf{i}\_{\mathbf{e}},t}^{\text{DA}} + P\_{\mathbf{i}\_{\mathbf{e}},t}^{\text{RT}} \right) + \sum\_{\mathbf{g}\mathbf{g}} \left( P\_{\mathbf{g}\mathbf{g},t}^{\text{DA}} + P\_{\mathbf{g}\mathbf{g},t}^{\text{RT}} \right) + \left( P\_{\mathbf{chp},t}^{\text{DA}} + P\_{\mathbf{chp},t}^{\text{RT}} \right) + \sum\_{j \in \mathcal{J}} \omega\_{j,t}^{\text{RT}} + \\ \sum\_{\mathbf{i}\_{\mathbf{e}} \in L\_{\mathbf{e}}} f\_{I\_{\mathbf{e},t}}^{\text{DA,ini}} + \sum\_{\mathbf{k}\_{\mathbf{e}} \in K\_{\mathbf{e}}} l\_{\mathbf{k}\_{\mathbf{e}},t}^{\text{shell}} = \sum\_{\mathbf{k}\_{\mathbf{e}} \in K\_{\mathbf{e}}} P\_{\mathbf{k}\_{\mathbf{e}},t}^{\text{d}} + \sum\_{\mathbf{b}\_{\mathbf{e}} \in B\_{\mathbf{e}}} P\_{\mathbf{b}\_{\mathbf{e}},t}^{\text{EB,RT}} + \sum\_{l\_{\mathbf{e}} \in L\_{\mathbf{e}}} f\_{I\_{\mathbf{e},t}}^{\text{DA,ini}} + P\_{\mathbf{s}\mathbf{p},t}^{\text{RT}} \end{aligned} \tag{34}$$

#### **3. Proposed Solution Method**

The proposed two-stage DR-JCCD problem cannot be solved easily due to the wind power uncertainty and the chance constraints. To find a solution to this problem, firstly, a Wasserstein distance-based ambiguity set is used to collect the uncertain wind power distributions information. Then, based on this ambiguity set, the day-ahead objective function is reformulated through linear decision rules by considering the decision variables' ambiguity in the worst-case expectation. Moreover, the joint chance constraints are transformed into tractable constraints through the Bonferroni approximation and the Worst-Conditional Value-at-Risk approximation. Furthermore, the Weymouth gas flow equation of the proposed model will increase computational burden due to its nonlinear and nonconvex nature. Therefore, the linear programming technique is utilized to solve the gas flow equation.

#### *3.1. Basic Formulation*

The proposed model is described as

$$\min\_{\mathbf{x}\in X} \mathbf{c}^{\prime}\mathbf{x} + \sup\_{\mathbf{P}\in D\_{\mathcal{J}}} E^{\mathbf{P}}[Q(\mathbf{x}, \mathcal{J})] \tag{35}$$

*s*.*t*. Ax, < *b* (36)

$$\mathbb{P}\{g(\bar{f})\} \ge 1 - \varepsilon^{\text{grid}} \tag{37}$$

$$\Pr\{h(P\_{\{\cdot\},t}^{\text{RT}}(\zeta))\} \ge 1 - \varepsilon^{\text{gen}} \tag{38}$$

$$\min\_{y} f^{\prime} y$$

$$\text{s.t. } Ex + Fy + G\mathcal{J} \le h \tag{40}$$

The objective function Equation (35) is to minimize the day-ahead operation cost and the expected cost caused by the energy adjustments Equation (1) combined. *x* is the decision including the energy output, reserve capacity, and adjustment of multiple energy devices. *D<sup>ζ</sup>* indicates the ambiguity set containing possible probability distribution P of wind data. Equations (36)–(38) present the day-ahead constraints. The real-time model is shown in Equations (39) and (40), where *f* represents the coefficient of decision variables in Equation (39).

#### *3.2. Wasserstein Distance-Based Ambiguity Set*

In order to estimate the probability distributions of wind power uncertainty, it is important to build an effective ambiguity set. Although potential probability distribution is uncertain, an enormous number of recorded historical data are accessible. Therefore, an empirical distribution P<sup>N</sup> = <sup>1</sup> <sup>N</sup> ∑ N k=1 *δζ* <sup>k</sup> can be considered as the approximate substitution for the true distribution P, where *δ<sup>ζ</sup>* <sup>k</sup> presents the Dirac measure on the wind power forecasting error sample *ζ* k [31]. To estimate the distance between the true distribution P and an empirical distribution PN, the Wasserstein distance is defined as follows.

Definition of the Wasserstein distance [32]: The Wasserstein distance *dw*(P1, P2) : *R <sup>W</sup>* <sup>×</sup> *<sup>R</sup> <sup>W</sup>* <sup>→</sup> *<sup>R</sup>* is defined via

$$d\_{\mathbf{W}}(\mathbf{P}\_1, \mathbf{P}\_2) = \inf \left\{ \int\_{\mathbb{R}^{\mathbf{w}} \times \mathbb{R}^{\mathbf{w}}} \|\zeta\_1 - \zeta\_2\| \prod (d\zeta\_1, d\zeta\_2) \right\} \tag{41}$$

where k*ζ*<sup>1</sup> − *ζ*2k is the distance between random variables *ζ*<sup>1</sup> and *ζ*2. Additionally, 1-norm is applied in this paper due to its superior numerical tractability. *M*(Ξ) denotes all probability measures of wind power uncertainty supported on the polyhedron Ξ = *<sup>ζ</sup>* <sup>∈</sup> *<sup>R</sup>*<sup>W</sup> : *<sup>H</sup><sup>ζ</sup>* <sup>≤</sup> *<sup>h</sup>* . The Wasserstein distance serves to establish an array of ambiguity sets. Every ambiguity set differs from the empirical distribution within the preset distance [33]:

$$D\_{\overline{\zeta}} \triangleq \{ \mathbf{P} \in M(\Xi) : \mathcal{W}(\mathbf{P}, \mathbf{P}\_N) \le \rho \}\tag{42}$$

#### *3.3. Reformulation of Objective Function*

It is complex and time-consuming to directly find the exact solution to DRO problems when the decision variables are coupled with random variables. Therefore, the use of LDRs [34], which is a typical approximate method that can deal with the coupling relationship between decision variables and uncertain parameters, is used to approximate the model [35]. In this context, the objective function Equation (35) can be reformulated as a conic program [32].

$$\begin{split} & \max\_{\mathbf{P} \in D\_{\zeta}} E^{P}[\gamma\_{c}^{T}(Y\_{0} + Y\zeta)] \\ &= \begin{cases} \max\_{\mathbf{P} \in D\_{\zeta}} \int\_{\Sigma} \gamma\_{c}^{T}(Y\_{0} + Y\zeta) \mathbf{P}(\zeta) d\zeta \\\ s.t. \frac{1}{N} \sum\_{i=1}^{N} \int\_{\Sigma} ||\zeta - \overset{\wedge}{\zeta\_{i}}|| \mathbf{P}\_{i}(d\zeta) \le \rho \quad \forall i \le N \\\ \min\_{\boldsymbol{\lambda}^{o}, \rho^{o}, \eta^{o}} \rho^{o} + \frac{1}{N} \sum\_{i=1}^{N} s\_{i}^{o} \\\ \text{s.t. } & \gamma\_{c}^{T}(Y\_{0} + Y\_{\xi i}^{\diamond}) + \gamma\_{i}^{oT}(h - H\_{\xi i}^{\diamond}) \le s\_{i}^{o} \quad \forall i \le N \\\ & \|H^{T} \gamma\_{i}^{o} - (Y\_{0} + Y)^{T} c\|\_{\*} \le \lambda^{o} \quad \forall i \le N \end{cases} \end{split} \tag{43}$$

where *γ o i* , *λ o* , and *s <sup>o</sup>* are auxiliary variables.

#### *3.4. Approximation of Joint Chance Constraints*

Joint chance constraints Equations (11) and (12) include a series of constraints of energy output adjustment and transmission lines separately. We consider the two joint chance constraints in a general form.

$$\mathbb{P}[A\_l \mathbf{x} + (B\_l \mathbf{Y} + \mathbf{C}\_l)\mathbb{Z} \le b\_l \text{ } \forall l \le L] \ge 1 - \varepsilon \tag{44}$$

where *l* represents the index of energy devices or transmission lines and is the total quantity of energy devices or transmission lines. *ε* is a predefined confidence level.

Then, the joint chance constraints can be divided into L individual chance constraints whose confidence level is *ε<sup>l</sup>* = *ε*/*L* by the Bonferroni conservative approximation [36]:

$$\begin{aligned} \min\_{\mathbf{P} \in D\_{\tilde{\zeta}}} & \mathbb{P}[A\_l \mathbf{x} + (B\_l Y + \mathbb{C}\_l) \tilde{\xi} \le b\_l] \ge 1 - \varepsilon\_l \\ &= \min\_{\mathbf{P} \in D\_{\tilde{\zeta}}} & \mathbb{P}[A\_l \mathbf{x} + (B\_l Y + \mathbb{C}\_l) \tilde{\xi} - b\_l \le 0 \mid \ge 1 - \varepsilon\_l \end{aligned} \tag{45}$$

The worst-case Conditional Value-at-Risk approximation can be used to transform Equation (45) [24] into

$$\begin{split} & \max\_{\mathbf{P} \in D\_{\tilde{\zeta}}} \mathbf{P} - \mathsf{CVaR}\_{\ell l} [A\_l \mathbf{x} (B\_l \mathbf{Y} + \mathsf{C}\_l) \boldsymbol{\zeta} - b\_l] \\ &= \max\_{\mathbf{P} \in D\_{\tilde{\zeta}}} \inf\_{\tau\_l} \left\{ \tau\_l + \frac{L}{\varepsilon} \mathrm{E}^{\mathbf{P}} [(A\_l \mathbf{x} (B\_l \mathbf{Y} + \mathsf{C}\_l) \boldsymbol{\zeta} - b\_l - \tau\_l)^{+}] \right\} \\ &= \max\_{\mathbf{P} \in D\_{\tilde{\zeta}}} \inf\_{\tau\_l} \left\{ \mathrm{E}^{\mathbf{P}} [\max \left\{ \tau\_l, \frac{L}{\varepsilon} (A\_l \mathbf{x} (B\_l \mathbf{Y} + \mathsf{C}\_l) \boldsymbol{\zeta} - b\_l) \right\}] \right\} \end{split} \tag{46}$$

which can be rewritten as s.t. ,

which can be rewritten as

$$\begin{array}{ll}\displaystyle\inf\_{\begin{subarray}{c}\mathsf{T}\_{l},\mathsf{A}\_{l},\mathsf{s}\_{l},\gamma\_{l}\\\mathsf{s},\mathsf{t},\ \mathsf{r}\leq s\_{il}\end{subarray}}\lambda\_{l}\rho+N^{-1}\sum\_{i=1}^{N}s\_{il}\\\displaystyle\mathsf{s},\mathsf{t},\ \mathsf{r}\leq s\_{il}\\\displaystyle\frac{L}{\varepsilon}(A\_{l}\mathsf{x}-b\_{l})+\frac{L}{\varepsilon}(B\_{l}Y+\mathsf{C}\_{l})\overset{\scriptstyle{\scriptstyle\rightarrow}}{\zeta}\_{i}+(1-\frac{L}{\varepsilon})\tau\_{l}\\\displaystyle+\gamma\_{il}^{T}(\mathsf{h}-\varPi\_{\xi\_{i}}^{\scriptstyle{\leftarrow}})\leq s\_{il}\\\displaystyle\|H^{T}\gamma\_{il}-\frac{L}{\varepsilon}(B\_{l}Y+\mathsf{C}\_{l})^{T}\|\_{\*}\leq\lambda\_{k}\\\displaystyle\forall i,\ \forall l\end{array}\tag{47}$$

(46)

#### *3.5. Reformulation of the Weymouth Gas Flow Equation 3.5. Reformulation of the Weymouth Gas Flow Equation*

*Energies* **2022**, *15*, x FOR PEER REVIEW 9 of 19

P

*CVaR A x B Y C b*

*N*

max inf [( ( ) ) ]

*l l l l ll <sup>D</sup>*

max inf [max , ( ( ) ) ]

*l ll l l <sup>D</sup>*

 

> 

*<sup>L</sup> E Ax BY C b*

*<sup>L</sup> E Ax BY C b*

*s i l*

max P [ ( ) ]

*ll l l l <sup>D</sup>*

P

1 ,,,

*l il*

*l il <sup>i</sup> <sup>s</sup>*

1

*N s*

P

P

inf

 

*l ll l*

*l*

*l*

P

The Weymouth gas flow Equation (16), applied to characterize the natural gas flow, is nonlinear and nonconvex. These properties make the optimization of natural gas system operation an NP-hard problem. The problem can be solved by mixed-integer linear programming techniques [37]. Among linear programming techniques, piecewise linear functions describing nonlinearities and binary variables can avoid local optima due to nonconvexities. The Weymouth gas flow Equation (16), applied to characterize the natural gas flow, is nonlinear and nonconvex. These properties make the optimization of natural gas system operation an NP‐hard problem. The problem can be solved by mixed‐integer linear pro‐ gramming techniques [37]. Among linear programming techniques, piecewise linearfunc‐ tions describing nonlinearities and binary variables can avoid local optima due to non‐ convexities.

Assuming that natural gas flows from pipe node *i*<sup>g</sup> to pipe node *j*g, the variable *ψ* 2 *ig*,*t* is replaced by the second-order conic *ψ* , *ig*,*t* , representing the pressure of the pipe node *i*g. Then, Equation (13) can be transformed into the constraint as follows: Assuming that natural gas flows from pipe node <sup>g</sup> *i* to pipe node <sup>g</sup>*j* , the variable 2 , *g i t* is replaced by the second‐order conic , , *g i t* , representing the pressure of the pipe node <sup>g</sup> *i* . Then, Equation (13) can be transformed into the constraint as follows:

$$w\_{i\_{\mathfrak{g}}j\_{\mathfrak{g}},t}^2 = \mathbb{C}\_{\text{ij}}^2 (\psi\_{i\_{\mathfrak{g}},t} - \psi\_{j\_{\mathfrak{g}},t}) \tag{48}$$

Next, the square term of each pipeline on the left hand of the equation can be piecewise approximated into m linear segments shown in Figure 2. Next, the square term of each pipeline on the left hand of the equation can be piece‐ wise approximated into m linear segments shown in Figure 2.

**Figure 2.** Approximation of a nonlinear function. **Figure 2.** Approximation of a nonlinear function.

When the direction of the gas flow is from node <sup>g</sup> *i* to node <sup>g</sup> *j* , the range of the gas flow , *g g wij t* is from 0 to ,max *g g wi j* , which is the maximum flow in the pipelines. The trans‐ formed model of <sup>2</sup> , *g g wij t* is shown as: When the direction of the gas flow is from node *i*<sup>g</sup> to node *j*g, the range of the gas flow *wi<sup>g</sup> <sup>j</sup>g*,*<sup>t</sup>* is from 0 to *wi<sup>g</sup> <sup>j</sup>g*,max, which is the maximum flow in the pipelines. The transformed model of *w* 2 *ig jg*,*t* is shown as:

$$\begin{aligned} 0 \le A^l\_{i\_\delta j\_\delta, t} &\le \delta^l\_{i\_\delta j\_\delta, t} d^l\_{i\_\delta j\_\delta} \quad (l=1) \\ \delta^l\_{i\_\delta j\_\delta, t} d^{l-1}\_{i\_\delta j\_\delta} &\le A^l\_{i\_\delta j\_\delta, t} \le \delta^l\_{i\_\delta j\_\delta, t} d^l\_{i\_\delta j\_\delta} \quad (l \ge 2) \\ \sum\_{i=1}^m \delta^l\_{i\_\delta j\_\delta, t} = 1 \quad \delta^l\_{i\_\delta j\_\delta, t} = 0, 1; \quad A^l\_{i\_\delta j\_\delta, t} = \omega\_{i\_\delta j\_\delta, t} \\ f^l\_{i\_\delta j\_\delta} = 0 \quad (l=1); \quad f^l\_{i\_\delta j\_\delta} = (d^{l-1}\_{i\_\delta j\_\delta})^2 \quad (l \ge 2) \\ \boldsymbol{k}^l\_{i\_\delta j\_\delta} = \frac{d^l\_{i\_\delta j\_\delta}}{f^l\_{i\_\delta j\_\delta} - l^{l-1}\_{i\_\delta j\_\delta}} \; l \ge 2 \\ d^l\_{i\_\delta j\_\delta} = \frac{\omega\_{i\_\delta j\_\delta} \max^l}{m} \\ \varphi^2\_{\min} \le \psi\_{i\_\delta, t}^r \le \varphi^2\_{\max} \end{aligned} \tag{49}$$

Finally, the Weymouth gas flow equation is transformed into the form as: Finally, the Weymouth gas flow equation is transformed into the form as:

gg gg gg

0 ( 1)

*A dl*

, , 1

*l ll ij t ij t ij ll l ll ij t ij ij t ij t ij <sup>m</sup> ll l*

g g

*l l i j*

g gg g g

*f l*

,max

*l i j l*

*d*

g

*f*

*i j l i j*

*i j*

g

*i t*

*m*

2, 2 min , max

gg gg gg gg gg

( 2)

*d A dl*

 

1 0,1;

*f l fd l*

2

*l*

, ,,

gg gg gg

1

*l ll ij ij ij*

gg gg gg gg

 

*ij t ij t ij t ij t*

, , ,,

0 ( 1); ( ) ( 2)

1 2

 

(49)

*A*

*Energies* **2022**, *15*, x FOR PEER REVIEW 10 of 19

g g

g g

*i j*

*k*

*d*

1

*i*

$$\mathbf{C}\_{i\_\mathbf{g}j\_\mathbf{g}}^2 (\boldsymbol{\psi}\_{i\_\mathbf{g},t}^\prime - \boldsymbol{\psi}\_{j\_\mathbf{g},t}^\prime) = \sum\_{i=1}^m \left[ (\boldsymbol{A}\_{i\_\mathbf{g}j\_\mathbf{g},t}^l - \boldsymbol{d}\_{i\_\mathbf{g}j\_\mathbf{g}}^{l-1}) \boldsymbol{k}\_{i\_\mathbf{g}j\_\mathbf{g}}^l + \boldsymbol{f}\_{i\_\mathbf{g}j\_\mathbf{g}}^l \boldsymbol{\delta}\_{i\_\mathbf{g}j\_\mathbf{g},t}^l \right] \tag{50}$$

Therefore, the two-stage distributionally robust joint chance-constrained dispatch model considering wind power uncertainty can be converted into a mixed-integer conic reformulation with Equations (1) and (30) as the objective function and Equations (2)–(5), (7)–(9), (11)–(15), (17)–(29), (31)–(34), (43), (47), and (50) as the constraints, which can be solved directly by calling Gurobi solver under Matlab. Therefore, the two‐stage distributionally robust joint chance‐constrained dispatch model considering wind power uncertainty can be converted into a mixed‐integer conic reformulation with Equations (1) and (30) as the objective function and Equations (2)–(5), (7)–(9), (11)–(15), (17)–(29), (31)–(34), (43), (47), and (50) as the constraints, which can be solved directly by calling Gurobi solver under Matlab.

#### **4. Case Study 4. Case Study**

As shown in Figure 3, a regionally integrated energy system for electricity, gas, and heating comprising a modified IEEE 24-bus system, 20 natural gas nodes, and 6 heat nodes, was used to test the validity of the proposed DR-JCCD model. The detailed information for the generators and natural gas sources are given in Appendix B. The power and the gas subsystems have three coupled points: buses 13 and 23 in the power subsystem are severally connected with natural gas nodes 3 and 19 via two gas-fired generators, and a CHP at bus 22 is linked with gas node 6. The heat demand is satisfied by the gas turbine at power bus 22, and an electric boiler supplied by power bus 13. As shown in Figure 3, a regionally integrated energy system for electricity, gas, and heating comprising a modified IEEE 24‐bus system, 20 natural gas nodes, and 6 heat nodes, was used to test the validity of the proposed DR‐JCCD model. The detailed infor‐ mation for the generators and natural gas sources are given in Appendix B. The power and the gas subsystems have three coupled points: buses 13 and 23 in the power subsys‐ tem are severally connected with natural gas nodes 3 and 19 via two gas‐fired generators, and a CHP at bus 22 is linked with gas node 6. The heat demand is satisfied by the gas turbine at power bus 22, and an electric boiler supplied by power bus 13.

**Figure 3.** Modified 24‐electricity bus, 20‐natural gas node, and 6‐heat node system for regional IES. **Figure 3.** Modified 24-electricity bus, 20-natural gas node, and 6-heat node system for regional IES.

#### *4.1. Robust Performance with Different Sample Sizes 4.1. Robust Performance with Different Sample Sizes*

In general, by using another dataset diverse from the experimental one, out-of-sample performance is a helpful tool for assessing the robustness of the optimal schedules [37]. Additionally, sampling errors that comes from limited historical data may cause poor out-of-sample performance in operation. In this condition, the empirical evidence based on out-of-sample forecasting errors is used to measure the robustness of the model in this paper.

As illustrated in Figure 4, an unwise decision that ignores ambiguity (by setting *<sup>ρ</sup>* <sup>=</sup> 0) has a large out-of-sample size, which is the maximum cost (\$5.37 <sup>×</sup> <sup>10</sup><sup>7</sup> ). Therefore, this unwise decision is costlier than a more advanced decision that takes the ambiguity of uncertain wind power into account by setting an appropriate distance *ρ*. The largest difference in cost is up to \$5.1 <sup>×</sup> <sup>10</sup><sup>5</sup> . In short, the distance *ρ* precisely regulates the conservativeness of the optimal decision. A considerable distance will make optimal decisions more independent of the characteristics of the historical data and offer stronger robustness to energy adjustments and reserve policies.

paper.

uncertain wind power into account by setting an appropriate distance

ference in cost is up to \$ <sup>5</sup> 5 1 10 . . In short, the distance

to energy adjustments and reserve policies.

**Figure 4.** The acquisition cost of out‐of‐sample size influenced by training sample sizes and Was‐ serstein distance. **Figure 4.** The acquisition cost of out-of-sample size influenced by training sample sizes and Wasserstein distance.

In general, by using another dataset diverse from the experimental one, out‐of‐sam‐ ple performance is a helpful tool for assessing the robustness of the optimal schedules [37]. Additionally, sampling errors that comes from limited historical data may cause poor out‐of‐sample performance in operation. In this condition, the empirical evidence based on out‐of‐sample forecasting errors is used to measure the robustness of the model in this

As illustrated in Figure 4, an unwise decision that ignores ambiguity (by setting

precisely regulates the conserv‐

. The largest dif‐

 0 ) has a large out‐of‐sample size, which is the maximum cost (\$ <sup>7</sup> 5 37 10 . ). Therefore, this unwise decision is costlier than a more advanced decision that takes the ambiguity of

ativeness of the optimal decision. A considerable distance will make optimal decisions more independent of the characteristics of the historical data and offer strongerrobustness

Furthermore, a model considering uncertain variable ambiguity is more competitive for larger sample sizes N. For example, as shown in Figure 4, the acquisition cost of out‐ of‐sample sizes is higher when using an N = 20 training sample (\$ <sup>7</sup> 5 35 10 . ) versus an N = 100 training sample (\$ <sup>7</sup> 5 32 10 . ). This is due to the fact that a smaller sample size results in a poorer robustness choice in the first stage, necessitating a higher cost for adjustment at the second stage. It is now obvious that the out‐of‐sample acquisition cost decreases gradually as the training sample size increases. It further implies that, with adequate data support, the proposed solution is fairly robust to wind uncertainty. Furthermore, a model considering uncertain variable ambiguity is more competitive for larger sample sizes N. For example, as shown in Figure 4, the acquisition cost of out-ofsample sizes is higher when using an N = 20 training sample (\$5.35 <sup>×</sup> <sup>10</sup><sup>7</sup> ) versus an N = 100 training sample (\$5.32 <sup>×</sup> <sup>10</sup><sup>7</sup> ). This is due to the fact that a smaller sample size results in a poorer robustness choice in the first stage, necessitating a higher cost for adjustment at the second stage. It is now obvious that the out-of-sample acquisition cost decreases gradually as the training sample size increases. It further implies that, with adequate data support, the proposed solution is fairly robust to wind uncertainty.

#### *4.2. The Influence of Different Confidence Levels*

*4.2. The Influence of Different Confidence Levels* Because risk levels of renewable energy uncertainty can be estimated by confidence level or the parameter , we investigate the influence of a series of confidence levels on the total cost. The whole system cost with different parameters 1 , including 99%, 95%, 90%, 75%, 65%, 55%, and 45%, are summarized in Figure 5. Because risk levels of renewable energy uncertainty can be estimated by confidence level or the parameter *ε*, we investigate the influence of a series of confidence levels on the total cost. The whole system cost with different parameters 1 − *ε*, including 99%, 95%, 90%, 75%, 65%, 55%, and 45%, are summarized in Figure 5. *Energies* **2022**, *15*, x FOR PEER REVIEW 12 of 19

**Figure 5.** The impact of different confidence levels on the operational cost. **Figure 5.** The impact of different confidence levels on the operational cost.

With the same , when the increases from 0.01 to 0.55, the marginal cost gradu‐ ally decreases from \$ <sup>7</sup> 5 39 10 . to \$ <sup>7</sup> 5 32 10 . , as shown in Figure 5. The largest marginal cost difference is \$ <sup>5</sup> 7 10 among the listed confidence levels. This is because the lower the decision‐risk maker's tolerance, the greater the necessity for the reserve to balance the unpredictability of wind turbine output, and the greater the amount of natural gas con‐ sumed. Meanwhile, a small value represents a low tolerance level of risk‐taking. It means that as the level of confidence diminishes, the marginal cost decreases. In particu‐ lar, when the confidence level change from 99% to 95%, the marginal cost dramatically decreases by \$ <sup>5</sup> 5 10 , which accounts for 71% of the largest marginal cost difference. Therefore, it will cost more to achieve more a reliable operation of the system. It acts as a reminder to decision‐makers that in practice, they should choose an adequate confidence With the same *ρ*, when the *ε* increases from 0.01 to 0.55, the marginal cost gradually decreases from \$5.39 <sup>×</sup> <sup>10</sup><sup>7</sup> to \$5.32 <sup>×</sup> <sup>10</sup><sup>7</sup> , as shown in Figure 5. The largest marginal cost difference is \$7 <sup>×</sup> <sup>10</sup><sup>5</sup> among the listed confidence levels. This is because the lower the decision-risk maker's tolerance, the greater the necessity for the reserve to balance the unpredictability of wind turbine output, and the greater the amount of natural gas consumed. Meanwhile, a small *ε* value represents a low tolerance level of risk-taking. It means that as the level of confidence diminishes, the marginal cost decreases. In particular, when the confidence level change from 99% to 95%, the marginal cost dramatically decreases by \$5 <sup>×</sup> <sup>10</sup><sup>5</sup> , which accounts for 71% of the largest marginal cost difference. Therefore, it will cost more to achieve more a reliable operation of the system. It acts as a reminder to decision-makers that in practice, they should choose an adequate confidence level to avoid the significant costs associated with high-reliability standards.

To assess the effectiveness of the proposed DR‐JCCD in balancing robustness and economy, this subsection compares the operating costs among the DR‐JCCD approach,

As shown in Figure 6, the total expected costs of the DR‐JCCD model increase as the confidence level increases. The expected costs of the DR‐JCCD model are always between the RO and SP, regardless of how the expenses change. Meanwhile, the RO approach has a higher total expected cost than its counterparts by at least \$ <sup>4</sup> 1 10 and at most \$ <sup>5</sup> 6 4 10 . due to overly conservative decisions on energy reserve and dispatch. In particular, when the Wasserstein distance is zero, the total expected cost of DR‐JCCD is equal to the antic‐ ipated expenses of SP. According to the definition of the Wasserstein distance, the ambi‐ guity set contains only empirical probability distributions of wind power on this condi‐ tion. When the Wasserstein radius approaches infinity and the confidence level reaches 100%, DR‐JCCD almost degrades into RO, where the cost of DR‐JCCD is only 0.02% lower than the cost of RO. This is because the former has the tendency to contain all probability

> **3 10**

, RO, and SP.

**5.32 0.5**

**65% 0.3**

**RO DR‐JCCD SP**

**Figure 6.** Comparison of total expected costs among DR‐JCCD with different

**75% 90% 0.2 95% 99% <sup>0</sup>**

**45% 0.4**

RO, and SP methods.

distributions.

**<sup>10</sup><sup>7</sup>**

**5.33**

**55%**

**5.34 5.35**

**Total expected cost**

**5.36 5.37**

**(\$)**

**5.38 5.39**

level to avoid the significant costs associated with high‐reliability standards.

#### *4.3. Comparisons among DR-JCCD, RO, and SP 4.3. Comparisons among DR‐JCCD, RO, and SP* To assess the effectiveness of the proposed DR‐JCCD in balancing robustness and

**0.4 0.3 0.2 0.1 <sup>0</sup>**

, when the

**<sup>3</sup> 45% <sup>10</sup>**

**Figure 5.** The impact of different confidence levels on the operational cost.

With the same

sumed. Meanwhile, a small

*Energies* **2022**, *15*, x FOR PEER REVIEW 12 of 19

**10<sup>7</sup>**

**Total expected cost (\$)**

**5.32 5.33 5.34 5.35 5.36 5.37 5.38 5.39**

To assess the effectiveness of the proposed DR-JCCD in balancing robustness and economy, this subsection compares the operating costs among the DR-JCCD approach, RO, and SP methods. economy, this subsection compares the operating costs among the DR‐JCCD approach, RO, and SP methods. As shown in Figure 6, the total expected costs of the DR‐JCCD model increase as the

ally decreases from \$ <sup>7</sup> 5 39 10 . to \$ <sup>7</sup> 5 32 10 . , as shown in Figure 5. The largest marginal cost difference is \$ <sup>5</sup> 7 10 among the listed confidence levels. This is because the lower the decision‐risk maker's tolerance, the greater the necessity for the reserve to balance the unpredictability of wind turbine output, and the greater the amount of natural gas con‐

means that as the level of confidence diminishes, the marginal cost decreases. In particu‐ lar, when the confidence level change from 99% to 95%, the marginal cost dramatically decreases by \$ <sup>5</sup> 5 10 , which accounts for 71% of the largest marginal cost difference. Therefore, it will cost more to achieve more a reliable operation of the system. It acts as a reminder to decision‐makers that in practice, they should choose an adequate confidence

level to avoid the significant costs associated with high‐reliability standards.

**99% 95% 90% 75% 65% 55%**

increases from 0.01 to 0.55, the marginal cost gradu‐

value represents a low tolerance level of risk‐taking. It

As shown in Figure 6, the total expected costs of the DR-JCCD model increase as the confidence level increases. The expected costs of the DR-JCCD model are always between the RO and SP, regardless of how the expenses change. Meanwhile, the RO approach has a higher total expected cost than its counterparts by at least \$1 <sup>×</sup> <sup>10</sup><sup>4</sup> and at most \$6.4 <sup>×</sup> <sup>10</sup><sup>5</sup> due to overly conservative decisions on energy reserve and dispatch. In particular, when the Wasserstein distance is zero, the total expected cost of DR-JCCD is equal to the anticipated expenses of SP. According to the definition of the Wasserstein distance, the ambiguity set contains only empirical probability distributions of wind power on this condition. When the Wasserstein radius approaches infinity and the confidence level reaches 100%, DR-JCCD almost degrades into RO, where the cost of DR-JCCD is only 0.02% lower than the cost of RO. This is because the former has the tendency to contain all probability distributions. confidence level increases. The expected costs of the DR‐JCCD model are always between the RO and SP, regardless of how the expenses change. Meanwhile, the RO approach has a higher total expected cost than its counterparts by at least \$ <sup>4</sup> 1 10 and at most \$ <sup>5</sup> 6 4 10 . due to overly conservative decisions on energy reserve and dispatch. In particular, when the Wasserstein distance is zero, the total expected cost of DR‐JCCD is equal to the antic‐ ipated expenses of SP. According to the definition of the Wasserstein distance, the ambi‐ guity set contains only empirical probability distributions of wind power on this condi‐ tion. When the Wasserstein radius approaches infinity and the confidence level reaches 100%, DR‐JCCD almost degrades into RO, where the cost of DR‐JCCD is only 0.02% lower than the cost of RO. This is because the former has the tendency to contain all probability distributions.

**Figure Figure 6. 6.** Comparison Comparison of total expected costs among DR-JCCD with different of total expected costs among DR‐JCCD with different , RO, and SP. *ε*, RO, and SP.

Moreover, the acquisition cost of out-of-sample size of the DR-JCCD model decreases with the increase in training sample size. However, this cost is always between RO and SP, as illustrated in Figure 7. Moreover, the acquisition cost of out‐of‐sample size of the DR‐JCCD model decreases with the increase in training sample size. However, this cost is always between RO and SP, as illustrated in Figure 7.

**Figure 7.** Comparison acquisition cost of out‐of‐sample size among DR‐JCCD with different sizes of training samples, RO, and SP. **Figure 7.** Comparison acquisition cost of out-of-sample size among DR-JCCD with different sizes of training samples, RO, and SP.

Moreover, the acquisition cost of out‐of‐sample size of the DR‐JCCD model decreases with the increasing training sample size. However, this cost is always in‐between RO and SP, as illustrated in Figure 7. Moreover, the acquisition cost of out-of-sample size of the DR-JCCD model decreases with the increasing training sample size. However, this cost is always in-between RO and SP, as illustrated in Figure 7.

From a different perspective, RO, which dispatches and reserves more energy in the worst case of wind power generation, has the most robustness of the three techniques. However, taking accurate dispatch and reserve into account, the proposed approach has lower operational costs than RO. This is due to the fact that it bases its decisions on the worst‐case probability distribution of wind power generation, which means more infor‐ From a different perspective, RO, which dispatches and reserves more energy in the worst case of wind power generation, has the most robustness of the three techniques. However, taking accurate dispatch and reserve into account, the proposed approach has lower operational costs than RO. This is due to the fact that it bases its decisions on the worst-case probability distribution of wind power generation, which means more

mation on uncertain wind power considered. It also has a higher energy reserve than the

Different confidence levels will result in various solutions in multi‐energy manage‐

generations and CHP in the gas system. Note that positive values in Figure 8 represent the result of day‐ahead operation, while negative values represent the result of real‐time

is chose to show the energy mutual assistance effect of gas‐fired

**Time**

A small amount of natural gas is transformed into power at 1:00–3:00 and 24:00 ow‐ ing to the large output of wind power, as illustrated in Figure 8. At this time, the results

**1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Day‐ahead GasGenNeed GasGenNeed Adjustment**

**Figure 8.** Natural gas consumed in day‐head and real‐time operation.

SP to cope with the wind power uncertainty.

ment. Here, 1 95%

dispatch.

**Gas(Mm3)**

**0 0.1 0.2 0.3 0.4**

**0 4**. **0 3**. **0 2**. **0 1**. **Real‐time GasGenNeed**

*4.4. Analysis of Energy Conversions in Energy Balance*

information on uncertain wind power considered. It also has a higher energy reserve than the SP to cope with the wind power uncertainty. mation on uncertain wind power considered. It also has a higher energy reserve than the SP to cope with the wind power uncertainty.

Moreover, the acquisition cost of out‐of‐sample size of the DR‐JCCD model decreases with the increase in training sample size. However, this cost is always between RO and

**1 1.5 2 2.5 3 3.5 4**

**Figure 7.** Comparison acquisition cost of out‐of‐sample size among DR‐JCCD with different sizes of

Moreover, the acquisition cost of out‐of‐sample size of the DR‐JCCD model decreases with the increasing training sample size. However, this cost is always in‐between RO and

From a different perspective, RO, which dispatches and reserves more energy in the worst case of wind power generation, has the most robustness of the three techniques. However, taking accurate dispatch and reserve into account, the proposed approach has lower operational costs than RO. This is due to the fact that it bases its decisions on the worst‐case probability distribution of wind power generation, which means more infor‐

**SP RO DRO‐1 DRO‐2 DRO‐3 DRO‐4**

To sum up, by employing partial distributional information, the DR-JCCD approach realizes that the robustness of the model is well-balanced with its economy. To sum up, by employing partial distributional information, the DR‐JCCD approach realizes that the robustness of the model is well‐balanced with its economy.

#### *4.4. Analysis of Energy Conversions in Energy Balance 4.4. Analysis of Energy Conversions in Energy Balance*

*Energies* **2022**, *15*, x FOR PEER REVIEW 13 of 19

**0 50 100 150 200 250 300 size of samples**

**106**

**1.16 1.165 1.17 1.175 1.18 1.185 1.19 1.195 1.2 1.205**

SP, as illustrated in Figure 7.

**<sup>0</sup> 50 100 150 200 250 300 5.31**

training samples, RO, and SP.

SP, as illustrated in Figure 7.

**Total expected cost (\$)**

**<sup>6</sup> <sup>10</sup><sup>7</sup>**

**5.32 5.33 5.34 5.35 5.36 5.37 <sup>10</sup> <sup>7</sup>**

Different confidence levels will result in various solutions in multi-energy management. Here, 1 − *ε* = 95% is chose to show the energy mutual assistance effect of gasfired generations and CHP in the gas system. Note that positive values in Figure 8 represent the result of day-ahead operation, while negative values represent the result of real-time dispatch. Different confidence levels will result in various solutions in multi‐energy manage‐ ment. Here, 1 95% is chose to show the energy mutual assistance effect of gas‐fired generations and CHP in the gas system. Note that positive values in Figure 8 represent the result of day‐ahead operation, while negative values represent the result of real‐time dispatch.

**Figure 8.** Natural gas consumed in day‐head and real‐time operation. **Figure 8.** Natural gas consumed in day-head and real-time operation.

A small amount of natural gas is transformed into power at 1:00–3:00 and 24:00 ow‐ ing to the large output of wind power, as illustrated in Figure 8. At this time, the results A small amount of natural gas is transformed into power at 1:00–3:00 and 24:00 owing to the large output of wind power, as illustrated in Figure 8. At this time, the results of the day-ahead operation will not be adjusted in real-time. Because the influence of the expectation of wind power real-time deviation at the first stage is considered, the adjustment may always take a cut action to revise the results over day-ahead dispatch between 4:00 and 23:00, where the maximum adjustment amount accounts for 19.83% of the day-ahead dispatch. It can reduce unnecessary costs and verify the robustness of the decisions made at the first stage. Meanwhile, by CHP and gas-fired units, the natural gas system effectively supports the power grid under wind power uncertainty.

The output of the heating system is shown in Figure 9. Additionally, the difference between the heating source and heating demand is precisely the transmission loss, as shown in Figure 9a,b. According to Figure 9c, the CHP and electric boiler convert the electricity in order to meet the heating demand, where CHP is the dominant heating source with 83% of the heating capacity and the electric boiler assists with 17% of the heating capacity. When the wind power is abundant at 1:00–5:00 and 23:00–24:00, the electric boiler transforms more wind power to heat, enhancing the IES capacity of wind power utilization in IES.

#### *4.5. Effect of Gas-Fired Generations on Electric Peak Shaving*

Peak shaving in a multi-energy system proactively adjusts actions of energy utilization to broaden energy sources or reduce short-term multi-energy demand at peak periods. As a device of multi-energy cooperation, gas-fired generations have a positive effect on Electric Peak shaving, as illustrated in Figure 10. At the peak of power consumption, gasfired generation systems quickly adjust their output to meet the power demand, thereby reducing the regulatory burden of the power system. Especially at 11:00, the supply of gas-fired generators reaches its maximum, accounting for 22.26% of the power demand. This demonstrates the advantages of multi-energy cooperation in the proposed model. This is particularly the case when the regulation resource is limited or the regulation cost is too high in a subsystem.

effectively supports the power grid under wind power uncertainty.

**Figure 9.** The output of the heating system. (**a**) The heat balance of the heating system. (**b**) The trans‐ mission loss of heating system. (**c**) The output of heat sources. **Figure 9.** The output of the heating system. (**a**) The heat balance of the heating system. (**b**) The transmission loss of heating system. (**c**) The output of heat sources. *Energies* **2022**, *15*, x FOR PEER REVIEW 15 of 19

of the day‐ahead operation will not be adjusted in real‐time. Because the influence of the expectation of wind power real‐time deviation at the first stage is considered, the adjust‐ ment may always take a cut action to revise the results over day‐ahead dispatch between 4:00 and 23:00, where the maximum adjustment amount accounts for 19.83% of the day‐ ahead dispatch. It can reduce unnecessary costs and verify the robustness of the decisions made at the first stage. Meanwhile, by CHP and gas‐fired units, the natural gas system

The output of the heating system is shown in Figure 9. Additionally, the difference between the heating source and heating demand is precisely the transmission loss, as shown in Figure 9a,b. According to Figure 9c, the CHP and electric boiler convert the elec‐ tricity in order to meet the heating demand, where CHP is the dominant heating source with 83% of the heating capacity and the electric boiler assists with 17% of the heating capacity. When the wind power is abundant at 1:00–5:00 and 23:00–24:00, the electric boiler transforms more wind power to heat, enhancing the IES capacity of wind power

tion cost is too high in a subsystem. **Figure 10.** Actions of gas‐fired generations in Electric Peak‐shaving. **Figure 10.**Actions of gas-fired generations in Electric Peak-shaving.

#### **5. Conclusions 5. Conclusions**

utilization in IES.

A two‐stage distributionally robust joint chance‐constrained dispatch model for elec‐ tricity–gas–heat IES with wind power uncertainty is investigated in this paper. The wind power generation uncertainty is captured in the model by employing the worst‐case prob‐ ability distributional information in an ambiguity set based on Wasserstein distance. In light of the operational risk caused by wind power uncertainty, the joint chance con‐ straints ensure that multiple safety conditions are met simultaneously with a high confi‐ dence level. Next, by the linear decision rules and linear incremental method, the problem is reformulated as a mixed‐integer tractable optimization issue. The effectiveness of the proposed model is corroborated on an electricity–gas–heating regional integrated energy system with a modified IEEE 24‐bus system, with 20 natural gas nodes and 6 heat nodes. Notably, the proposed DR‐JCCD method can pay 1.3% less than the RO method and achieve a more robust out‐of‐sample performance than the SP approach at risk, which is a shaving of 22.3% on the acquisition cost of out‐of‐sample size. Hence, the proposed model achieves a good balance between economy and robustness. Furthermore, the higher the risk preference of the decision‐maker, the cheaper the operating cost of the optimization solutions will be, but this will also lessen the robustness of this scheme. To A two-stage distributionally robust joint chance-constrained dispatch model for electricity–gas–heat IES with wind power uncertainty is investigated in this paper. The wind power generation uncertainty is captured in the model by employing the worst-case probability distributional information in an ambiguity set based on Wasserstein distance. In light of the operational risk caused by wind power uncertainty, the joint chance constraints ensure that multiple safety conditions are met simultaneously with a high confidence level. Next, by the linear decision rules and linear incremental method, the problem is reformulated as a mixed-integer tractable optimization issue. The effectiveness of the proposed model is corroborated on an electricity–gas–heating regional integrated energy system with a modified IEEE 24-bus system, with 20 natural gas nodes and 6 heat nodes. Notably, the proposed DR-JCCD method can pay 1.3% less than the RO method and achieve a more robust out-of-sample performance than the SP approach at risk, which is a shaving of 22.3% on the acquisition cost of out-of-sample size. Hence, the proposed model achieves a good balance between economy and robustness. Furthermore, the higher the risk preference of the decision-maker, the cheaper the operating cost of the optimization solutions will be, but this will also lessen the robustness of this scheme. To put it differently, the DR-JCCD method can provide decision-makers with information on cost and risk.

put it differently, the DR‐JCCD method can provide decision‐makers with information on cost and risk. With the increase in the number of wind power data, the statistical wind power char‐ acteristics are closer to the true wind power distribution. However, a large number of data With the increase in the number of wind power data, the statistical wind power characteristics are closer to the true wind power distribution. However, a large number of data sets will bring about a calculation issue. Therefore, it is necessary to consider effective

sets will bring about a calculation issue. Therefore, it is necessary to consider effective

multi‐energy demand response. In addition, multiple uncertainties can be considered in

**Author Contributions:** Conceptualization, H.L.; methodology, H.L. and Z.F.; software, Z.F.; valida‐ tion, H.L., H.X., N.W.; formal analysis, H.X.; investigation, Z.F.; resources, Z.F. and H.X.; data cura‐ tion, Z.F.; writing—original draft preparation, Z.F. and H.X.; writing—review and editing, H.X.; visualization, H.L. and N.W.; supervision, H.L. and N.W. All authors have read and agreed to the

**Funding:** This research was funded in part by the National Key Research and Development Pro‐ gram of China (Grant No. 2019YFE0118000), in part by the National Natural Science Foundation of China (Grant No. 61963003), and in part by the Natural Science Foundation for Distinguished Young

**Informed Consent Statement:** Written informed consent has been obtained from the patient(s) to

future work, such as various uncertain energy supply and energy demand.

published version of the manuscript.

publish this paper.

Scholars of Guangxi (Grant No. 2018GXNSFFA281006). **Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

scene reduction technology in future work. For more flexible energy management, the IES can take more functional interdependent coupling devices into account and integrate a multi-energy demand response. In addition, multiple uncertainties can be considered in future work, such as various uncertain energy supply and energy demand.

**Author Contributions:** Conceptualization, H.L.; methodology, H.L. and Z.F.; software, Z.F.; validation, H.L., H.X., N.W.; formal analysis, H.X.; investigation, Z.F.; resources, Z.F. and H.X.; data curation, Z.F.; writing—original draft preparation, Z.F. and H.X.; writing—review and editing, H.X.; visualization, H.L. and N.W.; supervision, H.L. and N.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded in part by the National Key Research and Development Program of China (Grant No. 2019YFE0118000), in part by the National Natural Science Foundation of China (Grant No. 61963003), and in part by the Natural Science Foundation for Distinguished Young Scholars of Guangxi (Grant No. 2018GXNSFFA281006).

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

**Informed Consent Statement:** Written informed consent has been obtained from the patient(s) to publish this paper.

**Data Availability Statement:** Not applicable.

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

#### **Appendix A**

**Table A1.** Nomenclature.



#### **Table A1.** *Cont.*

#### **Appendix B**

**Table A2.** Parameters of traditional units, gas-fired units and CHPs.


Where type 0\1\2 represents traditional units\gas-fired units\CHPs.


**Table A3.** Parameters of natural gas sources.

#### **References**

