Low-Carbon Operational Optimization of Integrated Electricity-Heat-Gas Systems Considering Regulation Capability

1 Introduction

1.1 Motivation

As of 2024, China’s renewable energy generation reached approximately 3.47 trillion kilowatt-hours, accounting for around 35% of the country’s total electricity production [1]. With the development of new energy such as wind and solar energy, limited by location factors, wind and solar resources are mainly concentrated in the three northern regions. The load demand is concentrated in the coastal areas of central and southeastern China, leading to difficulties in consumption and prominent issues of wind and solar power abandonment [2]. At the end of 2024, the National Development and Reform Commission and the National Energy Administration jointly issued the Implementation Plan for the Special Action on Enhancing Power System Regulation Capability (2025–2027) [3]. The plan outlines a series of key tasks aimed at improving the planning and development of regulation resources, optimizing their dispatch mechanisms, and refining their participation in market-based systems. By strengthening the regulatory flexibility and dispatchability of the power system, the initiative seeks to support the effective integration of over 200 GW of newly installed renewable energy capacity annually during 2025–2027. It sets a clear target for 2027: the national renewable energy utilization rate shall not fall below 90%. However, relying solely on the power system’s inherent regulation capability is insufficient to achieve the local integration of large-scale renewable energy. It is therefore imperative to establish a multi-energy system encompassing diverse energy forms, enabling coordinated operation between the power system and other energy infrastructures [4]. Hydrogen, as a clean and efficient secondary energy source, offers high energy density, large capacity, long lifespan, and ease of storage and transmission. By coupling and complementing different forms of energy, such as electricity, heat, cold, and gas—a hydrogen-containing coupled energy system can be established. This system can significantly improve energy utilization efficiency and promote clean energy consumption [5]. Consequently, coupling systems integrating multiple energy sources like electricity, heat, and gas have become a key focus of research.

1.2 Literature Review

In the electric-thermal-gas coupling system, electrolytic cells and hydrogen fuel cells serve as key conversion and coupling equipment, necessitating detailed research on their operation. Regarding electrolytic cells, Ma et al. [6] constructed a physical model of a PEM electrolyzer by simulating the performance of coupled thermal fluid, matter transfer and electrochemical processes. Dang et al. [7] compared and summarizes the theoretical mathematical models of PEM electrolyzers and other types of electrolyzers. It analyzes the impact of various control strategies and shutdown procedures on transient behavior, and explains the transient phenomena observed during the shutdown process based on the theoretical principles of hydrogen-oxygen fuel cells. Frensch et al. [8] provided impedance spectra and current-voltage characteristic curves of electrolytic cells under various current densities and over time, indicating that the voltage efficiency of electrolytic cells continuously decreases with operating time.

Concerning hydrogen fuel cells, Yang et al. [9] proposed accounting for the cogeneration characteristics of hydrogen fuel cells to mitigate the limitations of cogeneration units whose electricity output is heat-determined. Han et al. [10] highlighted a significant response delay when hydrogen fuel cells start at room temperature. From this literature review, it can be concluded that the efficiency of electrolytic cells dynamically varies with external conditions and component aging, while hydrogen fuel cells exhibit delayed response characteristics. However, existing research on optimizing hydrogen-containing energy systems typically assumes a fixed efficiency for electrolytic cells, while allowing hydrogen fuel cells to respond in real-time, thus neglecting both the dynamic efficiency changes of electrolytic cells and the delayed response of hydrogen fuel cells.

With the increasing coupling characteristics of multiple loads such as electricity, heat, cooling, and gas, traditional demand response has evolved into a comprehensive demand response. Dong et al. [11,12] developed an operational optimization model for comprehensive energy systems considering a single demand response type. Building upon this, Wang et al. [13] proposes a synergistic planning method for an integrated energy system with hydrogen storage taking into account the coupled use of electric-thermal energy, which effectively reduces the system carbon emission and improves the comprehensive energy efficiency level. Li et al. [14] established a demand response model incorporating both price-based and incentive-based loads, analyzing the behavioral characteristics of residential participants in demand response programs.

Existing demand response research effectively models various types of responses but still has two main limitations. First, it primarily focuses on electrical and thermal energy systems without simultaneously considering the regulatory capabilities of electrical, thermal, and gas energy systems. In hydrogen energy systems, flexibility arises from the integration of hydrogen fuel cells, electrolytic cells, and hydrogen blending pipelines, allowing the system to act as a schedulable resource. Second, demand response pricing generally adopts fixed call prices, which inadequately reflect the system’s willingness to adjust the demand scale.

Optimizing operation to maximize the benefits of electric, thermal, and gas coupling systems remains a key research area. Wang et al. [15,16] studied a coupled system incorporating renewable distributed generation and multiple energy loads. To mitigate uncertainty in renewable output, they developed an operational optimization model for participation in day-ahead and real-time electricity markets. Extending beyond single-market participation, Wang et al. [17] proposed a day-ahead economic scheduling model for a Renewable Energy-dominated Hybrid Integrated Electricity-Hydrogen Energy System (REHIES), aiming to accommodate market demand through enhanced coordination between electricity and hydrogen energy. Li et al. [18] proposed an innovative Electro-Hydrogen Regional Energy System and developed a cost-minimizing scheduling model that incorporates a stepped carbon pricing mechanism and the operational characteristics of hydrogen fuel cells, aiming to enhance system flexibility and renewable energy utilization.

