An Exergy-Enhanced Improved IGDT-Based Optimal Scheduling Model for Electricity–Hydrogen Urban Integrated Energy Systems

Article Content

Abstract

Urban integrated energy systems (UIESs) play a critical role in facilitating low-carbon and high-efficiency energy transitions. However, existing scheduling strategies predominantly focus on energy quantity and cost, often neglecting the heterogeneity of energy quality across electricity, heat, gas, and hydrogen. This paper presents an exergy-enhanced stochastic optimization framework for the optimal scheduling of electricity–hydrogen urban integrated energy systems (EHUIESs) under multiple uncertainties. By incorporating exergy efficiency evaluation into a Stochastic Optimization–Improved Information Gap Decision Theory (SOI-IGDT) framework, the model dynamically balances economic cost with thermodynamic performance. A penalty-based iterative mechanism is introduced to track exergy deviations and guide the system toward higher energy quality. The proposed approach accounts for uncertainties in renewable output, load variation, and Hydrogen-enriched compressed natural gas (HCNG) combustion. Case studies based on a 186-bus UIES coupled with a 20-node HCNG network show that the method improves exergy efficiency by up to 2.18% while maintaining cost robustness across varying confidence levels. These results underscore the significance of integrating exergy into real-time robust optimization for resilient and high-quality energy scheduling.

Keywords:

exergy efficiency; stochastic optimization; SOI-IGDT; multiple uncertainties; integrated energy systems

1. Introduction

Integrated energy systems (IESs) are rapidly evolving to meet the increasing demands of urban decarbonization, the electrification of end-use sectors, and the growing need for multi-energy coordination and interconnection [1]. By coupling various energy carriers such as electricity, heating, cooling, natural gas, and hydrogen, these systems offer a promising pathway to enhance overall energy efficiency and system flexibility [2].

However, conventional IES scheduling models are primarily designed to optimize energy “quantity” or minimize operational costs [3,4]. While effective from an economic standpoint, these models largely overlook the inherent differences in thermodynamic quality among energy carriers—such as the high-grade nature of electricity compared to low-grade thermal energy, or the variable quality of hydrogen-rich gases [5,6]. Treating diverse forms of energy as equivalent in optimization models, despite their distinct thermodynamic potentials, may lead to misleading conclusions. Relying solely on quantity-based performance indicators is insufficient to capture the true effectiveness and rationality of energy utilization in such complex systems [7,8].

To address this gap, the concept of exergy has been increasingly adopted in system analysis. Exergy is defined as the maximum useful work obtainable from a given form of energy when brought into equilibrium with its environment. It provides a unified and physically grounded metric for assessing the quality and usability of energy across different carriers [9]. By quantifying the inherent irreversibilities in energy conversion, transport, and storage processes, exergy-based models can reveal hidden inefficiencies that are often obscured in energy-only assessments [10,11,12]. This enables deeper insights into where energy quality losses occur and how they can be minimized [13,14]. Despite its analytical value, exergy has mostly been used retrospectively as an evaluation tool in existing studies, rather than as a forward-looking input for real-time operational decision-making [15,16]. This separation between evaluation and control limits the ability of IES to adaptively optimize operations for thermodynamic performance, especially under dynamic and uncertain conditions [17,18,19].

In practice, the operation of IES is subject to various uncertainties, including the stochastic fluctuations of renewable generation (e.g., wind and solar), time-varying energy demand, and the complex and unstable combustion characteristics of hydrogen-rich natural gas in gas turbines. These uncertainties pose significant challenges for effective real-time scheduling [20,21]. Moreover, several recent studies have explored uncertainty-aware scheduling in integrated energy systems. For example, Zhang et al. developed a fuzzy-stochastic dispatch model to manage the uncertainty of wind power and load under a carbon trading mechanism [22]. Another study proposed a robust chance-constrained model to optimize supply-demand scheduling in power markets, ensuring reserve requirements under confidence-level constraints [23]. More recently, a hybrid stochastic and chance-constrained bidding strategy was introduced [24], incorporating distributed renewable generation and multi-energy load uncertainties into a day-ahead electricity market framework. While these methods enhance robustness under uncertain conditions, they primarily focus on economic objectives and market coordination. Few of them consider the coupling effects of multiple uncertainty sources, such as renewables and fuel characteristics, and none of them have integrated thermodynamic performance indicators such as exergy efficiency into the scheduling framework. This underscores the need for more adaptive and exergy-aware optimization frameworks under compound uncertainty environments.

While the Information Gap Decision Theory (IGDT) has emerged as a powerful tool for decision-making under severe uncertainty, offering a structured risk-averse optimization framework, traditional implementations have focused exclusively on the degradation of economic objectives [25]. They fail to account for exergy-related performance metrics, which are essential for ensuring high-quality energy utilization in IES [26]. In view of the current research status, IGDT still has some limitations. First, the description of uncertain variables is often oversimplified and fails to capture their inherent randomness. Second, the deviation factor must be manually preset, which introduces subjectivity into the modeling process. Third, conventional IGDT models primarily focus on economic robustness while neglecting thermodynamic characteristics, making them less effective in exergy-sensitive systems such as EHUIES. Therefore, there is a need to improve the IGDT framework by integrating thermodynamic performance indicators and more accurate representations of uncertainty.

To overcome these challenges, this paper proposes a novel optimization framework that explicitly incorporates exergy efficiency into the scheduling process. Termed Stochastic Optimization–Improved IGDT (SOI-IGDT), the framework integrates exergy indicators into a dynamic, uncertainty-aware optimization model. Specifically, an enhanced IGDT method with exergy-oriented features is developed for constructing an optimal scheduling model for electricity–hydrogen urban integrated energy systems (EHUIESs) under multiple uncertainties. In this model, cost robustness and exergy efficiency tracking are jointly managed through a deviation penalty mechanism, which is iteratively updated based on the system’s fluctuating states. This ensures that the system not only withstands external uncertainties but also maintains high thermodynamic performance across the entire scheduling horizon.

A case study based on a real-world urban electricity–hydrogen–natural gas integrated energy system is conducted to validate the proposed SOI-IGDT approach. Simulation results demonstrate that the model significantly improves energy efficiency while maintaining operational resilience across multiple confidence levels, highlighting its suitability for deployment in intelligent and sustainable urban energy infrastructures.

2. Materials and Methods

2.1. Structure of Electric-Hybrid Hydrogen-Gas Urban Integrated Energy System

In order to reasonably construct the EHUIES low-carbon optimization scheduling model, this section introduces the EHUIES structure in detail, among which the natural gas–hydrogen coupling unit and the co-firing low-carbon multi-energy production module are the core parts of an EHUIES, as shown in Figure 1.

Figure 1. Structure of an EHUIES.

As illustrated in Figure 1, the urban distribution network and the hydrogen-enriched compressed natural gas network (HCNGN) serve as key hubs for hydrogen production and blending within the EHUIES. In the natural gas–hydrogen coupling unit, the water electrolysis system enables the conversion of electrical energy into hydrogen. The produced hydrogen is then blended with natural gas from the gas supply source via the HCNGN. A hydrogen storage tank, located in the hydrogen storage unit, acts as a critical energy buffer, allowing the system to decouple hydrogen production from usage. The hydrogen blending process is regulated by the on/off status of the control switch SH.

The co-combustion low-carbon multi-energy generation module comprises transformers, hydrogen-blended gas turbines (HBGTs), adsorption chillers, and waste heat recovery boilers. Each HBGT consists of an air compressor, combustion chamber, and gas turbine (including its generator section). The HCNG fuel is supplied via a pressure regulating station, which ensures that the gas mixture meets the required compositional specifications. Waste heat from the HBGTs is recovered and utilized by the adsorption chillers and WHRBs, enhancing system energy efficiency.

The energy supply section of EHUIES, connected to the urban distribution network, integrates wind power, photovoltaic (PV) power, and thermal power sources. The energy conversion units include electric boilers and electric chillers, while battery energy storage systems serve as the primary means of electric energy storage. On the demand side, EHUIES supports four types of loads: electricity, heating, cooling, and gas.

2.2. Low-Carbon Optimization Dispatch Model for Electricity–Hydrogen–Natural Gas Hybrid Urban Integrated Energy System

This section introduces the low-carbon optimal scheduling model of EHUIES, as depicted in Figure 1, focusing on the formulation of the objective function and the associated system constraints.

2.2.1. Optimization Goals

The low-carbon optimal scheduling model of the EHUIES aims to minimize the total system operating cost F, which serves as the objective function. The total cost F consists of multiple components, including the operation cost of thermal power units F_pg, the operation cost of water electrolysis units F_we, wind power curtailment cost F_aw, carbon emission cost F_ce, natural gas supply cost from the external source F_ng, HCNG purchase cost of the user unit F_hg, electricity purchase cost of the user unit F_pe, operation and maintenance cost of system equipment F_om, and the cost associated with the hydrogen storage unit F_hs. The variable μ_hs denotes the operating state of the hydrogen storage unit, which is primarily determined by the on/off status of the control switch S_H. This variable plays a key role in subsequent scenario analyses. The objective function is expressed as follows:

min 𝐹 = 𝐹 pg + 𝐹 we + 𝐹 aw + 𝐹 ce + 𝐹 ng + 𝐹 hg + 𝐹 pe + 𝐹 om + 𝜇 hs 𝐹 hs, 𝜇 hs = {0, S H is open 1, S H is closed

(1)

The operation cost of thermal power generation units, denoted as F_pg, can be calculated as follows:

𝐹 pg = \sum 𝑡 = 1 𝑇 \sum 𝑔 = 1 𝐺 (𝑎 𝑔 𝑃 2 𝑔, 𝑡 + 𝑏 𝑔 𝑃 𝑔, 𝑡 + 𝑐 𝑔)

(2)

where T is the scheduling time horizon, G is the total number of thermal power units, a_g, b_g, c_g are the cost coefficients of thermal power unit g, and P_g_,t is the power output of unit g at time t.

The operation cost of the water electrolysis units, denoted as F_we, is expressed as follows:

𝐹 we = 𝜁 we \sum 𝑡 = 1 𝑇 \sum 𝑛 = 1 𝑁 𝑃 we, 𝑛, 𝑡

(3)

where N is the number of water electrolysis units, ζ_we is the operation cost coefficient of the water electrolysis unit, P_we,n,t is the power consumption of electrolysis unit n at time t.

During the wind power integration process, the associated wind curtailment cost, denoted as F_aw, is calculated as follows:

𝐹 aw = 𝛿 aw \sum 𝑡 = 1 𝑇 \sum 𝑤 = 1 𝑊 (𝑃 ˜ 𝑤, 𝑡 - 𝑃 𝑤, 𝑡)

(4)

where N is the number of wind farms, δ_aw is the wind curtailment price, ${\tilde{P}}_{w, t}$ is the forecasted power of wind farm w at time t, P_w_, is the actual utilized power of wind farm w at time t.

The carbon emission cost, denoted as F_ce, is calculated as follows:

𝐹 ce = 𝜉 ce \sum 𝑡 = 1 𝑇 \sum 𝑚 = 1 𝑀 𝑄 cab, 𝑚, 𝑡

(5)

where M is the number of HBGTs; ξ_ce is the carbon tax price.

Since the natural gas source supplies natural gas to the HCNGN, which is then blended with the injected hydrogen in a non-fixed ratio before being sold as HCNG to user units, the natural gas supply cost F_ng and the user-side HCNG purchase cost F_hg are calculated as follows:

𝐹 ng = 𝛾 ng \sum 𝑡 = 1 𝑇 \sum 𝑠 = 1 𝑆 𝑞 𝑠, 𝑡 Δ 𝑡

(6)

𝐹 hg = \sum 𝑡 = 1 𝑇 \sum 𝑚 = 1 𝑀 𝛾 hg, 𝑡 𝑉 hg, 𝑚, 𝑡

(7)

where M is the number of natural gas sources, γ_ng is the natural gas price, γ_hg,t is the price of HCNG at time t, Δt is the scheduling time interval, and q_s,t is the natural gas supply flow rate from source s at time t.

The electricity purchase cost of the user unit, denoted as F_pe, is calculated as follows:

𝐹 pe = \sum 𝑡 = 1 𝑇 \sum 𝑙 = 1 𝐿 𝛾 pe, 𝑡 𝑃 pe, 𝑙, 𝑡

(8)

where L is the number of urban terminal user units, γ_pe,t is the price of electricity sold by the urban distribution network to user units at time t, and P_pe,l,t is the power purchased by the l-th user unit at time t.

The operation and maintenance cost F_om consists of the costs associated with the operation and maintenance of electrical energy storage systems, HBGTs, electric chillers, adsorption chillers, electric boilers, and waste heat recovery boilers. It is expressed as follows:

𝐹 om = \sum 𝑡 = 1 𝑇 ⎡⎣⎢ 𝜑 eh 𝑃 eh, 𝑡 + 𝜑 ht \sum 𝑚 = 1 𝑀 𝑃 ht, 𝑚, 𝑡 + 𝜑 ec 𝑃 ec, 𝑡 + 𝜑 rh 𝐻 rh, in, 𝑡 + 𝜑 ac 𝐻 ac, in, 𝑡 + 𝜑 bat \sum 𝑢 = 1 𝑈 (𝑃 ch, 𝑢, 𝑡 + 𝑃 dis, 𝑢, 𝑡) ⎤⎦⎥

(9)

where U is the number of electrical energy storage units, φ_eh, φ_ht, φ_ec, φ_rh, φ_ac, and φ_bat represent the operation and maintenance cost coefficients of the electric boiler, electric chiller, waste heat recovery boiler, adsorption chiller, HBGT, and electrical energy storage unit, respectively; P_eh,t is the electric power consumed by the electric boiler at time t, P_ec,t is the electric power consumed by the electric chiller at time t, and P_ch,u,t and P_dis,u,t are the charging and discharging power of the u-th electrical energy storage unit at time t, respectively. H_rh,in,t is the thermal input power to the waste heat recovery boiler at time t, and H_ac,in,t is the thermal input power to the adsorption chiller at time t.

The cost of the hydrogen storage unit F_hs primarily consists of the hydrogen storage investment cost and the hydrogen storage operational cost, which can be expressed as follows:

𝐹 hs = 𝑐 invest, hst 𝑁 hst  𝑖 𝑛 𝑣 𝑒 𝑠 𝑡 𝑚 𝑒 𝑛 𝑡 𝑐 𝑜 𝑠 𝑡 + 𝛼 hs \sum 𝑡 = 1 𝑇 \sum ℎ = 1 𝑁 h (𝑞 in hs, ℎ, 𝑡 + 𝑞 out hs, ℎ, 𝑡) Δ 𝑡  𝑜 𝑝 𝑒 𝑟 𝑎 𝑡 𝑖 𝑜 𝑛 𝑎 𝑙 𝑐 𝑜 𝑠 𝑡

(10)

where N_h is the number of hydrogen storage units, c_invest,hst is the price of a hydrogen storage tank, N_hst is the number of hydrogen storage tanks, α_hs is the hydrogen storage cost coefficient, $q_{hs, h, t}^{in}$ and $q_{hs, h, t}^{out}$ are the hydrogen inflow and outflow rates of the h-th hydrogen storage unit at time t, respectively.

2.2.2. Constraints

Operating Constraints of Urban Distribution Systems

(1): Thermal Power Unit Operation Constraints

The operation of thermal power units is subject to power output constraints, which are specifically expressed as follows:

𝑃 𝑔, min \leq 𝑃 𝑔, 𝑡 \leq 𝑃 𝑔, max

(11)

The ramp rate constraints for the operation of thermal power units can be expressed as follows:

{𝑃 𝑔, 𝑡 - 1 - 𝑃 𝑔, 𝑡 \leq 𝑟 ̂ d, 𝑔 Δ 𝑡 𝑃 𝑔, 𝑡 - 𝑃 𝑔, 𝑡 - 1 \leq 𝑟 ̂ u, 𝑔 Δ 𝑡

(12)

where ${\hat{r}}_{d, g}$ and ${\hat{r}}_{u, g}$ are the downward-climbing rate and upward-climbing rate of the g-th thermal power unit, respectively.

(2): Operational Constraints of Wind Farms

In the deterministic day-ahead optimal scheduling of EHUIES, it is assumed that the actual power output of the wind farm does not exceed its forecasted power, which can be expressed as follows:

0 \leq 𝑃 𝑤, 𝑡 \leq 𝑃 ˜ 𝑤, 𝑡

(13)

(3): Power Balance Constraint

Taking the urban distribution network as the core framework for electrical energy, the operation of the water electrolysis units, thermal power units, wind farms, and photovoltaic power stations are comprehensively considered, along with the electricity purchasing behavior of some user units and the electrical load of the urban distribution network. The power balance of the urban distribution network is calculated according to the following equation:

\sum 𝑔 = 1 𝐺 𝑃 𝑔, 𝑡 + \sum 𝑤 = 1 𝑊 𝑃 𝑤, 𝑡 + 𝑃 UD vps, 𝑡 = \sum 𝑙 = 1 𝐿 𝑃 ep, 𝑙, 𝑡 + \sum 𝑛 = 1 𝑁 𝑃 we, 𝑛, 𝑡 + 𝑃 UD load, 𝑡

(14)

where $P_{vps, t}^{UD}$ is the predicted power generation of photovoltaic power stations in urban distribution network at time t, $P_{load, t}^{UD}$ is electric load of urban distribution network at time t, which is given by $P_{load, t}^{UD} = \sum_{d = 1}^{N_{d}} P_{load, d, t}$ , where N_d is the number of load nodes, and P_load,d,t is the electrical load at load node d in the urban distribution network at time t.

(4): Transmission Capacity Constraints of Urban Distribution Network Lines

For the power flow constraints in the urban distribution network, the transmission capacity limits of each power line must be considered to ensure compliance with the actual operating conditions of the urban distribution network. These constraints can be expressed as follows:

- 𝑃 tc 𝑥, max \leq \sum 𝑔 = 1 𝐺 𝓁 𝑥, 𝑔 𝑃 𝑔, 𝑡 + \sum 𝑤 = 1 𝑊 𝓁 𝑥, 𝑤 𝑃 𝑤, 𝑡 + 𝓁 𝑥, vps 𝑃 UD vps, 𝑡 - \sum 𝑑 = 1 𝑁 d 𝓁 𝑥, 𝑑 𝑃 load, 𝑑, 𝑡 - \sum 𝑙 = 1 𝐿 𝓁 𝑥, 𝑙 𝑃 ep, 𝑙, 𝑡 - \sum 𝑛 = 1 𝑁 𝓁 𝑥, 𝑛 𝑃 we, 𝑛, 𝑡 \leq 𝑃 tc 𝑥, max

(15)

where $P_{x, \max}^{tc}$ is the upper limit of transmission capacity of line x, $l_{x, g}$ , $l_{x, w}$ , $l_{x, vps}$ , $l_{x, d}$ , $l_{x, l}$ and $l_{x, n}$ are the transfer distribution factors of the g-th thermal power unit, the w-th wind farm, the photovoltaic power station, the load node d, the l-th user unit, and the n-th water electrolysis device on the line x.

Operating Constraints of the Hydrogen-Enriched Natural Gas Network

Assuming that HCNG flows in the HCNGN under constant temperature and variable pressure, the inlet temperature of the HBGT is derived from the node outlet pressure (after secondary regulation) using the ideal gas law. Considering HCNG component tracking, the HCNGN must satisfy the following constraints:

(1): Pipeline Flow–Pressure Relationship Constraint

Assuming HCNG flows directionally, the pressure–flow relationship in the gas network is formulated as follows:

𝑞 𝑖 𝑗, 𝑡 | 𝑞 𝑖 𝑗, 𝑡 | = π 2 𝑅 64 \cdot (𝑇 st 𝑝 st) 2 \cdot 𝑝 2 𝑖 , 𝑡 - 𝑝 2 𝑗 , 𝑡 𝜌 𝑖 , 𝑡 𝐿 𝑖 𝑗 𝑇 hg 𝛷 𝑖 𝑗 𝑓 𝑖 𝑗 \cdot 𝐷 5 𝑖 𝑗

(16)

where q_ij,t is the pipeline flow between gas network nodes i and j at time t, p_i_,t is the pressure of gas network node i at time t, T_st is HCNG temperature under standard conditions, p_st is HCNG pressure under standard conditions, D_ij, L_ij, and f_ij are the pipe diameter, pipe length, and friction coefficient between gas network nodes i and j, respectively. T_hg is HCNG pipeline temperature, Φ_ij is the compression factor between gas grid nodes i and j.

Within safety limits, the HCNGN must satisfy both pipeline flow and node pressure constraints, as shown below.

- 𝑞 𝑖 𝑗, max \leq 𝑞 𝑖 𝑗, 𝑡 \leq 𝑞 𝑖 𝑗, max

(17)

𝑝 𝑖, min \leq 𝑝 𝑖, 𝑡 \leq 𝑝 𝑖, max

(18)

(2): Energy Balance Constraint at Gas Network Nodes

In the HCNGN, the pipeline flow into node i is defined as positive, as shown in Figure 2.

Figure 2. Positive direction of pipeline flow.

Any gas network node balances energy flow (power) and is related to the gas source points, hydrogen blending points, and other adjacent gas network nodes. Additionally, the reference flow direction in the HCNGN pipelines needs to be defined, where Lq(j,i) represents the set of all pipelines with reference direction from node j to node i, and Lq(i,j) represents the set of all pipelines with reference direction from node i to node j. Therefore, the energy balance at a gas network node in the HCNGN can be calculated using the following equation:

𝜒 H 2, 𝑖 𝑞 H 2, 𝑖, 𝑡 𝐻 ̂ gas, 𝑧 + 𝜒 𝑖, 𝑠 𝑞 𝑠, 𝑡 𝐻 ̂ ng, 𝑠 + \sum 𝑖 𝑗 \in 𝐿 q (𝑗, 𝑖) 𝑢 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝐺 ̂ hg, 𝑗, 𝑡 - \sum 𝑖 𝑗 \in 𝐿 q (𝑖, 𝑗) 𝑢 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝐺 ̂ hg, 𝑖, 𝑡 + \sum 𝑖 𝑗 \in 𝐿 q (𝑖, 𝑗) 𝛿 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝐺 ̂ hg, 𝑗, 𝑡 - \sum 𝑖 𝑗 \in 𝐿 q (𝑗, 𝑖) 𝛿 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝐺 ̂ hg, 𝑖, 𝑡 = 𝐸 d, 𝑖, 𝑡 + 𝐸 ht, 𝑖, 𝑡, 𝑧 \in 𝑧 ̂ H 2

(19)

where $χ_{H_{2}, i}$ is the hydrogen mixing status of gas network node i, which can be 0 or 1, ${\hat{H}}_{gas, z}$ is high calorific value of gas component z, ${\hat{H}}_{ng, s}$ is the high calorific value of natural gas from the s-th gas source, ${\hat{G}}_{hg, i, t}$ is the HCNG calorific value of gas grid node i at time t, E_d,i,t and E_ht,i,t are the conventional gas load and HBGT load of gas grid node i at time t, respectively. u_ij_,t and δ_ij_,t are the directional indicator variables for the actual pipeline flow between gas network nodes i and j at time t, representing the positive and negative flow directions, respectively. Their calculation formulas are as follows:

⎧⎩⎨ 𝜀 𝑖 𝑗, 𝑡 = sgn (𝑞 𝑖 𝑗, 𝑡) 𝑢 𝑖 𝑗, 𝑡 = | 𝜀 𝑖 𝑗 , 𝑡 | + 𝜀 𝑖 𝑗 , 𝑡 2 𝛿 𝑖 𝑗, 𝑡 = | 𝜀 𝑖 𝑗 , 𝑡 | - 𝜀 𝑖 𝑗 , 𝑡 2

(20)

where sgn(·) is the sign function, ε_ij,t is the actual pipeline flow direction between gas network nodes i and j at time t.

The supply flow constraint of natural gas sources is expressed as follows:

𝑞 𝑠, min \leq 𝑞 𝑠, 𝑡 \leq 𝑞 𝑠, max

(21)

(3): Pipeline Average Pressure Constraint

Based on the Soave–Redlich–Kwong equation of state, its form is further improved to describe the average pipeline pressure, which is expressed as follows:

⎧⎩⎨ 𝑝 av, 𝑖 𝑗, 𝑡 = 2 3 \cdot (𝑝 𝑖, 𝑡 + 𝑝 𝑗, 𝑡 - 𝑝 𝑖 , 𝑡 𝑝 𝑗 , 𝑡 𝑝 𝑖 , 𝑡 + 𝑝 𝑗 , 𝑡) 𝛷 3 𝑖 𝑗 - 𝛷 2 𝑖 𝑗 + 𝛷 𝑖 𝑗 (𝐴 𝑖 𝑗, 𝑡 - 𝐵 𝑖 𝑗, 𝑡 - 𝐵 2 𝑖 𝑗, 𝑡) - 𝐴 𝑖 𝑗, 𝑡 𝐵 𝑖 𝑗, 𝑡 = 0 𝐴 𝑖 𝑗, 𝑡 = 0.42747 \cdot 𝑝 av , 𝑖 𝑗 , 𝑡 𝑇 2 hg \cdot ⎛⎝⎜⎜ \sum 𝑧 \in 𝑍 𝑁 𝑧 , 𝑖 , 𝑡 𝑇 c , 𝑧 𝜙 𝑧 𝑝 c , 𝑧 ---\sqrt ⎞⎠⎟⎟ 2 𝐵 𝑖 𝑗, 𝑡 = 0.08664 \cdot 𝑝 av , 𝑖 𝑗 , 𝑡 𝑇 hg \cdot \sum 𝑧 \in 𝑍 𝑁 𝑧 , 𝑖 , 𝑡 𝑇 c , 𝑧 𝑝 c , 𝑧 𝜙 𝑧 = 1 + (0.48 + 1.574 𝜔 𝑧 - 0.176 𝜔 2 𝑧) \cdot ⎛⎝⎜⎜ 1 - 𝑇 hg 𝑇 c , 𝑧 ----\sqrt ⎞⎠⎟⎟

(22)

where p_av,ij,t is average pipeline pressure between gas grid nodes i and j at time t, T_c,z and p_c,z are critical temperature and critical pressure of gas component z, Z_ij is the pipeline compression coefficient between gas network nodes i and j, A_ij_,t and B_ij_,t are the HCNG state parameters between gas grid nodes i and j at time t, ϕ_z is the eccentricity coefficient conversion parameter of gas component z.

(4): Calculation of Gas Component Molar Fractions

During the HCNGN balancing process, the HCNG composition varies at each gas network node. Based on the node flow balance, the molar fraction of gas component z at node i and time t, N_z_,i,t can be calculated using the following equation:

𝑁 𝑧, 𝑖, 𝑡 = ⎛⎝⎜⎜⎜⎜ 𝜒 H 2, 𝑖 𝑁 H 2, 𝑧, 𝑖, 𝑡 𝑞 H 2, 𝑖, 𝑡 + 𝜒 𝑖, 𝑠 𝑁 ng, 𝑧, 𝑖, 𝑡 𝑞 𝑠, 𝑡 + \sum 𝑖 𝑗 \in 𝐿 q (𝑗, 𝑖) 𝑢 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝑁 𝑧, 𝑗, 𝑡 - \sum 𝑖 𝑗 \in 𝐿 q (𝑖, 𝑗) 𝑢 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝑁 𝑧, 𝑖, 𝑡 + \sum 𝑖 𝑗 \in 𝐿 q (𝑖, 𝑗) 𝛿 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝑁 𝑧, 𝑗, 𝑡 - \sum 𝑖 𝑗 \in 𝐿 q (𝑗, 𝑖) 𝛿 𝑖 𝑗, 𝑡 𝑞 𝑖 𝑗, 𝑡 𝑁 𝑧, 𝑖, 𝑡 / ⎛⎝⎜⎜⎜⎜ 𝜒 H 2, 𝑖 𝑞 H 2, 𝑖, 𝑡 + 𝜒 𝑖, 𝑠 𝑞 𝑠, 𝑡 + \sum 𝑖 𝑗 \in 𝐿 q (𝑗, 𝑖) 𝑞 𝑖 𝑗, 𝑡 - \sum 𝑖 𝑗 \in 𝐿 q (𝑖, 𝑗) 𝑞 𝑖 𝑗, 𝑡 ⎞⎠⎟⎟⎟⎟

(23)

where $N_{H_{2}, z, i, t}$ is the molar fraction of gas component z in the hydrogen added to the gas grid node i at time t, $N_{ng, z, i, t}$ is the molar fraction of gas component z in the natural gas supplied by the gas source of gas grid node i at time t.

(5): Compressor and Pressure Regulating Station Operation Constraints

Compressors and pressure regulating stations assist HCNGN in maintaining safe and stable HCNG flow by performing secondary pressure regulation at gas network node outlets, ensuring suitable initial combustion conditions for HBGT. Their operational constraints are given by the following equations:

⎧⎩⎨ 𝑝 𝑗, 𝑡 = 𝑟 𝑖 𝑗, 𝑡 𝑝 𝑖, 𝑡 𝑇 in, 𝑖, 𝑡 = 𝜛 tf 𝑝 𝑖 , 𝑡 𝑅 𝑘 𝑖 , 𝑡

(24)

where r_ij_,t is the compression ratio between gas network nodes i and j at time t, ϖ_tf is the pressure conversion factor, and k_i_,t is the pressure regulation ratio of the outlet side of the gas network node i at time t.

Operational Constraints of Urban End-User Units

To accurately represent EHUIES operation, equipment classification and power balance are detailed for different urban end-user types based on Figure 1. User units are categorized into residential communities, university campuses, and industrial parks. Residential and university areas mainly include HBGT, photovoltaic systems, and electrical energy storage, focusing on electrical power balance. Industrial parks incorporate HBGT, PV units, electric boilers, chillers, waste heat recovery boilers, adsorption chillers, and energy storage, requiring electrical, thermal, and cooling power balance.

(1): Charging and Discharging Constraints of Electrical Energy Storage

𝐸 st, 𝑢, 𝑡 = 𝐸 st, 𝑢, 𝑡 - 1 + (𝜂 st 𝑃 ch, 𝑢, 𝑡 - 𝑃 dis , 𝑢 , 𝑡 𝜂 st) Δ 𝑡

(25)

0 \leq 𝐸 st, 𝑢, 𝑡 \leq 𝐸 st, 𝑢, max

(26)

where E_st,u,t is the capacity of the u-th energy storage at time t, and η_st is the conversion efficiency of the energy storage.

The state transition between charging and discharging of electrical energy storage systems is subject to a limited number of cycles. Moreover, charging and discharging cannot occur simultaneously, and the respective charging and discharging power must not exceed their allowable maximum values. The constraints can be expressed as follows:

⎧⎩⎨ \sum 𝑡 = 1 𝑇 | 𝑈 ch, 𝑢, 𝑡 - 𝑈 dis, 𝑢, 𝑡 | \leq 𝑁 bat, 𝑢 𝑈 dis, 𝑢, 𝑡 + 𝑈 ch, 𝑢, 𝑡 \leq 1 0 \leq 𝑃 ch, 𝑢, 𝑡 \leq 𝑈 ch, 𝑢, 𝑡 𝑃 ch, 𝑢, max 0 \leq 𝑃 dis, 𝑢, 𝑡 \leq 𝑈 dis, 𝑢, 𝑡 𝑃 dis, 𝑢, max

(27)

where U_ch,u,t and U_dis,u,t are the charging and discharging states of the u-th energy storage at time t, respectively, P_ch,u,max and P_dis,u,max are the maximum charging and discharging powers allowed by the u-th energy storage, respectively, and N_bat,u is the number of charging and discharging conversions of the u-th energy storage.

(2): Operational Constraints of Waste Heat Recovery Boilers and Adsorption Chillers

The input heat for both the waste heat recovery boiler and the adsorption chiller is supplied by the HBGT combustion chamber. Their collected heat power must not exceed their respective maximum allowable limits. The operational constraints are given as follows:

⎧⎩⎨ 𝐻 rh, 𝑡 = 𝜂 rb, h 𝐻 rh, in, 𝑡 𝛽 rb, h 𝐶 ac, 𝑡 = 𝜂 ar, c 𝐻 ac, in, 𝑡 𝛽 ar, c 0 \leq 𝐻 rh, in, 𝑡 \leq 𝐻 rh, in, max 0 \leq 𝐻 ac, in, 𝑡 \leq 𝐻 ac, in, max

(28)

where η_rb,h and η_ar,c are the heat collection efficiency of the waste heat recovery boiler and the cooling efficiency of the adsorption chiller, respectively, β_rb,h and β_ar,c are the heating coefficient and the cooling coefficient, respectively, H_rh,t is the heating power on the output side of the waste heat recovery boiler at time t, C_ac,t is the cooling power on the output side of the adsorption chiller at time t, H_rh,in,max is the maximum heat collection power allowed on the input side of the waste heat recovery boiler, and H_ac,in,max is the maximum heat collection power allowed on the input side of the adsorption chiller.

In industrial parks, the thermal power output of HBGT and the corresponding heat transfer balance can be expressed as follows:

𝐻 ht, 𝑚, 𝑡 = 𝜂 re 𝑃 ht, 𝑚, 𝑡, 𝑚 \in 𝐿 ̂ ht, 3

(29)

\sum 𝑚 \in 𝐿 ̂ ht, 3 𝐻 ht, 𝑚, 𝑡 = 𝐻 rh, in, 𝑡 + 𝐻 ac, in, 𝑡

(30)

where η_re is the HBGT heat conversion efficiency, and $H_{ht, m, t}$ is the heating power of the m-th HBGT at time t.

(3): Operational Constraints of Electric Boilers and Electric Chillers

Electric boilers and electric chillers use electricity as their input energy source, converting it into heat and cooling, respectively. Their operational constraints ensure that the power consumption does not exceed the allowable maximum values and are defined as follows.

⎧⎩⎨ 𝐻 eh, 𝑡 = 𝜂 eh 𝑃 eh, 𝑡 𝐶 ec, 𝑡 = 𝜂 ec 𝑃 ec, 𝑡 0 \leq 𝑃 eh, 𝑡 \leq 𝑃 eh, max 0 \leq 𝑃 ec, 𝑡 \leq 𝑃 ec, max

(31)

where η_eh is the heating efficiency of the electric boiler, η_ec is the cooling efficiency of the electric refrigerator, C_ec,t is the cooling power on the output side of the electric refrigerator at time t, and P_eh,max and P_ec,max are the maximum power consumption allowed by the electric boiler and the electric refrigerator, respectively.

(4): Multi-energy Balance Constraint

The power balance equations for different types of user units are expressed as follows:

\sum 𝑚 \in 𝐿 ̂ ht, 1 \cup 𝑚 \in 𝐿 ̂ ht, 2 𝑃 ht, 𝑚, 𝑡 + 𝑃 pe, 𝑙, 𝑡 + 𝑃 US v, 𝑙, 𝑡 + \sum 𝑢 \in 𝐿 ̂ es, 1 \cup 𝑢 \in 𝐿 ̂ es, 2 𝑃 dis, 𝑢, 𝑡 = 𝑃 US load, 𝑙, 𝑡 + \sum 𝑢 \in 𝐿 ̂ es, 1 \cup 𝑢 \in 𝐿 ̂ es, 2 𝑃 ch, 𝑢, 𝑡, 𝑙 \in {𝐿 ̂ us, 1, 𝐿 ̂ us, 2}

(32)

\sum 𝑚 \in 𝐿 ̂ ht, 3 𝑃 ht, 𝑚, 𝑡 + 𝑃 pe, 𝑙, 𝑡 + 𝑃 US v, 𝑙, 𝑡 + \sum 𝑢 \in 𝐿 ̂ es, 3 𝑃 dis, 𝑢, 𝑡 = 𝑃 US load, 𝑙, 𝑡 + \sum 𝑢 \in 𝐿 ̂ es, 3 𝑃 ch, 𝑢, 𝑡 + 𝑃 eh, 𝑡 + 𝑃 ec, 𝑡, 𝑙 \in 𝐿 ̂ us, 3

(33)

𝐻 rh, 𝑡 + 𝐻 eh, 𝑡 = 𝐻 load, 𝑡

(34)

𝐶 ec, 𝑡 + 𝐶 ac, 𝑡 = 𝐶 load, 𝑡

(35)

where $P_{v, l, t}^{US}$ is the predicted photovoltaic power of the l-th user unit at time t, $P_{load, l, t}^{US}$ is the electrical load of the l-th user unit at time t, H_load,t and C_load,t are the heat and cold loads of the industrial park at time t, respectively.

Coupled Operational Constraints

(1): HBGT operational constraints

The operational constraints and boundary conditions of the HBGT are as follows:

𝜆 wh, 𝑧, 𝑚, 𝑡 = {𝜆 𝑧, 𝑚, 𝑡, 0 < 𝑇 in, 𝑚, 𝑡 < 𝑇 com, 𝑧 𝜆 com, 𝑧, 𝑇 in, 𝑚, 𝑡 \geq 𝑇 com, 𝑧

(36)

𝑄 cab, 𝑚, 𝑡 = \sum 𝑧 \in 𝑍 ̂ AG 𝜆 wh, 𝑧, 𝑚, 𝑡 𝛼 nv, 𝑧 𝑁 in, 𝑧, 𝑚, 𝑡 𝑉 hg, 𝑚, 𝑡

(37)

𝑃 ht, 𝑚, 𝑡 = 𝜂 c \sum 𝑧 \in 𝑍 ̂ [𝜓 𝑧 𝑁 2 in, 𝑧, 𝑚, 𝑡 𝜌 hg , 𝑚 , 𝑡 𝑀 hg , 𝑚 , 𝑡 \cdot 𝑖 𝑧 2 𝑅 (𝑇 com, 𝑧 - 𝑇 in, 𝑚, 𝑡) + 𝑖 𝑧 2 𝑅 𝑇 in, 𝑚, 𝑡 𝜆 wh, 𝑧, 𝑚, 𝑡 𝜓 𝑧 𝑁 in, 𝑧, 𝑚, 𝑡] \cdot 𝑉 hg, 𝑚, 𝑡

(38)

⎧⎩⎨ 𝑁 𝑧, 𝑖, 𝑡 = 𝑁 in, 𝑧, 𝑚, 𝑡 𝜌 𝑖, 𝑡 = 𝜌 hg, 𝑚, 𝑡 𝑇 in, 𝑖, 𝑡 = 𝑇 in, 𝑚, 𝑡 𝑉 hg, 𝑖, 𝑡 = 𝑉 hg, 𝑚, 𝑡

(39)

0 \leq 𝑉 hg, 𝑚, 𝑡 \leq 𝑉 hg, 𝑚, max

(40)

where m is the number of HBGTs, t is the scheduling time, λ_z is the actual combustion conversion coefficient of gas component z before reaching steady-state reaction equilibrium, λ_com,z is the actual combustion conversion coefficient of gas component z after reaching steady-state reaction equilibrium, λ_wh,z is the actual combustion conversion coefficient of gas component z, Q_cab is the carbon emission, ${\hat{Z}}_{AG}$ is the set of alkane gas types containing carbon elements in HCNG, α_nv,z is the volume conversion coefficient, N_in,z is the initial molar fraction of gas component z, V_hg is the volume of HCNG consumed by HBGT, ψ_z is the molar fraction conversion coefficient, η_c is the energy conversion efficiency between mechanical energy and electrical energy, P_ht is the mixed combustion power generation of HBGT, ρ_hg is the initial density of HCNG in the combustion chamber, M_hg is the initial molar mass of HCNG in the combustion chamber, N_z_,i,t is the mole fraction of gas component z in gas network node i at time t, ρ_i_,t is the HCNG density of gas network node i at time t, T_in,i,t is the temperature of HCNG after pressure regulation in gas network node i at time t, and V_hg,i,t is the volume of HCNG supplied by gas network node i at time t.

(2): Operational Constraints of Water Electrolysis Units

The hydrogen output flow equation of the water electrolysis unit, its power consumption constraints, and the boundary condition constraints between the water electrolysis unit and the gas network nodes are given as follows:

𝑞 H 2, 𝑛, 𝑡 = 𝜂 we 𝑃 we , 𝑛 , 𝑡 𝜎 hc 𝐻 ̂ gas , 𝑧, 𝑧 \in 𝑧 ̂ H 2

(41)

0 \leq 𝑃 we, 𝑛, 𝑡 \leq 𝑃 we, 𝑛, max

(42)

𝑞 H 2, 𝑛, 𝑡 = 𝑞 H 2, 𝑖, 𝑡

(43)

where η_we is the efficiency of hydrogen production using water electrolysis, and σ_hc is the correction coefficient of the high calorific value of hydrogen.

(3): Hydrogen storage constraints

The real-time balance relationship of the hydrogen storage tank, the maximum hydrogen storage capacity constraint, and the constraints on the hydrogen inflow and outflow rates to/from the storage unit are expressed as follows:

𝐸 hs, ℎ, 𝑡 = 𝐸 hs, ℎ, 𝑡 - 1 + ⎛⎝⎜⎜ 𝜂 hs 𝑞 in hs, ℎ, 𝑡 - 𝑞 out hs , ℎ , 𝑡 𝜂 H 2 ⎞⎠⎟⎟ Δ 𝑡

(44)

0 \leq 𝐸 hs, ℎ, 𝑡 \leq 𝐸 hs, ℎ, max

(45)

⎧⎩⎨ 0 \leq 𝑞 in hs, ℎ, 𝑡 \leq 𝑈 in hs, ℎ, 𝑡 𝑞 in hs, ℎ, max 0 \leq 𝑞 out hs, ℎ, 𝑡 \leq 𝑈 out hs, ℎ, 𝑡 𝑞 out hs, ℎ, max 𝑈 in hs, ℎ, 𝑡 + 𝑈 out hs, ℎ, 𝑡 \leq 1

(46)

where $η_{hs}$ is the conversion efficiency of the hydrogen storage tank.

2.3. Multiple Uncertainty Analysis

2.3.1. Source/Load Side Uncertainty

Considering that there exist varying deviations between the actual and forecasted values of wind power generation and load demand at different time periods, the fluctuation range is described using the ratio between actual and forecasted values. The load demand is mainly divided into urban distribution network load and user-side load, which are formulated as follows:

⎧⎩⎨ 𝛺 𝑤, 𝑡 = 𝑃 𝑤 , 𝑡 𝑃 ˜ 𝑤 , 𝑡 \in [1 - 𝛿 ̂ 𝑤, 𝑡, 1 + 𝛿 ̂ 𝑤, 𝑡], 0 \leq 𝛿 ̂ 𝑤, 𝑡 \leq 1 𝛺 UD load, 𝑡 = 𝑃 UD load , 𝑡 𝑃 ˜ UD load , 𝑡 \in [1 - 𝛿 ̂ UD load, 𝑡, 1 + 𝛿 ̂ UD load, 𝑡], 0 \leq 𝛿 ̂ UD load, 𝑡 \leq 1 𝛺 US load, 𝑙, 𝑡 = 𝑃 US load , 𝑙 , 𝑡 𝑃 ˜ US load , 𝑙 , 𝑡 \in [1 - 𝛿 ̂ US load, 𝑙, 𝑡, 1 + 𝛿 ̂ US load, 𝑙, 𝑡], 0 \leq 𝛿 ̂ US load, 𝑙, 𝑡 \leq 1

(47)

where ${\hat{δ}}_{w, t}$ is the power deviation of the w-th wind farm at time t, ${\hat{δ}}_{load, t}^{UD}$ is the load demand deviation of the urban distribution network at time t, ${\hat{δ}}_{load, l, t}^{US}$ is the load demand deviation of the l-th user unit at time t, ${\tilde{P}}_{load, t}^{UD}$ is the predicted value of the urban distribution network load at time t, ${\tilde{P}}_{load, l, t}^{US}$ is the predicted value of the load of the l-th user unit at time t.

The uncertainties of wind power and load demand exhibit distinct stochastic characteristics, which can be described using either symmetric or asymmetric probability distributions. For example, wind power is directly influenced by wind speed, whose stochastic nature approximately follows a Weibull distribution. In contrast, various types of load demand are typically modeled using a normal distribution [27]. It is worth noting that the parameters involved in the probability density functions can be determined by fitting the historical data of wind speed and load (or load forecast errors).

2.3.2. Combustion Reaction Uncertainty

Considering that the HCNG studied in this paper includes multi-component tracking—mainly methane, ethane, propane, and hydrogen—these components exhibit different combustion characteristics due to variations in heating value and physical properties. However, it is difficult to systematically quantify these combustion characteristics in existing research, which mostly relies on numerical simulations. Moreover, influenced by the internal environmental conditions of the HBGT combustion chamber, the outcome of each combustion reaction within the chamber is subject to variability and randomness [28]. Given the complexity and the number of influencing factors involved in HCNG combustion within HBGT, it is challenging to construct an explicit expression for the combustion reaction.

Therefore, to systematically quantify the uncertainty in combustion reactions, this study applies stochastic optimization theory to approximate the “combustion reaction region” of different gas components in HCNG based on available reference data, thus characterizing the uncertain operational characteristics of HBGT. Assuming that the combustion uncertainty is independent of the HBGT unit index m, and considering that the intake temperature of the HBGT determines the combustion process, this section reflects the randomness and fluctuation of the combustion outcome through variations in the actual steady-state combustion conversion coefficient of gas component z. Instead of assuming it as a constant, a combustion reaction margin $ϑ_{sm, z, t}$ is introduced to describe the non-constant steady-state combustion conversion coefficient. Accordingly, the improved calculation formula for the actual combustion conversion coefficient of component z is expressed as follows:

𝜆 ̂ wh, 𝑧, 𝑚, 𝑡 = ⎧⎩⎨ 𝜗 sm, 𝑧, 𝑡 \cdot 𝜆 𝑧 , 𝑚 , 𝑡 𝜆 com , 𝑧, 0 < 𝑇 in, 𝑚, 𝑡 < 𝑇 com, 𝑧 𝜗 sm, 𝑧, 𝑡, 𝑇 in, 𝑚, 𝑡 \geq 𝑇 com, 𝑧

(48)

Additionally, based on existing reference data for combustion reactions, experimental reference values for the steady-state combustion reaction margin of gas component zzz are obtained. However, at different time points, the actual steady-state combustion reaction margin of component z may deviate from these experimental reference values. This deviation is described using a fluctuation range, as shown in Equation (49).

(1 - 𝛿 ̂ co, 𝑧, 𝑡) 𝜗 ˜ sm, 𝑧, 𝑡 \leq 𝜗 sm, 𝑧, 𝑡 \leq (1 + 𝛿 ̂ co, 𝑧, 𝑡) 𝜗 ˜ sm, 𝑧, 𝑡

(49)

where ${\hat{δ}}_{co, z, t}$ is the deviation of the combustion reaction margin of gas component z at time t.

The stochastic characteristics of the combustion reaction are assumed to approximately follow a normal distribution. Given the limited number of actual reference data samples within the same time dimension, it is difficult to directly determine the parameters of the probability density function. Therefore, within a reasonable error margin, an approximate normal distribution sampling method based on experimental reference data is employed to expand the dataset. New data sample sets are generated via random sampling, and the mean

{\hat{μ}}_{co, z, t}

and variance

{\hat{σ}}_{co, z, t}^{2}

of the corresponding normal distribution

Γ_{co, z, t}

at each time step are determined. The mean

{\hat{μ}}_{co, z, t}

is regarded as the experimental reference value for the steady-state combustion reaction margin, the specific implementation process is shown in Figure 3.

Figure 3. Implementation process of combustion reaction uncertainty description.

2.3.3. Categorical Confidence Interval Description

Based on the stochastic characteristics of the aforementioned uncertainties, categorical probabilistic confidence intervals can be employed to correct the fluctuation deviations and determine the deviation bounds under different confidence levels. Given a confidence level τ, there exist an uncertainty variable matrix

ψ

and its corresponding predicted (reference) value matrix

\tilde{ψ}

, with the reference value

\tilde{ψ}

taken as the mean. The probability density function is then used to define the confidence interval boundaries according to different values of τ_c. As illustrated in Figure 4, the categorical confidence interval of the uncertainty variable is bounded by γ_up (upper bound) and γ_down (lower bound).

Figure 4. Classification probability confidence intervals for uncertain variables.

As shown in Figure 4, under different confidence levels τ_c, the uncertainty variables may follow symmetric or asymmetric probability distributions due to their distinct stochastic characteristics. Therefore, it is necessary to determine the minimum distance between γ_up, γ_down, and $\tilde{ψ}$ in order to compute the deviation boundaries under multiple uncertainties. The specific calculation formulas are given as follows:

⎧⎩⎨ 𝛿 ̂ 𝑤, 𝑡 = min [ 𝑃 ˜ 𝑤 , 𝑡 - 𝛾 down ( 𝑃 ˜ 𝑤 , 𝑡 , 𝜏 ) , 𝛾 up ( 𝑃 ˜ 𝑤 , 𝑡 , 𝜏 ) - 𝑃 ˜ 𝑤 , 𝑡 ] 𝑃 ˜ 𝑤 , 𝑡 𝛿 ̂ UD load, 𝑡 = min [ 𝑃 ˜ UD load , 𝑡 - 𝛾 down ( 𝑃 ˜ UD load , 𝑡 , 𝜏 ) , 𝛾 up ( 𝑃 ˜ UD load , 𝑡 , 𝜏 ) - 𝑃 ˜ UD load , 𝑡 ] 𝑃 ˜ UD load , 𝑡 𝛿 ̂ US load, 𝑙, 𝑡 = min [ 𝑃 ˜ US load , 𝑙 , 𝑡 - 𝛾 down ( 𝑃 ˜ US load , 𝑙 , 𝑡 , 𝜏 ) , 𝛾 up ( 𝑃 ˜ US load , 𝑙 , 𝑡 , 𝜏 ) - 𝑃 ˜ US load , 𝑙 , 𝑡 ] 𝑃 ˜ US load , 𝑙 , 𝑡 𝛿 ̂ co, 𝑧, 𝑡 = min [ 𝜗 ˜ sm , 𝑧 , 𝑡 - 𝛾 down ( 𝜗 ˜ sm , 𝑧 , 𝑡 , 𝜏 ) , 𝛾 up ( 𝜗 ˜ sm , 𝑧 , 𝑡 , 𝜏 ) - 𝜗 ˜ sm , 𝑧 , 𝑡 ] 𝜗 ˜ sm , 𝑧 , 𝑡

(50)

3. SOI-IGDT Considering Exergy Boost

3.1. Overview of Exergy Efficiency

Energy is a measure of the state transformation of a system, while exergy, as the usable portion of energy, is a thermodynamic parameter proposed by Rant in 1956. Exergy is defined as the maximum useful work that can be obtained when a system undergoes a reversible thermodynamic process from its initial state to a state in equilibrium with the environment, and it can be regarded as a form of available energy [29]. The stepwise degradation of energy quality in a system can be interpreted as a gradual decline in its ability to perform work.

Since only reversible processes can theoretically achieve complete energy conversion, exergy represents the maximum useful work output or minimum useful work input that can theoretically be achieved. Conversely, the portion of energy that cannot be converted into exergy is referred to as anergy. Any form of energy ${\hat{E}}_{s}$ can be decomposed into exergy ${\hat{E}}_{x}$ and anergy ${\hat{A}}_{n}$ , and is mathematically expressed as follows:

𝐸 ̂ s = 𝐸 ̂ x + 𝐴 ̂ n

(51)

The exergy of a system can be specifically categorized into input exergy ${\hat{E}}_{x, in}$ and output exergy ${\hat{E}}_{x, out}$ . The relationship between them can be expressed by the following equation:

𝐸 ̂ x, in = 𝐸 ̂ x, out + 𝐸 ̂ x, loss

(52)

where ${\hat{E}}_{x, loss}$ is the exergy loss in the system.

From the perspective of different energy supply sources, exergy efficiency is adopted as an evaluation metric to quantify the utilization level of useful energy in the overall system or its subsystems. It helps identify weak links in effective energy use and supports the formulation of rational urban integrated energy system planning, operational optimization, and scheduling strategies—ultimately promoting the system’s development toward “high-exergy” performance. The specific expression is defined as follows:

𝜂 ̂ ex = 𝐸 ̂ x , out 𝐸 ̂ x , in

(53)

For the optimization scheduling of EHUIES considering multiple uncertainties, the system output exergy reflects the energy utilization on the user side, while the system input exergy represents the energy supply on the supply side. A higher exergy efficiency indicates a lower proportion of exergy loss relative to the input exergy during the supply process, suggesting that the system is closer to achieving high-exergy performance.

3.2. Overview of SOI-IGDT

The traditional IGDT does not account for the probability distributions of uncertain variables, and instead only considers the maximum fluctuation range within the [0, 1] interval. This leads to a coarse and oversimplified characterization of uncertainty. Moreover, it applies the same robustness coefficient ζ across all time periods as the boundary condition for optimization, which results in a “temporal symmetry of deviation.” This assumption is rather idealized and may deviate significantly from the actual fluctuation intervals. As demonstrated in Section 2.3.3 on the classification of confidence intervals for multiple uncertainties, the probabilistic distribution characteristics of uncertainties at different time periods can constrain the effective fluctuation range under traditional IGDT. Based on this, the conventional IGDT exhibits the following limitations:

(1): The deviation factor in traditional IGDT requires setting based on the decision-maker’s experience, which is highly subjective and lacks risk assessment of the model optimization results.
(2): When considering the stochastic characteristics of multiple uncertainties, the maximum fluctuation range obtained using traditional IGDT optimization tends to be overly conservative and may not accurately reflect the actual fluctuation range.
(3): The traditional IGDT’s maximum boundary condition is linearly consistent (symmetric) across time, using a global deviation parameter ζ as the optimization target. This leads to an overly conservative robust solution that lacks effective integration with actual uncertainty fluctuations. Moreover, the impact of maximum deviations on total system cost varies over time, causing inconsistency in the objective function’s direction. When actual uncertainty bounds at different times are known, using the original ζ as the optimization target may result in infeasibility, increased nonlinearity, and reduced solution reliability and efficiency.

Therefore, to address the above shortcomings when considering the stochastic characteristics of multiple uncertainties, this section proposes a SOI-IGDT by integrating IGDT with stochastic optimization theory. The specific improvements are as follows:

(1): First, by applying chance constraints from stochastic optimization theory, the traditional IGDT deviation factor is removed. The constraint satisfaction probability Pr is guaranteed to be no less than the given confidence level τ. Thus, the original IGDT constraint on the optimization objective can be reformulated as follows:

max 𝐹 (𝜿 ̂, 𝝍) \leq 𝐹 ̂ et = (1 + 𝜙 d) 𝐹 ac \to Pr {max 𝐹 (𝜿 ̂, 𝝍) \leq 𝐹 ac} \geq 𝜏

(54)

Accordingly, considering the stochastic characteristics of multiple uncertainties, the actual boundary conditions of the uncertainty variables $ψ$ need to be classified and adjusted based on the confidence level τ. Combining the content from Section 2.3.3 and using the boundary-solving method in Equation (50), the actual boundary conditions of multiple uncertainties are corrected as follows:

𝝍 \in 𝛷 ̂ (𝝍 ˜, 𝜏) = {𝝍 | 𝝍 - (𝝍 ˜, 𝝉) \leq 𝝍 \leq 𝝍 + (𝝍 ˜, 𝝉)}

(55)

where $\hat{Φ}$ is the corrected confidence interval of the uncertainty variable, $ψ^{-}$ is the lower bound of the corrected confidence interval, and $ψ^{+}$ is the upper bound of the corrected confidence interval.

(2): Secondly, the original deviation ζ is a global optimization variable with consistent boundary conditions across all time periods, only satisfying a unified maximum deviation. This approach is unsuitable for scenarios where uncertainty varies over time. Since deviations at different times impact the system differently, the optimization must consider the time-varying effects of deviations on total system cost. Therefore, the model accounts for the impact of actual deviations at each time step during optimization. Furthermore, under multiple uncertainties, it is necessary to decompose multiple objective deviations by time periods and establish an uncertainty deviation matrix ${\hat{ζ}}_{r}$ with time scale t. Using the time-varying deviation $ζ_{r, t}$ allows better tracking of deviation boundary changes considering the uncertainty probability distributions at different times. The improved uncertainty variable fluctuation range is expressed as follows:

𝜻 ̂ 𝑟 = [𝜁 𝑟, 1, 𝜁 𝑟, 2, \dots, 𝜁 𝑟, 𝑡], \forall 𝑟 \in [1, 𝑅 uv], 𝑡 \in [1, 𝑇 vd, 𝑟]

(56)

𝝍 \in 𝛷 ̂ (𝝍 ˜, 𝜏) = {𝝍 | 𝝍 - (𝝍 ˜, 𝜏) \leq 𝝍 \leq 𝝍 + (𝝍 ˜, 𝜏)}

(57)

where $\hat{ζ}$ is the time-varying deviation matrix for multiple uncertainties, R_uv is the number of uncertainty variables, T_vd,r is the effective time period of the forecast (reference) value for the r-th uncertainty variable, determined by whether its value is zero at time t.

Uncertainty probability distributions affect the variation of deviation boundaries at different time points, making it difficult to accurately determine both boundaries of the fluctuation interval. Therefore, a time-varying deviation direction indicator is introduced into the constraints to reflect the actual fluctuation direction of BBB around AAA at each time, eliminating the subjective assumption of a single deviation direction in traditional IGDT. In practice, while managing uncertainty risks and satisfying Equation (54), the maximum fluctuation intervals obtained at different times exhibit mutual coordination, avoiding the one-sided degradation observed in conventional IGDT. To capture the chance-constrained risk-averse feature, based on the expected cost under chance constraints, a relaxed directional condition for uncertainty fluctuation intervals at each time is established. This formulates a deviation maximization problem that ensures the expected cost does not exceed a given threshold at a specified confidence level, thereby removing the max operator in Equation (54) and expanding the solution space of the model. The formulation is expressed as follows:

⎧⎩⎨ (𝒖 ̂ + - 𝒖 ̂ -) 𝝍 + (𝒖 ̂ + - 𝒖 ̂ -) 𝝍 ˜ = 𝜻 ̂ ⊙ 𝝍 ˜ 𝝌 ̂ = sign (𝝍 - 𝝍 ˜), 𝒖 ̂ + = | 𝝌 ̂ | + 𝝌 ̂ 2, 𝒖 ̂ - = | 𝝌 ̂ | - 𝝌 ̂ 2

(58)

where $\hat{χ}$ , ${\hat{u}}_{+}$ , and ${\hat{u}}_{-}$ denote the directional indicator matrix of uncertainty fluctuations, the positive direction indicator, and the negative direction indicator, respectively. $⊙$ represents the Hadamard product operator (element-wise matrix multiplication).

(3): In traditional IGDT, the optimization of multiple uncertainties relies on subjectively assigned priority levels, aiming to maximize deviation under idealized conditions without considering the actual probability distributions. However, in practical system optimization, multiple uncertainties are often correlated and constrained by system balance and equipment coupling. When the true probabilistic characteristics are considered, the feasible region is further restricted by the adjusted confidence intervals. As these intervals vary significantly over time, the predefined global priority scheme may hinder the optimization process. Therefore, without incorporating the actual distributions of uncertainties, the traditional IGDT approach tends to be overly conservative. Therefore, based on the specific improvements introduced in Step (2), $\forall r \in [1, R_{uv}]$ , vor the time-varying deviation matrix ${\hat{ζ}}_{r}$ of the r-th uncertainty variable, the time-varying deviations $ζ_{r, t}$ at different time steps are simultaneously driven toward maximization. Meanwhile, the priority order of ${\hat{ζ}}_{r}$ will be determined within the optimization objective to better reflect real-world conditions. Therefore, in this section, the ∞-norm of the matrix is employed in the optimization objective to achieve maximization. By approximating the mean of time-varying deviations over their effective time periods, this approach ensures consistency with the form of the original optimization objective. In addition, a reciprocal matrix ${\tilde{T}}_{vd} = {[1 / T_{vd, 1}, 1 / T_{vd, 2}, \dots, 1 / T_{vd, r}]}^{Τ}$ of the effective time periods for multiple uncertainties is constructed, and the improved optimization objective is then formulated as follows:

max 𝜿 ̂ 𝜁 CM = ‖ 𝑯 ̂ ‖ \infty, 𝑯 ̂ = 𝜻 ̂ ⊙ 𝑻 ˜ vd

(59)

where ζ_CM is the actual classified mean deviation of the uncertainty variables, and $\hat{H}$ is the time-varying mean judgment matrix of multiple uncertainties.

Therefore, the system optimization model based on SOI-IGDT is specifically formulated as follows:

⎧⎩⎨ max 𝜿 ̂ 𝜁 CM = ‖ 𝑯 ̂ ‖ \infty, 𝑯 ̂ = 𝜻 ̂ ⊙ 𝑻 ˜ vd s . t . Pr {𝐹 (𝜿 ̂, 𝝍) \leq 𝐹 ac} \geq 𝜏 ℎ ̂ (𝜿 ̂, 𝝍) = 0, 𝑔 ̂ (𝜿 ̂, 𝝍) \leq 0 𝝍 \in 𝛷 ̂ (𝝍 ˜, 𝜏) = {𝝍 | 𝝍 - (𝝍 ˜, 𝜏) \leq 𝝍 \leq 𝝍 + (𝝍 ˜, 𝜏)} (𝒖 ̂ + - 𝒖 ̂ -) 𝝍 + (𝒖 ̂ + - 𝒖 ̂ -) 𝝍 ˜ = 𝜻 ̂ ⊙ 𝝍 ˜ 𝝌 ̂ = sign (𝝍 - 𝝍 ˜), 𝒖 ̂ + = | 𝝌 ̂ | + 𝝌 ̂ 2, 𝒖 ̂ - = | 𝝌 ̂ | - 𝝌 ̂ 2

(60)

3.3. System Optimization Scheduling Model Based on SOI-IGDT

Based on the aforementioned basic structure and optimization model of EHUIES, this section establishes an EHUIES optimization scheduling model considering multiple uncertainties, based on SOI-IGDT. The model primarily includes the optimization objective and the corresponding constraints.

(1): Optimization Objective

Combining the multi-uncertainty analysis in Section 2.3, $ψ$ considers the uncertainties and corresponding stochastic characteristics of wind power, load, and combustion reactions. Additionally, the effective time period for wind power, load, and combustion reactions is set as T, i.e., $T_{vd, r} = T$ . Therefore, the system model’s optimization objective can be expressed as follows:

max 𝜿 ̂ 𝜁 CM = ‖ 𝑯 ̂ ‖ \infty, 𝑯 ̂ = 𝜻 ̂ ⊙ 𝑻 ˜ vd, 𝜻 ̂ = ⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢ 𝜻 ̂ 𝑤 𝜻 ̂ UD load 𝜻 ̂ US load, 𝑙 𝜻 ̂ co, 𝑧 ⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

(61)

where ${\hat{ζ}}_{w}$ , ${\hat{ζ}}_{load}^{UD}$ , ${\hat{ζ}}_{load, l}^{US}$ , and ${\hat{ζ}}_{co, z}$ represent the time-varying deviation matrices of the actual power output of the w-th wind farm, the urban distribution network load, the load of the l-th user unit, and the time-varying deviation matrix of the actual steady-state combustion reaction margin of gas component z, respectively.

(2): Constraints

By setting the confidence level τ, more accurate and reasonable scheduling strategies can be formulated for the EHUIES optimal operation. Here, $F_{ac}$ represents the deterministic model solution results, and $ψ = {[P_{w}, P_{load}^{UD}, P_{load, l}^{US}, ϑ_{sm, z}]}^{Τ}$ is established as time-dependent with respect to the scheduling horizon. Accordingly, the system’s total operation cost risk assessment constraint based on classified confidence can be expressed as follows.

Pr {𝐹 (𝜿 ̂, 𝑷 𝑤, 𝑷 UD load, 𝑷 US load, 𝑙, 𝝑 sm, 𝑧) \leq 𝐹 ac} \geq 𝜏

(62)

The boundary conditions for multiple uncertain variables can be specifically expressed by the following formula:

⎧⎩⎨ (𝒖 ̂ + - 𝒖 ̂ -) 𝝍 + (𝒖 ̂ + - 𝒖 ̂ -) 𝝍 ˜ = 𝜻 ̂ ⊙ 𝝍 ˜, 𝝍 ˜ = [𝑷 ˜ 𝑤, 𝑷 ˜ UD load, 𝑷 ˜ US load, 𝑙, 𝝑 ˜ sm, 𝑧] Τ 𝝌 ̂ = sign (𝝍 - 𝝍 ˜), 𝒖 ̂ + = | 𝝌 ̂ | + 𝝌 ̂ 2, 𝒖 ̂ - = | 𝝌 ̂ | - 𝝌 ̂ 2

(63)

In addition, other related EHUIES constraints also include those from Equations (11) to (50). In summary,

\forall τ_{c} \in [{\hat{τ}}_{\min}, {\hat{τ}}_{\max}]

c \in ℤ_{+}

, the EHUIES optimization scheduling process based on SOI-IGDT is shown in Figure 5.

Figure 5. Flowchart of optimal scheduling of EHUIES based on SOI-IGDT.

3.4. Model Solution of SOI-IGDT Considering Exergy Boost

The previous sections have introduced the definition of exergy efficiency, the formulation of the SOI-IGDT framework, and its application to the scheduling optimization model of EHUIES under uncertainty. Building upon this foundation, this section aims to incorporate an exergy enhancement mechanism to further improve the system’s thermodynamic performance. Specifically, a penalty-based approach is designed to embed exergy efficiency improvement into the original model. The motivation for introducing exergy-based optimization lies in the need to improve the quality of energy utilization in integrated energy systems. While traditional scheduling approaches often focus on energy quantity and cost minimization, they tend to overlook energy grade differences among electricity, heat, cold, and gas. Enhancing exergy efficiency allows for a more refined and effective use of energy resources, which is essential for achieving high-performance, low-carbon multi-energy system operation.

To achieve efficient utilization of electricity, cooling, heating, HCNG, and other energy forms in EHUIES, this section introduces an exergy efficiency evaluation indicator into the aforementioned model to further enhance the system’s exergy. However, since the existing SOI-IGDT does not effectively facilitate this process, improvements to the original SOI-IGDT are necessary. The existing approach for handling the exergy efficiency optimization objective under multi-objective optimization frameworks based on IGDT and stochastic optimization incorporates the maximization of exergy efficiency from the deterministic model into the constraints. Specifically, it ensures that, under different confidence levels, the exergy efficiency

{\hat{η}}_{x}

is not lower than the deterministic model’s exergy efficiency

{\hat{η}}_{ex, ac}

. Specifically, when

\forall τ_{c} \in [{\hat{τ}}_{\min}, {\hat{τ}}_{\max}]

, it follows that

{\hat{η}}_{ex} \geq {\hat{η}}_{ex, ac}

holds, thereby ensuring that the system exergy efficiency does not fall below the expected target value. It is evident that the feasible solution space for exergy efficiency optimization depends on both

{\hat{η}}_{ex, ac}

and

F (\hat{κ}, ψ)

. However, since the primary optimization objective is to maximize the deviation margin, this does not necessarily guarantee condition

{\hat{η}}_{ex} ≫ {\hat{η}}_{ex, ac}

, and the resulting

{\hat{η}}_{ex}

will represent a compromised solution after system optimization. Consequently, in the process of pursuing maximum exergy enhancement, the aforementioned method may exhibit randomness and limitations. Based on the aforementioned limitations, the following improvement strategies are proposed:

(1): The exergy efficiency is incorporated as an evaluation metric within the constraint framework, replacing the expected target value ${\hat{η}}_{ex, ac}$ derived from the deterministic model with a new assessment target ${\hat{η}}_{ex, st}$ . During the optimization process, the discrepancy between the actual exergy performance ${\hat{η}}_{ex}$ and the target value ${\hat{η}}_{ex, ac}$ is analyzed. A penalty function is constructed in the form of a first-order term using the method of Lagrange multipliers, while the second-order term is neglected to reduce computational complexity. Through iterative updates, this penalty mechanism gradually minimizes the deviation between ${\hat{η}}_{ex}$ and ${\hat{η}}_{ex, ac}$ , thereby enhancing the overall system exergy efficiency ${\hat{η}}_{x}$ .
(2): Due to the interdependent relationship between ${\hat{η}}_{ex}$ and $F (\hat{κ}, ψ)$ within the system, enhancing ${\hat{η}}_{x}$ induces corresponding adjustments in $F (\hat{κ}, ψ)$ , depending on the variation direction of ${\hat{η}}_{x}$ during the iterative process. Therefore, in the constraint formulation of $F (\hat{κ}, ψ)$ after equivalent deterministic transformation, a variation multiplier is introduced to dynamically track the iterative evolution of ${\hat{η}}_{ex}$ .

Firstly, the energy utilization within all regions of the EHUIES is analyzed, which primarily consists of the urban distribution network, the HCNG network, and urban terminal user units. To ensure consistency in the representation of variables, all energy forms are expressed in power units. A detailed examination of the input–output energy forms across different EHUIES regions reveals that the input side mainly includes natural gas sources, wind farms, photovoltaic power stations, user-side PV generation, electricity purchased by users, and discharging from electrical energy storage systems. The output side primarily comprises urban grid electrical loads, user-side electrical loads, energy storage charging, industrial park cooling loads, industrial park heating loads, and HCNG network gas loads. Accordingly, the exergy efficiency of the system can be expressed as follows:

𝜂 ̂ ex = \sum 𝑡 = 1 𝑇 ⎡⎣⎢ 𝜛 e 𝑃 UD load , 𝑡 + \sum 𝑙 = 1 𝐿 𝜛 e 𝑃 US load , 𝑙 , 𝑡 + \sum 𝑢 = 1 𝑈 𝜛 e 𝑃 ch , 𝑢 , 𝑡 </

Related Articles

Contact us