While existing research can optimize system operation across multiple markets through multi-energy coupling, it generally fails to consider the coupling correlation between prices across these markets.

1.3 Contributions of the Paper

On the basis of the above research, this article conducts a low-carbon operation optimization study of the electric-heat-gas coupling system (EHGCS) that incorporates regulation capability. Compared with existing studies, this article presents the following innovations:

• (1)Device-Level Dynamic Modeling

Unlike previous works that assume constant efficiency or instantaneous response of hydrogen-related components, this study incorporates a dynamic hydrogen production efficiency model for proton exchange membrane (PEM) electrolyzers, accounting for environmental variations and degradation over time. Additionally, a delay-response penalty cost is introduced for hydrogen fuel cells, reflecting the temporal inertia in electricity-to-gas-to-electricity conversion. This dual modeling approach enhances the accuracy and responsiveness of system-level optimization.

• (2)Coupled Regulation and Pricing Mechanism

While most prior studies focus on single energy carriers or assume static regulation capacity, this work constructs a regulation capability model that captures the coupling effects among electric, thermal, and gas subsystems. By integrating the operational flexibility of hydrogen pipelines, electrolyzers, and fuel cells, the model enables a more comprehensive and realistic quantification of system regulation potential. Furthermore, a dynamic pricing mechanism based on multi-energy contribution coefficients is proposed to reflect real-time supply–demand interactions, which has been largely overlooked in existing literature.

• (3)Market-Aware Optimization with Price Dependency

Existing models often treat energy and carbon markets independently or adopt simplistic correlation assumptions. In contrast, this study applies a kernel density estimation–copula approach to quantify nonlinear dependencies between electricity and carbon prices. This allows the optimization framework to dynamically respond to market fluctuations and improve resource allocation under joint market signals, leading to a more robust and profitable low-carbon operation strategy.

1.4 Paper Organization

The remainder of this paper is organized as follows. Section 2 presents the design of the operational framework for the EHGCS and constructs the equipment output model, as well as a regulation capability model that accounts for the coupling of electrical, thermal, and gas energy subsystems across multiple systems. Section 3 builds upon this foundation to develop a low-carbon operation optimization model based on the kernel density copula method, with the objectives of maximizing EHGCS net profit and clean energy consumption. Section 4 provides a numerical analysis using a case study of a specific park, offering practical insights for optimizing the operation of multi-energy coupling systems.

2 EHGCS Framework Design and Output Modeling

2.1 System Operation Framework

The EHGCS described in this article consists of coupled subsystems of electrical energy, thermal energy, and gas energy. The electricity subsystem meets load demand through wind power generation, photovoltaic power generation, batteries, and hydrogen fuel cells. The gas energy subsystem comprises two types: hydrogen energy and natural gas. Hydrogen energy is primarily supplied by electrolyzer units and hydrogen storage tanks to satisfy hydrogen load demand. Natural gas is partially produced via carbon capture equipment that captures carbon dioxide from gas boilers and hydrogen from hydrogen storage tanks. This captured gas is converted into natural gas through methanation equipment to meet natural gas load demand. Additionally, to balance economic and low-carbon objectives, surplus hydrogen energy can be transported via natural gas pipelines to natural gas plants and hydrogen refueling stations in other regions. The heat energy subsystem supplies thermal load by utilizing waste heat generated from electrolyzers, hydrogen fuel cells, and gas boilers. The overall system framework is illustrated in Fig. 1.

Figure 1: System coupling framework.

2.2 Equipment Output Model

Based on the system coupling framework, model the equipment involved in the electric, gas, and heat energy subsystems.

2.2.1 Electric Energy Subsystem

• (1)Wind power and photovoltaic

The output of wind power generation in the electrical energy subsystem is closely related to wind speed, while photovoltaic power generation output depends on light intensity. The Weibull distribution and Beta distribution are commonly used to model the outputs of wind power and photovoltaic power, respectively [19]. This is expressed as shown in Eq. (1): fptWT=α2WTα1WTptWTα1WTα2WT−1e∫−ptWTα1WTα2WTfptpv=Γα1pv+α2pvΓα1pvΓα2pvptpvpmaxpvα1pv−11−ptpvpmaxpvα2pv−1(1) where p t WT and p t pv represents the output of wind power generation and photovoltaic power generation at time t, respectively; α 1 WT and α 2 WT represent a Weibull distribution parameter; α 1 WT denotes a scale parameter; α 2 WT indicates a shape parameter; α 1 pv and α 2 pv mean the Beta distribution parameter for photovoltaic power generation; and p max pv is the maximum output of photovoltaic power generation.

• (2)Battery (BA)

The power state values of the battery in the electric energy system should remain equal throughout a cycle of use, as shown in Eq. (2). The operating life of a battery is an important parameter, and its operating life model is specifically shown in Eq. (3): Esoc,tba=Esoc,t−Δtba+θchabaPcha,tbaΔt−Pdis,tbaθdisbaΔt(2) Nmaxbad=fd=N0bad−αbaθndba,d=ndbafd(3) where E s o c , t b a and E s o c , t − Δ t b a represent the power state values of the battery at time t and time t − Δ t , respectively; θ c h a b a and θ d i s b a are the charging and discharging efficiencies of the battery, respectively; P c h a , t b a and P d i s , t b a are the charging and discharging power of the battery at time t, respectively; d denotes the depth of charge and discharge of the battery; N max b a d indicates the maximum number of discharges of the battery when the discharge depth is d ; N 0 b a represents the number of charging and discharging times when the battery device is fully charged; α b a denotes a characteristic parameter; and θ n d b a , d is the battery loss rate when the discharge depth is d .

• (3)Hydrogen fuel cell

PEM hydrogen fuel cells offer advantages such as cleanliness, zero pollution, and high energy conversion efficiency. Therefore, in this article, the term “hydrogen fuel cells” specifically refers to PEM hydrogen fuel cells. When the fuel cell parameters and external environment remain constant, the power output can be indirectly controlled by adjusting the hydrogen input rate and the circulating water flow rate. The output power and efficiency of hydrogen fuel cells are detailed in Eq. (4): PtHFC=2FvH2HFCMH2Vnst−Vact−Vom−VncθtHFC=PtHFC/Qt,inHFC·θtG2E(4) where P t H F C , Q t , i n H F C and θ t H F C are the power generation, hydrogen consumption, and efficiency of the hydrogen fuel cell at time t, respectively; θ t G 2 E is the hydrogen to electricity conversion efficiency; F is the Faraday coefficient; v H 2 H F C is the rate at which hydrogen gas flows into the hydrogen fuel cell; M H 2 is the molar mass of hydrogen gas is the molar mass of hydrogen gas; and V n s t , V a c t , V o m and V n c represents Nernst voltage, activation overvoltage, Ohmic overvoltage, and concentration overvoltage, respectively; Reference [20] provides a comprehensive electrochemical model based on experimental data, which quantifies the contributions of various overpotentials under different operating conditions. Their work serves as the theoretical foundation and parameter reference for constructing the voltage model adopted in this study, as shown in Eq. (5): Vnst=1.229−0.85×10−3×Tin−Tno+4.3×10−5×TinInρH2+0.5ρo2Vact=α1HFC+α2HFCTin+α3HFCTinInOo2+α4HFCTinIndenVom=denRinVnc=0.016In1−den/denmax(5) where T i n and T n o represent the internal stack temperature and standard temperature of hydrogen fuel cells, respectively; ρ H 2 and ρ o 2 denote hydrogen partial pressure and oxygen partial pressure, respectively; α 1 H F C , α 2 H F C , α 3 H F C and α 4 H F C are the activation overvoltage function coefficients; O o 2 is the oxygen concentration at the cathode reaction interface; and d e n and d e n max are the current density and current limit density, respectively.

PEM hydrogen fuel cells require a heating process of several tens of minutes from start-up to optimal operating state, with delayed response characteristics. The specific delay coefficient is shown in Eq. (6): Tdd=Tddmin−Tddexeζt+Tddexζ=0.589×10−3−0.709×10−7IHFC+0.554×10−8IHFCInden−0.427×10−9IHFCden(6) where T d d is the temperature of the fuel cell stack; T d d min and T d d e x represent the minimum and optimal operating temperatures of the fuel cell stack, respectively; and ζ is the delay response coefficient of the hydrogen fuel cell.

2.2.2 Gas Energy Subsystem

• (1)Electrolytic cell

The application of PEM electrolytic cells is relatively mature; therefore, the electrolytic cell discussed in this article refers to a PEM electrolytic cell. Hydrogen production depends on the efficiency of the electrolytic cell and the electric-to-hydrogen conversion efficiency. The efficiency of the electrolytic cell is the product of voltage efficiency and current efficiency, as expressed in Eq. (7): QtEC=Pt,inECθtECθtE2GθtEC=θtIθtV(7) where Q t E C and P t , i n E C are the hydrogen production and power consumption of the electrolytic cell at time t, respectively; θ t E C is the efficiency of the electrolytic cell at time t; θ t E 2 G is electric to hydrogen conversion efficiency at time t; and θ t I and θ t V represent the current efficiency and voltage efficiency of the electrolytic cell at time t, respectively.

After long-term operation, PEM electrolytic cells experience component aging and thinning of the PEM. This leads to increased hydrogen permeation and a decline in dynamic current efficiency. Additionally, electrolyte conductivity decreases, causing a rise in Ohmic resistance and an increase in Ohmic voltage. Consequently, the dynamic voltage efficiency of the electrolytic cell decreases as operating time increases. These dynamic changes in current efficiency and voltage efficiency together cause the overall efficiency of the electrolytic cell to vary over time, as described in Eq. (8): θtI=1−2FαdfECαdsECΔρH2ECdenEChM,tECθtV=VtnVtn+aT0ECαanECαchECFIndenECden_anECγan+denEChM_anECδan+hM_neECδne+denEChM0EC2hM,tECb0(8) where α d f E C and α d s E C represent the diffusion coefficient and solubility coefficient of hydrogen gas, respectively; Δ ρ H 2 E C represents the hydrogen partial pressure difference on both sides of the electrolytic cell membrane; d e n E C is the current density of the electrolytic cell; h M , t E C denotes the thickness of the electrolytic cell membrane at time t; V t n is the thermal neutral voltage; T 0 E C indicates the rated operating temperature; α a n E C presents the anode charge transfer coefficient; α c h E C is the stoichiometric coefficient; d e n _ a n E C is the conversion current density of the anode; γ a n is the ratio of anode active area; h M _ a n E C and h M _ n e E C are the thickness of the anode and cathode films, respectively; δ a n and δ n e represents the conductivity of the anode plate and the cathode plate, respectively; h M 0 E C means the initial thickness of the PEM; and b 0 is the initial conductivity of the electrolyte.

• (2)Compressed hydrogen storage system

The compressed hydrogen system consists of a compressor and a hydrogen storage tank. Its operating principle involves compressing hydrogen produced by the electrolytic cell and storing it in the hydrogen storage tank to meet hydrogen demand. However, the compression process consumes power, which is specifically described in Eq. (9): p¯losscm=acmTincmαhrcmαhrcm−1ρoutcmρincmαhrcmαhrcm−1−1ploss,tcm=p¯losscmvtcm3.6×106αwrcmMH2(9) where p ¯ l o s s c m is the power consumption per unit molar mass of hydrogen compressed by the compressor; a c m is the ideal gas constant; T i n c m is the working temperature of the compressor; α h r c m is the heat capacity ratio of compressed gas; ρ i n c m and ρ o u t c m represent the pressure at the inlet and outlet, respectively; p l o s s , t c m is the total power consumption of the compressor at time t; v t c m is the rate at which hydrogen gas enters the compressor at time t; and α w r c m is the working efficiency of the compressor.

Hydrogen, having the smallest relative molecular weight, is more prone to expansion under high pressure compared to other gases. Therefore, accurately characterizing the pressure within hydrogen storage tanks is essential. The relationship between pressure and hydrogen storage capacity is expressed in Eq. (10): ρtHST=fQtHST=QtHSTacmTH2MH2VHST−a1HSTQtHST−a2HSTQtHSTMH2VHST2(10) where ρ t H S T is the pressure of the hydrogen storage tank at time t; Q t H S T is the capacity of the hydrogen storage tank at time t; T H 2 is the gas temperature; V H S T is the volume of the hydrogen storage tank; a 1 H S T and a 2 H S T are the normal form coefficients.

The temporal relationship between the hydrogen storage capacities of two adjacent hydrogen storage tanks is expressed in Eq. (11): QtHST=Qt−ΔtHST+θchaHSTQcha,tHSTΔt−Qdis,tHSTθchaHSTΔt(11) where Q t − Δ t H S T is the capacity of the hydrogen storage tank at time t − Δ t ; θ c h a H S T and θ c h a H S T are the energy charging and discharging efficiencies of the hydrogen storage tank, respectively; and Q c h a , t H S T and Q d i s , t H S T represent the energy of hydrogen storage tank charging and discharging at time t, respectively.

• (3)Methanation device

In the methanation unit, hydrogen gas produced by the electrolyzer or released from the hydrogen storage tank reacts with carbon dioxide captured by the carbon capture equipment to generate natural gas under the action of a catalyst. The consumption of hydrogen and carbon dioxide, as well as the production of natural gas in the methanation unit, are detailed in Eq. (12): Qt,H2md=4Mt,co2mdρH2ρco2αlossmdQt,CH4md=Mt,co2mdρCH4ρco2αlossmd(12) where Q t , H 2 m d is the amount of hydrogen consumed by the methanation unit; M t , c o 2 m d is the mass of carbon dioxide consumed by the methanation unit; ρ H 2 , ρ c o 2 and ρ C H 4 represent the densities of hydrogen, carbon dioxide, and natural gas, respectively; and α l o s s m d is the efficiency of the methanation unit.

• (4)Hydrogen doped natural gas pipeline

Hydrogen-doped natural gas pipelines enable long-distance hydrogen transport without additional investment. Assuming hydrogen is uniformly distributed within the natural gas pipeline, the specific ratio of hydrogen to the calorific value of hydrogen-doped natural gas is expressed in Eq. (13): ηtH2=vtH2/vtH2+vtNGHtmix=ηtH2HH2+1−ηtH2HNG(13) where η t H 2 is the hydrogen doping ratio at time t; v t H 2 and v t N G are the flow rate of hydrogen and natural gases in the natural gas pipeline; H t m i x , H H 2 and H N G represent the calorific value of mixed hydrogen natural, hydrogen, and natural gases, respectively.

2.2.3 Heat Energy Subsystem

• (1)Waste heat utilization

Waste heat utilization refers to recovering waste heat from flue gases of electrolytic cells and hydrogen fuel cells using waste heat recovery devices to meet users’ thermal load demands. The specific amount of waste heat utilized is shown in Eq. (14): HtEC=Pt,inEC·θtEC·θtE2HHtHFC=Qt,inHFC·θtHFC·θtG2HHtrh=θtrhHtEC+HtHFC(14) where H t E C and H t H F C represent the residual heat generated by the electrolytic cell and hydrogen fuel cell at time t, respectively; H t r h is the overall utilization of waste heat at time t; θ t E 2 H is the heat generation efficiency of the electrolytic cell; θ t G 2 H is the heat generation efficiency of hydrogen fuel cells; and θ t r h is the efficiency of waste heat utilization.

• (2)Gas boiler

The gas boiler generates heat energy by consuming natural gas produced by the methanation equipment. The heat output of the gas boiler depends on the natural gas consumption, the boiler’s efficiency, and the calorific value of the natural gas, as shown in Eq. (15): HtGB=Qt,inGBθtGBHNG(15) where H t G B is the heat energy generated by the gas boiler at time t; Q t , i n G B is the natural gas consumed by the gas boiler at time t; and θ t G B is the efficiency of the gas boiler.

• (3)Carbon capture equipment

The carbon capture equipment absorbs carbon dioxide emitted by gas boilers and supplies it to the methanation units for further use. The amount of carbon dioxide captured by the equipment is expressed in Eq. (16): Gco2,tccs=HtGB·g¯co2GB·αccs(16) where G c o 2 , t c c s is the carbon dioxide captured by the carbon capture equipment at time t; g ¯ c o 2 G B is the carbon dioxide emissions per unit from gas boilers; and α c c s is the efficiency of carbon capture equipment.

The energy consumption of the carbon capture system includes both operational and fixed energy consumption, as shown in Eq. (17): Ptccs=Gco2,tccsp¯co2ccs+PBccs(17) where P t c c s is the total energy consumption of the carbon capture equipment at time t; p ¯ c o 2 c c s is the operating energy consumption per unit of carbon dioxide capture; and P B c c s is the fixed energy consumption of carbon capture equipment.

3 Analysis of EHGCS Flexible Regulation Capability

The EHGCS involves multiple energy flows and loads. This article analyzes the flexible regulation capabilities of the electric, heat, and gas energy subsystems, as well as the overall coupled system.

3.1 Electric Energy Subsystem Regulation Capability

In the energy subsystems, regulation capability varies due to users’ subjective willingness and the price difference between peak and off-peak periods. When the price difference is below a certain threshold, users show low enthusiasm for participating in regulation, creating a “dead zone” of regulation participation. Conversely, if the price difference is too large, users tend to reach a saturation point to maintain a comfortable energy consumption experience, resulting in a “saturation zone” for participation. Additionally, the greater the users’ willingness to participate, the higher the likelihood of their engagement in regulation. Based on a combination of logistic and trapezoidal membership functions, a regulation capability model for the electric energy subsystem is constructed. The principle of this model is illustrated in Fig. 2:

Figure 2: Principle diagram of electrical energy subsystem regulation capability.

To account for user participation willingness, the average load transfer rate is calculated based on the load transfer rates of users with high and low willingness to participate. This average rate serves as a measure of the regulation capability of the energy subsystem. Using peak and valley periods as an example, the average load transfer rate during these periods is given by Eq. (18): λ¯fg=λmaxfg+λminfg20≤Δuelefg≤b1fgλminfg+λmaxfg+λminfg21+Δuelefg−b1fgb2fg−b1fgb1fg<Δuelefg≤b2fgλmaxfgΔuelefg>b2fg(18) where λ ¯ f g is the electricity load transfer rate during peak and valley periods; λ max f g and λ min f g are the load transfer rate during peak and valley periods under high and low participation willingness, respectively; b 1 f g is the peak valley difference threshold for the transition from nonparticipating areas to participating areas; b 2 f g is the peak valley difference threshold for the transition from the participating region to the saturated region; and Δ u e l e f g is the peak valley electricity price difference.

Using Eq. (18), the load transfer rates for the peak and valley periods are calculated as a and b, respectively. Based on these values, the regulation capability of the electrical energy subsystem is expressed in Eq. (19): ΔLtele=−λ¯fgLt,fele−λ¯fpLt,felet∈Tfλ¯fpLt,fele−λ¯pgLt,pelet∈Tpλ¯fgLt,fele+λ¯pgLt,pelet∈Tg(19) where Δ L t e l e is the regulation ability of the electric energy subsystem at time t; L t , f e l e and L t , p e l e represent the electricity load demand during peak and normal periods at time t, respectively; and T f , T p and T g denote peak, flat, and valley periods, respectively.

3.2 Heat Energy Subsystem Regulation Capability

Because indoor temperature can fluctuate within a certain range without compromising user comfort, the thermal load can be adjusted accordingly while still meeting comfort requirements. The magnitude of this adjustable thermal load represents the regulation capability of the thermal energy subsystem. The user’s indoor temperature model is shown in Eq. (20): Tt+1id_in=Ttid_ine−ΔtRidτid+RidLtheat+Ttid_out1−e−ΔtRidτid(20) where T t + 1 i d _ i n is the indoor temperature at time t + 1; R i d is the thermal resistance of the building; τ i d is the indoor heat capacity; L t h e a t is the indoor heating power at time t; and T t i d _ i n and T t i d _ o u t are the indoor and outdoor temperature at time t, respectively.

Based on Eq. (20), the user’s heat load model is obtained as shown in Eq. (21):

Let the acceptable range of indoor temperature variation for users be T min i d _ i n , T max i d _ i n , where indoor temperature is increased during low energy consumption periods and decreased during peak energy consumption periods. Based on this, the regulation capability of the thermal energy subsystem is derived by combining with Eq. (21), as shown in Eq. (22): ΔLtheat=Tmaxid_in−Tminid_inRid1−e−ΔtRidτid(22) where Δ L t h e a t is the regulating ability of the thermal energy subsystem at time t.

3.3 Gas Energy Subsystem Regulation Capability

To capture the time-delay characteristics inherent in the electricity-to-gas and gas-to-electricity conversion processes, this study incorporates a response lag parameter within the hydrogen fuel cell modeling. This approach reflects the dynamic delays associated with hydrogen production via electrolysis, storage, and reconversion to electricity. It enhances the representation of the gas subsystem’s dynamic behavior and allows for a more realistic assessment of its temporal effectiveness and stability in supporting power system regulation. The gas energy subsystem can achieve spatial energy transfer via hydrogen-doped natural gas pipelines during periods of low energy consumption. The volume of this spatial transfer represents the regulation capacity of the gas energy subsystem. This regulation capability depends on both the hydrogen blending capacity of the natural gas pipeline and the hydrogen production capacity of the electrolytic cell. The regulation capability model of the gas energy subsystem is detailed in Eq. (23): ΔLtH2=minQmaxEC−LtH2,QmaxGD×20%(23) where Δ L t H 2 is the regulating ability of the gas energy subsystem at time t; Q max E C is the maximum hydrogen production capacity of the electrolytic cell; L t H 2 represents the hydrogen load demand of the user at time t; and Q max G D is the maximum transportation capacity of the pipeline.

3.4 Multi-System Coupling Regulation Capability

The EHGCS enables the conversion between electricity, gas, and heat. Therefore, beyond considering the transfer capabilities of each individual subsystem, it also accounts for the substitutable regulation capability arising from multi-system coupling. This mainly includes electricity-gas load substitution and gas-heat load substitution, as shown in Eqs. (24) and (25): ΔLtE2G=ιE2G·μtele−utgasutgas·L˜adeleΔLtG2E=−θtE2G·ΔLtE2G(24) ΔLtG2H=ιG2H·μtgas−utheatutgas·L˜adheatΔLtH2G=−θtG2H·ΔLtG2H(25) where Δ L t E 2 G and Δ L t G 2 E represent the substitution responses of electrical load and gas load, respectively; ι E 2 G are the elasticity coefficients of electricity load and gas load influenced by price; μ t e l e and u t g a s denote the electricity price and gas price at time t; L ˜ a d e l e is the amount of electricity load that can be replaced before the response; θ t E 2 G is the conversion efficiency between electricity and gas; similarly, Δ L t G 2 H and Δ L t H 2 G represent the substitution responses of gas load and heat load; ι G 2 H is the elasticity coefficient of gas load and heat load affected by price; u t h e a t denotes the heat price L ˜ a d h e a t is the replaceable heat load before response; and θ t G 2 H is the conversion efficiency between gas and heat.

4 Optimization Modeling of EHGCS Low-Carbon Operation Considering Multi-System Regulation

After satisfying its internal diversified load demands, the EHGCS can participate in the joint energy-carbon market. When the system’s free carbon quota exceeds its actual carbon emissions, EHGCS can sell the surplus carbon allowances in the carbon market. Conversely, if emissions exceed the quota, it must purchase additional allowances. Similarly, when EHGCS’s energy supply surpasses demand, it can sell the excess energy in the energy market; otherwise, it can buy the shortfall energy to meet its needs. This interaction forms the joint energy-carbon market. The operating principle of EHGCS’s market participation is illustrated in Fig. 3.

Figure 3: EHGCS participation in market operation schematic.

4.1 Objective Function

To balance the economic and environmental benefits of EHGCS, the objective function aims to maximize both the system’s net profit and clean energy consumption.

4.1.1 Maximizing Net Income

EHGCS net revenue is defined as the difference between total revenue and total cost, as shown in Eq. (26): RnetEHGCS=∑t=1ttotalRtotal,tEHGCS−Ctotal,tEHGCS(26) where R n e t E H G C S represents the net income of EHGCS; R t o t a l , t E H G C S and C t o t a l , t E H G C S represent the total revenue and cost of EHGCS at time t.

• (1)Total revenue

The total revenue of EHGCS includes internal energy sales revenue R s a l e , t E H G C S , market participation revenue R c e m , t E H G C S , and hydrogen blending pipeline gas transmission revenue R h d , t E H G C S , as shown in Eq. (27):

The specific internal sales revenue is shown in Eq. (28): Rsale,tEHGCS=uteleLac,tele+utH2Lac,tH2+utNGLac,tNG+utheatLac,theat(28) where u t H 2 and u t N G are the price of hydrogen and natural gases at time t, respectively; and L a c , t e l e , L a c , t H 2 , L a c , t N G and L a c , t h e a t represent the electricity, hydrogen, natural gas, and heat sold to internal users of the system at time t, respectively.

The specific gas transmission benefits of hydrogen blending pipelines are shown in Eq. (29): Rhd,tEHGCS=Qhd,tH2·uhd,tH2+Qhd,tNG·uhd,tNG(29) where Q h d , t H 2 and Q h d , t N G are the hydrogen transmission rate and natural gas transmission rate of the hydrogen blending pipeline at time t, respectively; and u h d , t H 2 and u h d , t N G represent the hydrogen purchase price and natural gas purchase price of the recipient of the gas transmission at time t, respectively.

The benefits of participating in the market are shown in Eq. (30): Rcem,tEHGCS=Gsale,tcmusale,tcm,co2+Psale,temusale,tem,ele+Qsale,tem,H2usale,tem,H2+Qsale,tem,NGusale,tem,NG+Hsale,temusale,tem,heat(30) where G s a l e , t c m represents the carbon emission rights sold by the system in the carbon market at time t; u s a l e , t c m , c o 2 is the unit carbon emission rights price in the carbon market at time t; P s a l e , t e m , Q s a l e , t e m , H 2 , Q s a l e , t e m , N G and H s a l e , t e m represent the electricity, hydrogen, natural gas, and heat energy sold by the system in the energy market at time t, respectively; and u s a l e , t e m , e l e , u s a l e , t e m , H 2 , u s a l e , t e m , N G and u s a l e , t e m , h e a t are the price per unit of electricity, hydrogen, natural gas, and heat energy in the energy market at time t, respectively.

When the carbon market price rises, EHGCS tends to reduce carbon emissions by lowering gas boiler output, which decreases the system’s energy supply. If the system is in a shortage state, this shortage worsens, increasing market demand and causing energy prices to drop. Conversely, if the system has excess supply, the reduction causes energy prices to rise. This demonstrates a correlation between carbon market prices and energy market prices. In this paper, the kernel density estimation method-Copula method is used to characterize the nonlinear correlation between carbon market price and energy market (price). Firstly, the marginal distribution functions of various energy prices and carbon prices are obtained through the kernel density estimation method, so as to preserve their actual distribution characteristics; subsequently, the Copula function is introduced to construct a joint distribution model to establish the dependence structure between different energy prices and carbon prices. This joint distribution is used to simulate the joint fluctuation of various price variables in the scenario generation stage, realizing the multi-scenario inputs under the price linkage between carbon and energy markets. On this basis, the generated price scenarios are embedded into the EHGCS low-carbon operation optimization model, which is used to calculate the net benefit and uncertainty cost under multiple scenarios, and to control the tail risk through the CVaR indicator, so as to improve the model’s adaptability and robustness to market fluctuations. The historical sample vectors of carbon emission rights, electricity, hydrogen, natural gas, and thermal energy prices as u ˜ s a l e , 1 c m , c o 2 , u ˜ s a l e , 2 c m , c o 2 , ⋯ , u ˜ s a l e , t c m , c o 2 , ⋯ , u ˜ s a l e , t t o t a l c m , c o 2 , u ˜ s a l e , 1 e m , e l e , u ˜ s a l e , 2 e m , e l e , ⋯ , u ˜ s a l e , t e m , e l e , ⋯ , u ˜ s a l e , t t o t a l e m , e l e , u ˜ s a l e , 1 e m , H 2 , u ˜ s a l e , 2 e m , H 2 , ⋯ , u ˜ s a l e , t e m , H 2 , ⋯ , u ˜ s a l e , t t o t a l e m , H 2 , u ˜ s a l e , 1 e m , N G , u ˜ s a l e , 2 e m , N G , ⋯ , u ˜ s a l e , t e m , N G , ⋯ , u ˜ s a l e , t t o t a l e m , N G and u ˜ s a l e , 1 e m , h e a t , u ˜ s a l e , 2 e m , h e a t , ⋯ , u ˜ s a l e , t e m , h e a t , ⋯ , u ˜ s a l e , t t o t a l e m , h e a t , respectively. The kernel density estimation method is used to estimate the cumulative distribution function of these multivariate prices, as shown in Eq. (31): fusalecm,co2=1n????∑t=1ttotalusalecm,co2−u˜sale,tcm,co2????fusaleem,ele=1n????∑t=1ttotalusaleem,ele−u˜sale,tem,ele????fusaleem,H2=1n????∑t=1ttotalusaleem,H2−u˜sale,tem,H2????fusaleem,NG=1n????∑t=1ttotalusaleem,NG−u˜sale,tem,NG????fusaleem,heat=1n????∑t=1ttotalusaleem,heat−u˜sale,tem,heat????(31) where n is the smoothing coefficient; and ???? is the nuclear density.

• (2)Total cost

The total cost of EHGCS comprises of system uncertainty cost C r i s k , t E H G C S , unit operation cost C op , t E H G C S , regulation cost C DR , t E H G C S , delayed response penalty cost of hydrogen fuel cell C pun , t E H G C S , market operation cost C cem , t E H G C S , and transportation cost of hydrogen doped natural gas channel C hd , t E H G C S , as shown in Eq. (32):

In this study, uncertainty costs are characterized using Conditional Value-at-Risk (CVaR), a widely adopted risk measure that effectively captures the impact of extreme deviations caused by load fluctuations and renewable generation uncertainty. As a second-order risk metric, CVaR incorporates tail risks into the objective function by constraining the expected loss beyond a specified confidence level, thereby enhancing the model’s robustness and its ability to manage volatility. This enables the proposed optimization framework to explicitly account for the stochastic nature of renewable energy and demand. The detailed formulation can be found in reference [21]. The specific operating cost of the unit is shown in Eq. (33): Cop,tEHGCS=∑j∈JPtju¯opj+∑s∈SQtsu¯ops+HtGBu¯opGB+Gco2,tccsu¯opccs(33) where P t j and u ¯ o p j are the output and unit operating costs of the j-th device in the electrical energy subsystem at time t, respectively; J is a collection of partial equipment in the electrical energy subsystem, J ∈ W T , p v , b a , H F C ; Q t s and u ¯ o p s are the output and operating cost of the s-th device in the gas energy subsystem at time t, respectively; S is a collection of partial equipment in the gas energy subsystem, S ∈ E C , H S T , m d ; u ¯ o p G B and u ¯ o p c c s are the unit operating costs of gas boilers and carbon capture equipment.

When the system experiences a significant energy demand shortage, the adjustment price should increase accordingly. To reflect this, a contribution coefficient is incorporated into the adjustment cost to form a dynamic adjustment price, as shown in Eq. (34): CDR,tEHGCS=∑w∈WΔL⌢tw·u¯tDR,wu¯tDR,w=minc¯tDR,w×εtwmaxεtw×ϑDR,usale,tem,w(34) where Δ L ⌢ t w and u ¯ t D R , w are the adjustment amount and unit dynamic adjustment price for the w-th type of load, respectively; W is the load set, W ∈ e l e , H 2 , N G , H e a t ; c ¯ t D R , w is the adjustment cost for the w-th class of load users; ε t w is the supply-demand deviation at time t; ϑ D R is the benchmark contribution coefficient; and u s a l e , t e m , w is the selling price of the w-th load in the energy market.

The delayed response cost of hydrogen fuel cells is specifically represented in Eq. (35): Cpun,tEHGCS=ς·PtHFC·u¯pun,tHFC(35) where u ¯ p u n , t H F C is the penalty cost per unit of delayed response.

The specific market operating costs are shown in Eq. (36): Ccem,tEHGCS=Gbuy,tcmubuy,tcm,co2+Pbuy,temubuy,tem,ele+Qbuy,tem,H2ubuy,tem,H2+Qbuy,tem,NGubuy,tem,NG+Hbuy,temubuy,tem,heat(36) where G buy , t c m represents the carbon emission rights purchased by the system in the carbon market at time t; u buy , t c m , c o 2 is the purchase price of unit carbon emission rights in the carbon market at time t; P buy , t e m , Q buy , t e m , H 2 , Q buy , t e m , N G and H buy , t e m represent the electricity, hydrogen, natural gas, and heat energy purchased by the system in the energy market at time t, respectively; and u buy , t e m , e l e , u buy , t e m , H 2 , u buy , t e m , N G and u buy , t e m , h e a t represent the purchase prices of unit electricity, hydrogen, natural gas, and thermal energy in the energy market at time t.

4.1.2 Maximizing Clean Energy Consumption Rate

The maximization of the clean energy consumption rate aims to maximize the absorption of wind and solar power output within the system, as shown in Eq. (37): maxϕEHGCS=max∑t=1ttotalPtac,WT+Ptac,pvPtWT+Ptpvttotal(37) where ϕ E H G C S is the absorption rate of EHGCS while P t a c , W T and P t a c , p v represent the actual consumption of wind and photovoltaic power generation at time t, respectively.

4.2 Constraint Condition

The operational constraints of EHGCS primarily include unit operation constraints, power balance constraints, regulation capacity constraints, and transmission channel constraints.

4.2.1 Unit Operation Constraints

The primary constraint on unit operation is that the output of each piece of equipment must remain within its specified upper and lower limits, as shown in Eq. (38). Pminj≤Ptj≤Pmaxjj∈JQmins≤Qts≤Qmaxss∈SHminGB≤HtGB≤HmaxGBGco2,minccs≤Gco2,tccs≤Gco2,maxccsρ¯mincm≤ρoutcmρincm≤ρ¯maxcm(38) where P min j and P max j are the minimum and maximum output of the j-th device in the electrical energy subsystem, respectively; Q min s and Q max s are the minimum and maximum output of the s-th device in the hydrogen energy subsystem, respectively; H min G B and H max G B are the minimum and maximum output of the gas boiler, respectively; G c o 2 , min c c s and G c o 2 , max c c s represent the minimum and maximum values of carbon capture by carbon capture equipment, respectively; and ρ ¯ min c m and ρ ¯ max c m represent the minimum and maximum values of the compressor boost ratio, respectively.

4.2.2 Regulatory Ability Constraints

The regulation capacity constraint refers to the requirement that the actual regulation performed by multiple subsystems must not exceed their respective maximum regulation capacities, as shown in Eq. (39).

4.2.3 Transmission Channel Constraints

The transmission channel constraint refers to the upper and lower limits on both the energy exchanged between the EHGCS and the external market and the volume of hydrogen transported through hydrogen-doped natural gas pipelines, as shown in Eq. (40). Pminem≤Psale,tem/Pbuy,tem≤PmaxemQminem,H2≤Qsale,tem,H2/Qbuy,tem,H2≤Qmaxem,H2Qminem,NG≤Qsale,tem,NG/Qbuy,tem,NG≤Qmaxem,NGHminem≤Hsale,tem/Hbuy,tem≤Hmaxem0≤Qhd,tH2≤QmaxGD×20%0≤Qhd,tNG≤QmaxGD(40) where P min e m and P max e m represent the minimum and maximum values of electricity transmission in the energy market, respectively; Q min e m , H 2 and Q max e m , H 2 represent the minimum and maximum values of hydrogen energy transfer in the energy market, respectively; Q min e m , N G and Q max e m , N G represent the minimum and maximum values of natural gas transmission in the energy market, respectively; and H min e m and H max e m represent the minimum and maximum values of heat energy transfer in the energy market, respectively.

4.2.4 Power Balance Constraint

Power balance constraint refers to the power balance of the electrical, thermal, and gas subsystems, as shown in Eq. (41):