# Joint Optimization of Planning and Operation for Integrated Electrothermal and Hydrogen Coupling System

Authors: Xiaoou Liu
Journal: International Journal of Energy Research
 (2023)

---

## Abstract

In view of the rough modeling of hydrogen energy system, insufficient analysis of hydrogen energy utilization mode, and the lack of consideration of the multiple uncertain factors, this paper proposes a joint optimization method of planning and operation for integrated electrothermal and hydrogen coupling system (IEHS) based on multiverse optimizer algorithm (MVO) and interval linear programming method. Firstly, the frame structure of IEHS is constructed, and multiple utilization modes of hydrogen energy are analyzed. Secondly, considering the hydrogen physical characteristics, the operation model of high-pressure hydrogen storage device based on Van Deemter equation is established, and the linear model of compressor power consumption is built considering dynamic compression ratio. On this basis, the bilevel optimization model of planning and operation for IEHS is established considering multiuncertainty and multiple hydrogen energy utilization modes. The upper model realizes the capacity configuration of equipment with the goal of minimizing the total construction and operation cost within the service life. The lower model optimizes the output of each equipment with the goal of minimizing the operation cost. The example results show that MVO and interval linear programming method for solving the bilevel optimization model have the better convergence speed and convergence ability. Further, the multiple utilization modes of hydrogen energy can effectively reduce the construction and operation cost of IEHS, and accurately depicting the physical characteristics of the hydrogen storage system can increase the coordination and optimization ability of IEHS. Considering multiple uncertain factors, the optimal configuration scheme and operation cost of IEHS are more stable.


---

## Body

## 1. Introduction

In modern society, the energy crisis and environmental pollution are becoming more and more serious, and the transition of energy structure and consumption is imminent. As the most promising secondary energy in the 21st century, hydrogen energy (H2) is an important means to realize the low-carbon transformation of power system, such as in Refs. [1, 2]. Hydrogen energy can realize flexible conversion among various energy sources in Ref. [3] and has flexible regulation characteristics and great application potential in Ref. [4]. Hydrogen is produced through water electrolysis with renewable energy that can not only meet the end hydrogen energy demand but also realize the seasonal storage of large-scale energy in Ref. [5] and make significant economic value in Refs. [6, 7]. With the hydrogen energy utilization in a large scale and the increase of renewable energy penetration level, coupled electricity and hydrogen energy system has become a research hotspot in the current energy field, such as in Ref. [8]. Integrated energy system (IES) can realize the complementary advantages and cascade utilization of different energy sources through the mutual conversion of multiple energy sources and has outstanding effects in improving energy efficiency and the capacity of renewable energy utilization. Therefore, the capacity configuration and operation scheduling of integrated electrothermal and hydrogen coupling system (IEHS) are the key basis to improve the technical and economic benefits of IES, which has important theoretical research value and broad engineering application prospects.At present, the researches on the optimal operation of IEHS are mainly realized by the hydrogen fuel cell (HFC) mode. Gabrielli et al. [9] proposed a seasonal hydrogen storage model and evaluated the impact of HFC on the mid- and long-term operation of the system. Leander et al. [10] use electrolyzer, seasonal hydrogen storage device, and HFC to promote the capacity of renewable energy utilization and reduce the reserve capacity of equipment. However, the HFC mode has multiple energy conversion links, and the comprehensive efficiency is only 15%~22%, which can be deduced from Refs. [11, 12]. It is difficult to improve the efficiency of IEHS only by relying on the HFC mode. Therefore, it is necessary to deeply explore the multiple utilization modes of hydrogen energy and their correlation. Multiple hydrogen energy utilization modes are widely used, such as hydrogen injection in nature gas network (HING) in Refs. [13, 14], hydrogen energy supplied to hydrogen stations in Refs. [15, 16], cogeneration of HFC, and seasonal hydrogen storage. It further deepens the coupling relationship between hydrogen energy system, power system, nature gas system, and thermal system. How to optimize the multiple utilization modes of hydrogen energy and improve the operation efficiency of IEHS has become a key problem to be solved.Hydrogen storage device (HSD) is the key equipment to realize the mid- and long-term scheduling of IEHS. The accurate simulation of HSD operation characteristics is an important basis for the optimal operation of IEHS. Because hydrogen is the material with the lowest relative molecular mass, hydrogen physical characteristics have a significant impact on the hydrogen storage system. On the one hand, existing researches often describe the relationship between pressure and mass of HSD based on the perfect gas equation, such as in Refs. [17, 18]. The hydrogen physical characteristics under different gas pressures are quite different, and the deviation between the state of high-pressure hydrogen and the ideal gas state is up to 20%, such as in Ref. [19]. Ignoring this deviation will be difficult to accurately evaluate the safe operation state of HSD. On the other hand, hydrogen compressor requires more energy than other gases, such as in Ref. [20]. However, the existing researches often estimate the compressor power consumption based on fixed compression ratio, which is difficult to accurately calculate the operation efficiency of HSD.In addition, the uncertainty of load and renewable energy will have a negative impact on the capacity configuration of IEHS. When the uncertainty of renewable energy, load, and energy price are fully considered in the capacity configuration, the commonly used modeling methods mainly include random method, fuzzy method, robust method, and interval method. Zhou et al. [21] describe the uncertainty of load and renewable energy by using interval linear programming and random constraint programming, respectively, and establish an economic scheduling model based on interval linear programming with random constraint. Yao and Wang [22] consider the impact of renewable energy uncertainty on the system, obtain typical renewable energy scenarios by scenario analysis method, and carry out the optimal capacity configuration of the system. Considering the uncertainty of demand response, Ref. [23] used the fuzzy method to deal with the uncertainty. Zhang et al. [24] deal with the uncertainty of electrical load and heating load by a two-stage robust optimization model. Most of the above references only consider single uncertain factors and less consider the combined action of multiple uncertain factors. When using random method and fuzzy method to deal with uncertain problems, it is necessary to obtain accurate probability density distribution and membership relationship, which is difficult to realize. Robust method requires that constraints be satisfied in extreme cases, and the optimization results are relatively conservative.In view of the above problems, this paper puts forward the joint optimization method of planning and operation for integrated electrothermal and hydrogen coupling system (IEHS), as shown in Figure1. Firstly, the frame structure of IEHS is constructed, and multiple utilization modes of hydrogen energy are analyzed. In order to accurately describe the operation characteristics of IEHS, Part II analyzes the physical characteristics of hydrogen storage device and establishes the operation model of high-pressure hydrogen storage device (HHS) based on the van Deemter equation and the linear model of compressor power consumption considering dynamic compression ratio. Secondly, a joint optimization model of planning and operation for IEHS is established in Part III. The upper model carries out the equipment capacity configuration with the goal of minimizing the system construction and operation cost. The lower model optimizes the output of each equipment with the goal of minimizing the system operation cost, considering the multiple uncertain factors. On this basis, Part IV analyzes the influence of combined action of multiple uncertain factors on IEHS, and interval linear programming method is proposed to deal with the uncertainty of renewable energy, load and energy price. Further, the upper model is solved by the multiverse optimizer (MVO) algorithm, and the CPLEX solver is called to solve the bilevel optimization model. Finally, through simulation analysis, it is verified that the joint optimization method of planning and operation for IEHS has advantages in feasibility, robustness, economy, and environmental protection. Further, the multiple utilization modes of hydrogen energy can reduce the operation cost of IEHS, and accurately describing the physical characteristics of the hydrogen storage device is necessity to keep HSD running safely.Figure 1
Logical structure of this thesis.


## 2. Framework Model of Integrated Electrothermal and Hydrogen Coupling System

The structure of integrated electrothermal and hydrogen coupling system (IEHS) studied in this paper is shown in Figure2. The equipment in IEHS mainly consists of wind turbine (WT), photovoltaic (PV), power to hydrogen (P2H), gas-fired boiler (GB), combined heating and power unit (CHP), hydrogen fuel cell (HFC), heating energy storage device (HES), and hydrogen storage device (HSD). Energy input includes power purchase from the utility grid (UG) and nature gas purchase from the gas network. Energy output mainly includes electricity, thermal energy, and hydrogen energy. The terminal load includes electrical load (Eload), heating load (Hload), and hydrogen load (H2load).Figure 2
Frame structure of integrated electrothermal and hydrogen coupling system.In the frame structure of IEHS shown in Figure2, the hydrogen energy system includes the following links. Renewable energy power generation and P2H constitute the hydrogen production link. Referring to hydrogen energy storage technology in Ref. [6], hydrogen storage link consists of low-pressure hydrogen storage device (LHS), high-pressure hydrogen storage device (HHS), and the pressure conversion device, which is constituted by the compressor (CP) of HHS, CP of LHS, and reducing valve (RV). Hydrogen station, hydrogen fuel cell, and hydrogen injection in nature gas network (HING) constitute the hydrogen utilization link, including the following four utilization modes. In mode 1, hydrogen is consumed by hydrogen fuel cell to realize the heat and power generation, i.e., the HFC mode. In mode 2, hydrogen is injected into the natural gas network to form a mixed gas to supply the terminal gas appliances. In mode 3, hydrogen is pressurized by CP and injected into HHS for storage and then pressurized again by CP and supplied to the hydrogen station to meet the hydrogen demand of hydrogen fuel-cell vehicles. In mode 4, hydrogen is pressurized and injected into LHS for large-scale seasonal storage of hydrogen energy. The hydrogen energy supply and demand in modes 1~3 show typical seasonal characteristics, and the seasonal hydrogen storage in mode 4 can realize large-scale energy spatial and temporal transferring to balance supply and demand in modes 1~3. Multiple utilization modes of hydrogen energy can effectively improve the multienergy complementary and coordination of IEHS and ensure the balance between supply and demand of IEHS.In this section, considering the hydrogen physical characteristics, the model of hydrogen storage device is established based on the van Deemter equation. Considering the dynamic compression ratio, the model of compressor power consumption is constructed. In addition to the hydrogen storage link, IEHS also includes various hydrogen production equipment and hydrogen consumption equipment, which are detailed as follows.


### 2.1. Hydrogen Energy System

#### 2.1.1. Hydrogen Production Link

Alkaline electrolytic water technology is the most mature and widely used technology of power to hydrogen (P2H) at present. Alkaline electrolytic water equipment with the rated power of MW levels has been applied to practical projects. The pressure of hydrogen production by alkaline electrolytic water is about 0.1~1.6 MPa, and its operation model is shown in the following equations.(1)mk,tP2H=ηP2HPk,tP2H/HHVH,Bk,tP2HPminP2H≤Pk,tP2H≤Bk,tP2HPmaxP2H.In the formula,mk,tP2H is the amount of hydrogen production in season k at time t. This paper mainly considers three typical seasons, such as summer, transitional season, winter. k∈sum,tra,win. ηP2H is the efficiency of P2H. Pk,tP2H is the power consumption of P2H in season k at time t. HHVH is the higher heating value of hydrogen. PmaxP2H and PminP2H are the maximum and minimum power of P2H, respectively. Bk,tP2H is the state variable of P2H in season k at time t. When Bk,tP2H=1, the P2H is active. When Bk,tP2H=0, the P2H is at shutdown.In order to avoid irreversible damage to the electrode caused by frequent startup and shutdown of P2H, the constrain of startup and shutdown frequency of P2H is shown in the following equation.(2)365×∑k∈sum,tra,winNk∑t=1NtIk,tP2H,on+Ik,tP2H,off/2≤NP2H,cycle,(3)0≤Ik,tP2H,on+Ik,tP2H,off≤1.(4)Ik,tP2H,on−Ik,tP2H,off=Bk,tP2H−Bk,t−1P2H.In the formula,Nk is the proportion of season k in the whole year. Nt is the total number of hours in a day. Ik,tP2H,on and Ik,tP2H,off are, respectively, the start and stop action variables of P2H in season k at time t. If Ik,tP2H,on=1, it means that P2H starts at time t. If Ik,tP2H,off=1, P2H is shut down at time t. NP2H,cycle is the maximum frequency of starts and stops of P2H equipment in a year.


#### 2.1.2. Hydrogen Storage Link

(1) Hydrogen Storage Device. According to the pressure of hydrogen storage device, which can be divided into high-pressure hydrogen storage device (HHS) and low-pressure hydrogen storage device (LHS), the gas pressure of HHS can reach 20 MPa, while the gas pressure of LHS is generally 2~5 MPa, such as in Refs. [25, 26]. When the hydrogen pressure is less than 10 MPa, the perfect gas equation can approximately simulate the relationship between hydrogen pressure and mass, such as in Ref. [19]. The model of LHS based on the perfect gas equation is shown in the following equations.
(5)ptLHSVLHS=ntLHSRTH,(6)mtLHS=MHntLHS.In the formula,ptLHS is the gas pressure of LHS at time t. ntLHS is the amount of substance of hydrogen in LHS at time t. mtLHS is the mass of hydrogen in LHS at time t. R is the ideal gas constant. TH is the gas temperature. MH is the relative molecular weight of hydrogen. VLHS is the LHS volume. According to Equations (5) and (6), when temperature was constant, ptLHS was linear to mtLHS in LHS.Hydrogen, as the material with the lowest relative molecular mass, is more likely to expand under high pressure. The perfect gas equation cannot accurately describe the gravitation and repulsion between high-pressure hydrogen molecules, and it is difficult to accurately reflect the relationship between hydrogen mass and pressure in HHS, such as in Ref. [27]. The van Deemter equation can truly reflect the macroscopic physical properties of gas, such as in Ref. [28]. In order to accurately characterize the hydrogen physical characteristics in HHS, the model of HHS based on the van Deemter equation is shown in the following equation.
(7)ptHHS+antHHS/VHHS2VHHS−bntHHS=ntHHSRTH.In the formula,ptHHS is the gas pressure of HHS at time t. VHHS is the HHS volume. ntHHS is the amount of substance of hydrogen in HHS at time t. a and b are the correction of the gravitation and repulsion between high-pressure hydrogen molecules, respectively.According to Equation (7), the relationship between hydrogen pressure ntHHS and mass mtHHS in HHS is shown in Equation (3).
(8)ptHHS=fmtHHS=mtHHSRTH/MHVHHS−bmtHHS−amtHHS/MHVHHS.According to Equation (8), ptHHS and mtHHS in HHS show a nonlinear behavior. In order to simplify the calculation, the pressure piecewise linear equation is shown in Equation (9) based on the incremental piecewise linearization method, such as in Ref. [29].
(9)ptHHS≈fmp,1+∑i=1Np−1fmp,i+1−fmp,iθi,0≤θi≤1,(10)mtHHS=mp,1+∑i=1Np−1mp,i+1−mp,iθi,(11)θi+1≤ηi≤θi,ηi∈0,1.In the formula,Np is the total number of linearized segments. mp,i is the segment point. θi reflects the position of mtHHS in the i-th segment. ηi is a binary variable, and Equation (11) is used to ensure the continuity of subsection interval.During operation, LHS and HHS shall meet the time sequence equation of residual capacity in hydrogen storage device as Equation (12), safe operation constraint as Equation (13), hydrogen storage and hydrogen release constraint as Equation (14) and the constrain of hydrogen storage state and hydrogen release state as Equation (15).
(12)mt+1HS=mtHS−mtHS.in−mtHS,out,(13)pminHS≤ptHS≤pmaxHS,(14)BtHS,inmmaxHS,in≤mtHS,in≤00≤mtHS,out≤BtHS,outmmaxHS,out,(15)0≤BtHS,out+BtHS,in≤1.In the formula,mtHS is the residual capacity of hydrogen storage device at time t, HS∈LHS,HHS. ptHS is the gas pressure of hydrogen storage device at time t. mtHS,in and mtHS,out are the hydrogen storage rate and hydrogen release rate of hydrogen storage device at time t. mmaxHS,in and mmaxHS,out are the maximum hydrogen storage rate and hydrogen release rate. pminHS and pmaxHS are the minimum and maximum gas pressure of hydrogen storage device. BtHS,in and BtHS,out are the state variables of hydrogen storage and hydrogen release of hydrogen storage device at time t. When BtHS,in=1, the hydrogen storage device is in hydrogen storage state. When BtHS,out=1, the hydrogen storage device is in hydrogen release state.(2) Compressor. Accurate description of compressor power consumption is the key to calculate the energy efficiency of hydrogen storage system. Under the condition of isentropy, the steady-state power consumption of compressor is mainly related to compression ratio and air intake. Existing researches usually set the compression ratio to a fixed value, such as in Ref. [30]. In fact, the compression ratio is related to the inlet and outlet pressures of the compressor. Considering the dynamic change of compression ratio, the model of compressor power consumption is shown in the following equation.
(16)PtCP=WtCPmtin/3.6×106×ηcompMH,(17)WtCP≤RTinγpout/pinγ/γ−1−1/γ−1.In the formula,PtCP is the power consumption of compressor at time t. WtCP is the power consumption of the compressor to compress the unit molar mass of hydrogen at time t. mtin is the gas flow rate of the compressor at time t. γ is the heat capacity ratio of compressed gas. ηcomp is the efficiency of the compressor. pin and pout are the inlet and outlet pressures of the compressor, respectively, and pout/pin is the compression ratio.According to Equations (16) and (17), the compression ratio changes dynamically with the change of inlet and outlet pressures, which affects the compressor power consumption. Considering the dynamic compression ratio, the model of compressor power consumption contains the nonlinear term mWx, which is the product of two continuous variables. In order to simplify the calculation, the binary method and M method in Ref. [13] are used to linearize the nonlinear term as follows.
(18)Wx=Wxmin+y∑n=0p2nδnx,(19)0≤y∑n=0p2nδnx≤Wxmax−Wxmin,δnx∈0,1,(20)mWx=mWxmin+y∑n=0p2nτnx,(21)m−M1−δnx≤τnx≤m+M1−δnx,(22)−Mδnx≤τnx≤Mδnx.In the formula,m and Wx are continuous variables. δnx is a binary variable. τnx is an auxiliary variable. Equation (18) discretizes the continuous variable Wx, which is expressed by δnx. y is the resolution. Equations (20) and (22) linearize the nonlinear term mWx based on the M method.


#### 2.1.3. Hydrogen Utilization Link

(1) Hydrogen Fuel Cell. The gas pressure of hydrogen supplied to hydrogen fuel cell (HFC) is usually 0.03~0.06 MPa in Ref. [31], and the thermal energy can be recovered through the heat collection device to realize the cogeneration of HFC in Ref. [26], as shown in Equations (23) and (24).
(23)Pk,tHFC,e=ηeHFCmk,tHFCHHVH,(24)Pk,tHFC,h=ηhHFC1−ηeHFCmk,tHFCHHVH,(25)0≤Pk,tHFC,e≤Bk,tHFCPmaxHFC,e.In the formula,mk,tHFC is the hydrogen consumption of HFC in season k at time t. Pk,tHFC,e and Pk,tHFC,h are the electric power and thermal power of HFC in season k at time t, respectively. ηeHFC and ηhHFC are the power generation efficiency and heat generation efficiency of HFC. Bk,tHFC is the state variable of HFC in season k at time t. When Bk,tHFC=1, the HFC is active. When Bk,tHFC=0, the HFC is at shutdown.(2) Model of Hydrogen Injection in Natural Gas Network. Hydrogen injection in natural gas network can transmit hydrogen without increasing additional investment, which is one of the effective methods to realize large-scale transmission and supply of hydrogen, such as in Ref. [13]. Assuming that hydrogen is evenly distributed in the natural gas network, the proportion of hydrogen in the natural gas network is shown in the following equation.
(26)ωkH=Gk,tH,inject/Gk,tmix,Gk,tmix=Gk,tH,inject+Gk,tNG.In the formula,ωkH is the proportion of hydrogen injection in the nature gas network in season k. Due to the phenomenon of hydrogen embrittlement caused by adding more hydrogen, ωkH is generally not higher than 20%, such as in Ref. [13]. Gk,tmix is the flow of hydrogen, nature gas mixed gas in the natural gas network in season k at time t. Gk,tH,inject is the hydrogen injection flow in season k at time t. Gk,tNG is the natural gas flow in season k at time t.For the calorific value of hydrogen, nature gas mixed gas is related to the hydrogen injection ratio, as shown in the following equation.(27)HHVmix=ωkHHHVH+1−ωkHHHVNG.In formula,HHVmix is the higher heating value of hydrogen, nature gas mixed gas.


### 2.2. Electrothermal Coupling System

#### 2.2.1. Combined Heat and Power Generation Unit

Combined heat and power generation unit (CHP) and gas-fired boiler (GB) are important heat sources and power sources in IEHS. They have little impact on equipment operation under a certain hydrogen injection ratio. Assuming that the equipment operation efficiency is constant in Ref. [32], the operation model of CHP is shown in the following equations.
(28)Pk,tCHP,e=ηeCHPHHVmixGk,tCHP,Pk,tCHP,h=ηhCHPHHVmixGk,tCHP,0≤Pk,tCHP,e≤PmaxCHP,e.In the formula,Pk,tCHP,e and Pk,tCHP,h are the electric power and thermal power of CHP in season k at time t. ηeCHP and ηhCHP are the power generation efficiency and heat generation efficiency of CHP. Gk,tCHP is the gas consumption of CHP in season k at time t. PmaxCHP,e is the maximum electric power of CHP.


#### 2.2.2. Gas-Fired Boiler

(29)Pk,tGB,h=ηhGBHHVmixGk,tGB,0≤Pk,tGB,h≤PmaxGB,h.In the formula,Pk,tGB,h is the thermal power of GB in season k at time t. ηhGB is the heat generation efficiency of GB. Gk,tGB is the gas consumption of GB in season k at time t. PmaxGB,h is the maximum thermal power of GB.


#### 2.2.3. Heating Energy Storage Device

The operation model of heating energy storage device (HES) is similar to that of hydrogen storage device, as shown below.(30)St+1HES=StHES−PtHES,out/ηoutHES−PtHES,inηinHES,0.2×SinstHES≤StHES≤SinstHES,0.3×SinstHESBtHES,in≤PtHES,in≤00≤PtHES,out≤0.3×SinstHESBtHES,out,0≤BtHES,out+BtHES,in≤1.In the formula,StHES is the residual capacity of HES at time t. PtHES,in and PtHES,out refer to the heat storage and heat release of HES at time t. ηinHES and ηoutHES are the heat storage rate and heat release rate of HES. SinstHES is the installed capacity of HES. BtHES,in and BtHES,out are, respectively, the state of heat storage and heat release of HES at time t. When BtHES,in=1, HES is in the heat storage state, and when BtHES,out=1, HES is in the heat release state.


## 2.1. Hydrogen Energy System

### 2.1.1. Hydrogen Production Link

Alkaline electrolytic water technology is the most mature and widely used technology of power to hydrogen (P2H) at present. Alkaline electrolytic water equipment with the rated power of MW levels has been applied to practical projects. The pressure of hydrogen production by alkaline electrolytic water is about 0.1~1.6 MPa, and its operation model is shown in the following equations.(1)mk,tP2H=ηP2HPk,tP2H/HHVH,Bk,tP2HPminP2H≤Pk,tP2H≤Bk,tP2HPmaxP2H.In the formula,mk,tP2H is the amount of hydrogen production in season k at time t. This paper mainly considers three typical seasons, such as summer, transitional season, winter. k∈sum,tra,win. ηP2H is the efficiency of P2H. Pk,tP2H is the power consumption of P2H in season k at time t. HHVH is the higher heating value of hydrogen. PmaxP2H and PminP2H are the maximum and minimum power of P2H, respectively. Bk,tP2H is the state variable of P2H in season k at time t. When Bk,tP2H=1, the P2H is active. When Bk,tP2H=0, the P2H is at shutdown.In order to avoid irreversible damage to the electrode caused by frequent startup and shutdown of P2H, the constrain of startup and shutdown frequency of P2H is shown in the following equation.(2)365×∑k∈sum,tra,winNk∑t=1NtIk,tP2H,on+Ik,tP2H,off/2≤NP2H,cycle,(3)0≤Ik,tP2H,on+Ik,tP2H,off≤1.(4)Ik,tP2H,on−Ik,tP2H,off=Bk,tP2H−Bk,t−1P2H.In the formula,Nk is the proportion of season k in the whole year. Nt is the total number of hours in a day. Ik,tP2H,on and Ik,tP2H,off are, respectively, the start and stop action variables of P2H in season k at time t. If Ik,tP2H,on=1, it means that P2H starts at time t. If Ik,tP2H,off=1, P2H is shut down at time t. NP2H,cycle is the maximum frequency of starts and stops of P2H equipment in a year.


### 2.1.2. Hydrogen Storage Link

(1) Hydrogen Storage Device. According to the pressure of hydrogen storage device, which can be divided into high-pressure hydrogen storage device (HHS) and low-pressure hydrogen storage device (LHS), the gas pressure of HHS can reach 20 MPa, while the gas pressure of LHS is generally 2~5 MPa, such as in Refs. [25, 26]. When the hydrogen pressure is less than 10 MPa, the perfect gas equation can approximately simulate the relationship between hydrogen pressure and mass, such as in Ref. [19]. The model of LHS based on the perfect gas equation is shown in the following equations.
(5)ptLHSVLHS=ntLHSRTH,(6)mtLHS=MHntLHS.In the formula,ptLHS is the gas pressure of LHS at time t. ntLHS is the amount of substance of hydrogen in LHS at time t. mtLHS is the mass of hydrogen in LHS at time t. R is the ideal gas constant. TH is the gas temperature. MH is the relative molecular weight of hydrogen. VLHS is the LHS volume. According to Equations (5) and (6), when temperature was constant, ptLHS was linear to mtLHS in LHS.Hydrogen, as the material with the lowest relative molecular mass, is more likely to expand under high pressure. The perfect gas equation cannot accurately describe the gravitation and repulsion between high-pressure hydrogen molecules, and it is difficult to accurately reflect the relationship between hydrogen mass and pressure in HHS, such as in Ref. [27]. The van Deemter equation can truly reflect the macroscopic physical properties of gas, such as in Ref. [28]. In order to accurately characterize the hydrogen physical characteristics in HHS, the model of HHS based on the van Deemter equation is shown in the following equation.
(7)ptHHS+antHHS/VHHS2VHHS−bntHHS=ntHHSRTH.In the formula,ptHHS is the gas pressure of HHS at time t. VHHS is the HHS volume. ntHHS is the amount of substance of hydrogen in HHS at time t. a and b are the correction of the gravitation and repulsion between high-pressure hydrogen molecules, respectively.According to Equation (7), the relationship between hydrogen pressure ntHHS and mass mtHHS in HHS is shown in Equation (3).
(8)ptHHS=fmtHHS=mtHHSRTH/MHVHHS−bmtHHS−amtHHS/MHVHHS.According to Equation (8), ptHHS and mtHHS in HHS show a nonlinear behavior. In order to simplify the calculation, the pressure piecewise linear equation is shown in Equation (9) based on the incremental piecewise linearization method, such as in Ref. [29].
(9)ptHHS≈fmp,1+∑i=1Np−1fmp,i+1−fmp,iθi,0≤θi≤1,(10)mtHHS=mp,1+∑i=1Np−1mp,i+1−mp,iθi,(11)θi+1≤ηi≤θi,ηi∈0,1.In the formula,Np is the total number of linearized segments. mp,i is the segment point. θi reflects the position of mtHHS in the i-th segment. ηi is a binary variable, and Equation (11) is used to ensure the continuity of subsection interval.During operation, LHS and HHS shall meet the time sequence equation of residual capacity in hydrogen storage device as Equation (12), safe operation constraint as Equation (13), hydrogen storage and hydrogen release constraint as Equation (14) and the constrain of hydrogen storage state and hydrogen release state as Equation (15).
(12)mt+1HS=mtHS−mtHS.in−mtHS,out,(13)pminHS≤ptHS≤pmaxHS,(14)BtHS,inmmaxHS,in≤mtHS,in≤00≤mtHS,out≤BtHS,outmmaxHS,out,(15)0≤BtHS,out+BtHS,in≤1.In the formula,mtHS is the residual capacity of hydrogen storage device at time t, HS∈LHS,HHS. ptHS is the gas pressure of hydrogen storage device at time t. mtHS,in and mtHS,out are the hydrogen storage rate and hydrogen release rate of hydrogen storage device at time t. mmaxHS,in and mmaxHS,out are the maximum hydrogen storage rate and hydrogen release rate. pminHS and pmaxHS are the minimum and maximum gas pressure of hydrogen storage device. BtHS,in and BtHS,out are the state variables of hydrogen storage and hydrogen release of hydrogen storage device at time t. When BtHS,in=1, the hydrogen storage device is in hydrogen storage state. When BtHS,out=1, the hydrogen storage device is in hydrogen release state.(2) Compressor. Accurate description of compressor power consumption is the key to calculate the energy efficiency of hydrogen storage system. Under the condition of isentropy, the steady-state power consumption of compressor is mainly related to compression ratio and air intake. Existing researches usually set the compression ratio to a fixed value, such as in Ref. [30]. In fact, the compression ratio is related to the inlet and outlet pressures of the compressor. Considering the dynamic change of compression ratio, the model of compressor power consumption is shown in the following equation.
(16)PtCP=WtCPmtin/3.6×106×ηcompMH,(17)WtCP≤RTinγpout/pinγ/γ−1−1/γ−1.In the formula,PtCP is the power consumption of compressor at time t. WtCP is the power consumption of the compressor to compress the unit molar mass of hydrogen at time t. mtin is the gas flow rate of the compressor at time t. γ is the heat capacity ratio of compressed gas. ηcomp is the efficiency of the compressor. pin and pout are the inlet and outlet pressures of the compressor, respectively, and pout/pin is the compression ratio.According to Equations (16) and (17), the compression ratio changes dynamically with the change of inlet and outlet pressures, which affects the compressor power consumption. Considering the dynamic compression ratio, the model of compressor power consumption contains the nonlinear term mWx, which is the product of two continuous variables. In order to simplify the calculation, the binary method and M method in Ref. [13] are used to linearize the nonlinear term as follows.
(18)Wx=Wxmin+y∑n=0p2nδnx,(19)0≤y∑n=0p2nδnx≤Wxmax−Wxmin,δnx∈0,1,(20)mWx=mWxmin+y∑n=0p2nτnx,(21)m−M1−δnx≤τnx≤m+M1−δnx,(22)−Mδnx≤τnx≤Mδnx.In the formula,m and Wx are continuous variables. δnx is a binary variable. τnx is an auxiliary variable. Equation (18) discretizes the continuous variable Wx, which is expressed by δnx. y is the resolution. Equations (20) and (22) linearize the nonlinear term mWx based on the M method.


### 2.1.3. Hydrogen Utilization Link

(1) Hydrogen Fuel Cell. The gas pressure of hydrogen supplied to hydrogen fuel cell (HFC) is usually 0.03~0.06 MPa in Ref. [31], and the thermal energy can be recovered through the heat collection device to realize the cogeneration of HFC in Ref. [26], as shown in Equations (23) and (24).
(23)Pk,tHFC,e=ηeHFCmk,tHFCHHVH,(24)Pk,tHFC,h=ηhHFC1−ηeHFCmk,tHFCHHVH,(25)0≤Pk,tHFC,e≤Bk,tHFCPmaxHFC,e.In the formula,mk,tHFC is the hydrogen consumption of HFC in season k at time t. Pk,tHFC,e and Pk,tHFC,h are the electric power and thermal power of HFC in season k at time t, respectively. ηeHFC and ηhHFC are the power generation efficiency and heat generation efficiency of HFC. Bk,tHFC is the state variable of HFC in season k at time t. When Bk,tHFC=1, the HFC is active. When Bk,tHFC=0, the HFC is at shutdown.(2) Model of Hydrogen Injection in Natural Gas Network. Hydrogen injection in natural gas network can transmit hydrogen without increasing additional investment, which is one of the effective methods to realize large-scale transmission and supply of hydrogen, such as in Ref. [13]. Assuming that hydrogen is evenly distributed in the natural gas network, the proportion of hydrogen in the natural gas network is shown in the following equation.
(26)ωkH=Gk,tH,inject/Gk,tmix,Gk,tmix=Gk,tH,inject+Gk,tNG.In the formula,ωkH is the proportion of hydrogen injection in the nature gas network in season k. Due to the phenomenon of hydrogen embrittlement caused by adding more hydrogen, ωkH is generally not higher than 20%, such as in Ref. [13]. Gk,tmix is the flow of hydrogen, nature gas mixed gas in the natural gas network in season k at time t. Gk,tH,inject is the hydrogen injection flow in season k at time t. Gk,tNG is the natural gas flow in season k at time t.For the calorific value of hydrogen, nature gas mixed gas is related to the hydrogen injection ratio, as shown in the following equation.(27)HHVmix=ωkHHHVH+1−ωkHHHVNG.In formula,HHVmix is the higher heating value of hydrogen, nature gas mixed gas.


## 2.1.1. Hydrogen Production Link

Alkaline electrolytic water technology is the most mature and widely used technology of power to hydrogen (P2H) at present. Alkaline electrolytic water equipment with the rated power of MW levels has been applied to practical projects. The pressure of hydrogen production by alkaline electrolytic water is about 0.1~1.6 MPa, and its operation model is shown in the following equations.(1)mk,tP2H=ηP2HPk,tP2H/HHVH,Bk,tP2HPminP2H≤Pk,tP2H≤Bk,tP2HPmaxP2H.In the formula,mk,tP2H is the amount of hydrogen production in season k at time t. This paper mainly considers three typical seasons, such as summer, transitional season, winter. k∈sum,tra,win. ηP2H is the efficiency of P2H. Pk,tP2H is the power consumption of P2H in season k at time t. HHVH is the higher heating value of hydrogen. PmaxP2H and PminP2H are the maximum and minimum power of P2H, respectively. Bk,tP2H is the state variable of P2H in season k at time t. When Bk,tP2H=1, the P2H is active. When Bk,tP2H=0, the P2H is at shutdown.In order to avoid irreversible damage to the electrode caused by frequent startup and shutdown of P2H, the constrain of startup and shutdown frequency of P2H is shown in the following equation.(2)365×∑k∈sum,tra,winNk∑t=1NtIk,tP2H,on+Ik,tP2H,off/2≤NP2H,cycle,(3)0≤Ik,tP2H,on+Ik,tP2H,off≤1.(4)Ik,tP2H,on−Ik,tP2H,off=Bk,tP2H−Bk,t−1P2H.In the formula,Nk is the proportion of season k in the whole year. Nt is the total number of hours in a day. Ik,tP2H,on and Ik,tP2H,off are, respectively, the start and stop action variables of P2H in season k at time t. If Ik,tP2H,on=1, it means that P2H starts at time t. If Ik,tP2H,off=1, P2H is shut down at time t. NP2H,cycle is the maximum frequency of starts and stops of P2H equipment in a year.


## 2.1.2. Hydrogen Storage Link

(1) Hydrogen Storage Device. According to the pressure of hydrogen storage device, which can be divided into high-pressure hydrogen storage device (HHS) and low-pressure hydrogen storage device (LHS), the gas pressure of HHS can reach 20 MPa, while the gas pressure of LHS is generally 2~5 MPa, such as in Refs. [25, 26]. When the hydrogen pressure is less than 10 MPa, the perfect gas equation can approximately simulate the relationship between hydrogen pressure and mass, such as in Ref. [19]. The model of LHS based on the perfect gas equation is shown in the following equations.
(5)ptLHSVLHS=ntLHSRTH,(6)mtLHS=MHntLHS.In the formula,ptLHS is the gas pressure of LHS at time t. ntLHS is the amount of substance of hydrogen in LHS at time t. mtLHS is the mass of hydrogen in LHS at time t. R is the ideal gas constant. TH is the gas temperature. MH is the relative molecular weight of hydrogen. VLHS is the LHS volume. According to Equations (5) and (6), when temperature was constant, ptLHS was linear to mtLHS in LHS.Hydrogen, as the material with the lowest relative molecular mass, is more likely to expand under high pressure. The perfect gas equation cannot accurately describe the gravitation and repulsion between high-pressure hydrogen molecules, and it is difficult to accurately reflect the relationship between hydrogen mass and pressure in HHS, such as in Ref. [27]. The van Deemter equation can truly reflect the macroscopic physical properties of gas, such as in Ref. [28]. In order to accurately characterize the hydrogen physical characteristics in HHS, the model of HHS based on the van Deemter equation is shown in the following equation.
(7)ptHHS+antHHS/VHHS2VHHS−bntHHS=ntHHSRTH.In the formula,ptHHS is the gas pressure of HHS at time t. VHHS is the HHS volume. ntHHS is the amount of substance of hydrogen in HHS at time t. a and b are the correction of the gravitation and repulsion between high-pressure hydrogen molecules, respectively.According to Equation (7), the relationship between hydrogen pressure ntHHS and mass mtHHS in HHS is shown in Equation (3).
(8)ptHHS=fmtHHS=mtHHSRTH/MHVHHS−bmtHHS−amtHHS/MHVHHS.According to Equation (8), ptHHS and mtHHS in HHS show a nonlinear behavior. In order to simplify the calculation, the pressure piecewise linear equation is shown in Equation (9) based on the incremental piecewise linearization method, such as in Ref. [29].
(9)ptHHS≈fmp,1+∑i=1Np−1fmp,i+1−fmp,iθi,0≤θi≤1,(10)mtHHS=mp,1+∑i=1Np−1mp,i+1−mp,iθi,(11)θi+1≤ηi≤θi,ηi∈0,1.In the formula,Np is the total number of linearized segments. mp,i is the segment point. θi reflects the position of mtHHS in the i-th segment. ηi is a binary variable, and Equation (11) is used to ensure the continuity of subsection interval.During operation, LHS and HHS shall meet the time sequence equation of residual capacity in hydrogen storage device as Equation (12), safe operation constraint as Equation (13), hydrogen storage and hydrogen release constraint as Equation (14) and the constrain of hydrogen storage state and hydrogen release state as Equation (15).
(12)mt+1HS=mtHS−mtHS.in−mtHS,out,(13)pminHS≤ptHS≤pmaxHS,(14)BtHS,inmmaxHS,in≤mtHS,in≤00≤mtHS,out≤BtHS,outmmaxHS,out,(15)0≤BtHS,out+BtHS,in≤1.In the formula,mtHS is the residual capacity of hydrogen storage device at time t, HS∈LHS,HHS. ptHS is the gas pressure of hydrogen storage device at time t. mtHS,in and mtHS,out are the hydrogen storage rate and hydrogen release rate of hydrogen storage device at time t. mmaxHS,in and mmaxHS,out are the maximum hydrogen storage rate and hydrogen release rate. pminHS and pmaxHS are the minimum and maximum gas pressure of hydrogen storage device. BtHS,in and BtHS,out are the state variables of hydrogen storage and hydrogen release of hydrogen storage device at time t. When BtHS,in=1, the hydrogen storage device is in hydrogen storage state. When BtHS,out=1, the hydrogen storage device is in hydrogen release state.(2) Compressor. Accurate description of compressor power consumption is the key to calculate the energy efficiency of hydrogen storage system. Under the condition of isentropy, the steady-state power consumption of compressor is mainly related to compression ratio and air intake. Existing researches usually set the compression ratio to a fixed value, such as in Ref. [30]. In fact, the compression ratio is related to the inlet and outlet pressures of the compressor. Considering the dynamic change of compression ratio, the model of compressor power consumption is shown in the following equation.
(16)PtCP=WtCPmtin/3.6×106×ηcompMH,(17)WtCP≤RTinγpout/pinγ/γ−1−1/γ−1.In the formula,PtCP is the power consumption of compressor at time t. WtCP is the power consumption of the compressor to compress the unit molar mass of hydrogen at time t. mtin is the gas flow rate of the compressor at time t. γ is the heat capacity ratio of compressed gas. ηcomp is the efficiency of the compressor. pin and pout are the inlet and outlet pressures of the compressor, respectively, and pout/pin is the compression ratio.According to Equations (16) and (17), the compression ratio changes dynamically with the change of inlet and outlet pressures, which affects the compressor power consumption. Considering the dynamic compression ratio, the model of compressor power consumption contains the nonlinear term mWx, which is the product of two continuous variables. In order to simplify the calculation, the binary method and M method in Ref. [13] are used to linearize the nonlinear term as follows.
(18)Wx=Wxmin+y∑n=0p2nδnx,(19)0≤y∑n=0p2nδnx≤Wxmax−Wxmin,δnx∈0,1,(20)mWx=mWxmin+y∑n=0p2nτnx,(21)m−M1−δnx≤τnx≤m+M1−δnx,(22)−Mδnx≤τnx≤Mδnx.In the formula,m and Wx are continuous variables. δnx is a binary variable. τnx is an auxiliary variable. Equation (18) discretizes the continuous variable Wx, which is expressed by δnx. y is the resolution. Equations (20) and (22) linearize the nonlinear term mWx based on the M method.


## 2.1.3. Hydrogen Utilization Link

(1) Hydrogen Fuel Cell. The gas pressure of hydrogen supplied to hydrogen fuel cell (HFC) is usually 0.03~0.06 MPa in Ref. [31], and the thermal energy can be recovered through the heat collection device to realize the cogeneration of HFC in Ref. [26], as shown in Equations (23) and (24).
(23)Pk,tHFC,e=ηeHFCmk,tHFCHHVH,(24)Pk,tHFC,h=ηhHFC1−ηeHFCmk,tHFCHHVH,(25)0≤Pk,tHFC,e≤Bk,tHFCPmaxHFC,e.In the formula,mk,tHFC is the hydrogen consumption of HFC in season k at time t. Pk,tHFC,e and Pk,tHFC,h are the electric power and thermal power of HFC in season k at time t, respectively. ηeHFC and ηhHFC are the power generation efficiency and heat generation efficiency of HFC. Bk,tHFC is the state variable of HFC in season k at time t. When Bk,tHFC=1, the HFC is active. When Bk,tHFC=0, the HFC is at shutdown.(2) Model of Hydrogen Injection in Natural Gas Network. Hydrogen injection in natural gas network can transmit hydrogen without increasing additional investment, which is one of the effective methods to realize large-scale transmission and supply of hydrogen, such as in Ref. [13]. Assuming that hydrogen is evenly distributed in the natural gas network, the proportion of hydrogen in the natural gas network is shown in the following equation.
(26)ωkH=Gk,tH,inject/Gk,tmix,Gk,tmix=Gk,tH,inject+Gk,tNG.In the formula,ωkH is the proportion of hydrogen injection in the nature gas network in season k. Due to the phenomenon of hydrogen embrittlement caused by adding more hydrogen, ωkH is generally not higher than 20%, such as in Ref. [13]. Gk,tmix is the flow of hydrogen, nature gas mixed gas in the natural gas network in season k at time t. Gk,tH,inject is the hydrogen injection flow in season k at time t. Gk,tNG is the natural gas flow in season k at time t.For the calorific value of hydrogen, nature gas mixed gas is related to the hydrogen injection ratio, as shown in the following equation.(27)HHVmix=ωkHHHVH+1−ωkHHHVNG.In formula,HHVmix is the higher heating value of hydrogen, nature gas mixed gas.


## 2.2. Electrothermal Coupling System

### 2.2.1. Combined Heat and Power Generation Unit

Combined heat and power generation unit (CHP) and gas-fired boiler (GB) are important heat sources and power sources in IEHS. They have little impact on equipment operation under a certain hydrogen injection ratio. Assuming that the equipment operation efficiency is constant in Ref. [32], the operation model of CHP is shown in the following equations.
(28)Pk,tCHP,e=ηeCHPHHVmixGk,tCHP,Pk,tCHP,h=ηhCHPHHVmixGk,tCHP,0≤Pk,tCHP,e≤PmaxCHP,e.In the formula,Pk,tCHP,e and Pk,tCHP,h are the electric power and thermal power of CHP in season k at time t. ηeCHP and ηhCHP are the power generation efficiency and heat generation efficiency of CHP. Gk,tCHP is the gas consumption of CHP in season k at time t. PmaxCHP,e is the maximum electric power of CHP.


### 2.2.2. Gas-Fired Boiler

(29)Pk,tGB,h=ηhGBHHVmixGk,tGB,0≤Pk,tGB,h≤PmaxGB,h.In the formula,Pk,tGB,h is the thermal power of GB in season k at time t. ηhGB is the heat generation efficiency of GB. Gk,tGB is the gas consumption of GB in season k at time t. PmaxGB,h is the maximum thermal power of GB.


### 2.2.3. Heating Energy Storage Device

The operation model of heating energy storage device (HES) is similar to that of hydrogen storage device, as shown below.(30)St+1HES=StHES−PtHES,out/ηoutHES−PtHES,inηinHES,0.2×SinstHES≤StHES≤SinstHES,0.3×SinstHESBtHES,in≤PtHES,in≤00≤PtHES,out≤0.3×SinstHESBtHES,out,0≤BtHES,out+BtHES,in≤1.In the formula,StHES is the residual capacity of HES at time t. PtHES,in and PtHES,out refer to the heat storage and heat release of HES at time t. ηinHES and ηoutHES are the heat storage rate and heat release rate of HES. SinstHES is the installed capacity of HES. BtHES,in and BtHES,out are, respectively, the state of heat storage and heat release of HES at time t. When BtHES,in=1, HES is in the heat storage state, and when BtHES,out=1, HES is in the heat release state.


## 2.2.1. Combined Heat and Power Generation Unit

Combined heat and power generation unit (CHP) and gas-fired boiler (GB) are important heat sources and power sources in IEHS. They have little impact on equipment operation under a certain hydrogen injection ratio. Assuming that the equipment operation efficiency is constant in Ref. [32], the operation model of CHP is shown in the following equations.
(28)Pk,tCHP,e=ηeCHPHHVmixGk,tCHP,Pk,tCHP,h=ηhCHPHHVmixGk,tCHP,0≤Pk,tCHP,e≤PmaxCHP,e.In the formula,Pk,tCHP,e and Pk,tCHP,h are the electric power and thermal power of CHP in season k at time t. ηeCHP and ηhCHP are the power generation efficiency and heat generation efficiency of CHP. Gk,tCHP is the gas consumption of CHP in season k at time t. PmaxCHP,e is the maximum electric power of CHP.


## 2.2.2. Gas-Fired Boiler

(29)Pk,tGB,h=ηhGBHHVmixGk,tGB,0≤Pk,tGB,h≤PmaxGB,h.In the formula,Pk,tGB,h is the thermal power of GB in season k at time t. ηhGB is the heat generation efficiency of GB. Gk,tGB is the gas consumption of GB in season k at time t. PmaxGB,h is the maximum thermal power of GB.


## 2.2.3. Heating Energy Storage Device

The operation model of heating energy storage device (HES) is similar to that of hydrogen storage device, as shown below.(30)St+1HES=StHES−PtHES,out/ηoutHES−PtHES,inηinHES,0.2×SinstHES≤StHES≤SinstHES,0.3×SinstHESBtHES,in≤PtHES,in≤00≤PtHES,out≤0.3×SinstHESBtHES,out,0≤BtHES,out+BtHES,in≤1.In the formula,StHES is the residual capacity of HES at time t. PtHES,in and PtHES,out refer to the heat storage and heat release of HES at time t. ηinHES and ηoutHES are the heat storage rate and heat release rate of HES. SinstHES is the installed capacity of HES. BtHES,in and BtHES,out are, respectively, the state of heat storage and heat release of HES at time t. When BtHES,in=1, HES is in the heat storage state, and when BtHES,out=1, HES is in the heat release state.


## 3. Bilevel Optimization Model of Integrated Electrothermal and Hydrogen Coupling System

According to the frame structure of IEHS shown in Figure2, a bilevel optimization model is established in Part 3, considering multiple utilization modes of hydrogen energy and multiuncertainty, such as renewable energy, load, and energy price. The upper model is aimed at minimizing the total cost such as construction cost and operation cost within the service life of IEHS, and the lower model is aimed at minimizing the operation cost of IEHS. Firstly, the upper model randomly generates the configuration scheme of each equipment. According to each equipment capacity generated by the upper model, the lower model optimizes each equipment output, which is transmitted to the upper model. And then, the upper model modifies the configuration capacity of each equipment. Through the loop iteration of the upper and lower layers, the optimal planning and operation scheme of IEHS can be obtained.


### 3.1. Upper Model

The upper model is aimed at minimizing the total cost such as construction cost and operation cost within the service life of IEHS, as shown in Equation (31). Among them, the decision variable is the configuration capacity of each equipment.
(31)minCtot=Cope+Cinv.In the formula,Ctot is the total cost within the service life of IEHS. Cope is the operation cost of IEHS. Cinv is the total investment and construction cost of equipment in IEHS.
(32)Cinv=CCHPSinstCHP+CGBSinstGB+CP2HSinstP2H+CHFCSinstHFC+CHESSinstHES+CLHSSinstLHS+CHHSSinstHHS.In the formula,CCHP, CGB, CP2H, CHFC, CHES, CLHS, and CHHS are the unit construction cost of CHP, GB, P2H, HFC, HES, and hydrogen storage device, respectively. SinstCHP, SinstGB, SinstP2H, SinstHFC, SinstHES, SinstLHS, and SinstHHS are the installed capacity of CHP, GB, P2H, HFC, HES, and hydrogen storage device, respectively.


### 3.2. Lower Model

The lower model is aimed at minimizing the annual operation cost of the IEHS, including annual energy purchase costCbuy, carbon emission cost Cem, and operation and maintenance cost Com, as shown in Equation (33). Among them, the decision variable is each equipment output.
(33)minCope=365∑m=1N0m1+rm×∑k∈sum,tra,winNk∑t=1NtCbuy+Cem+Com,(34)Cbuy=Pk,tgridck,tgrid+Gk,tNGck,tNG+mk,tH,bck,tH,b,(35)Cem=cemτgridPk,tgrid+τNGHHVNGGk,tNG+τHmk,tH,b,(36)Com=comCHPPk,tCHP,e+comHFCPk,tHFC,e+comGBPk,tGB,h+comP2HPk,tP2H+comWTPk,tWT+comPVPk,tPV+comHESPk,gHES,out−Pk,gHES,in+comLHSHHVHmk,tLHS,out−mk,tLHS,in+comHHSHHVHmk,tHHS,out−mk,tHHS,in.In the formula,N0 is the planning period. r is the discount rate. Nk is the proportion of season kin the whole year. This paper mainly considers three typical seasons, such as summer, transitional season, and winter. k∈sum,tra,win. Nt is the total number of hours in a day. ck,tgrid, ck,tNG, and ck,tH,b are the unit prices of electricity purchased from the power grid, nature gas purchased from the gas network, and hydrogen purchased from the hydrogen energy market in season k at time t. Pk,tgrid is the exchange power between IEHS and power grid in season k at time t. When Pk,tgrid>0, IEHS purchases power from power grid. When Pk,tgrid<0, IEHS sells power to power grid. Gk,tNG is the nature gas volume purchased by IEHS from the gas network in season k at time t. mk,tH,b is the volume of hydrogen purchased by IEHS from the hydrogen energy market in season k at time t. cem is the unit price of carbon emission. τgrid, τNG, and τH is the unit carbon emissions of electricity purchase, nature gas purchase, and hydrogen purchase, respectively. comCHP, comHFC, comGB, comP2H, comWT, comPV, comHES, comLHS, and comHHS are unit operation and maintenance costs of CHP, HFC, GB, P2H, WT, PV, HES, and hydrogen storage device. Pk,tCHP,e, Pk,tHFC,e, Pk,tGB,h, Pk,tP2H, Pk,tWT, and Pk,tPV are the output power of CHP, HFC, GB, P2H, WT, and PV, respectively. For the absolute value of nonlinear term in Equation (34), auxiliary variable P¯k,tgrid is introduced to linearize it, and constraint such as P¯k,tgrid≥±Pk,tgrid is added.


### 3.3. Constraint Condition

The bilevel optimization of IEHS shall meet the energy balance constraints of electricity, gas, heat, and hydrogen, as detailed below.(1)Electrical Power Balance ConstraintThe power system model in this paper adopts the linear power flow model in Ref. [33], and the nodal power balance constraint is shown in the following equation.
(37)Pk,tL,e=Pk,tgrid+Pk,tWT+Pk,tPV+Pk,tCHP,e+Pk,tHFC,e−Pk,tP2H.In the formula,Pk,tL,e is the electrical load in season k at time t.(2)Natural Gas Flow Balance Constraint(38)Gk,tmix=Gk,tCHP+Gk,tGBGminNG≤Gk,tmix≤GmaxNG.In the formula,GmaxNG and GminNG are the upper and lower limit of natural gas pipeline flow rate.(3)Thermal Energy Balance Constraint(39)Pk,tL,h=Pk,tCHP,h+Pk,tHFC,h+Pk,tGB,h.In the formula,Pk,tL,h is the heating load in season k at time t.(4)Hydrogen Energy Balance Constraint(40)mk,tL,H+mk,tHFC+mk,tH,inject−mk,tHHS,in−mk,tLHS,in=mk,tP2H+mk,tH,b+mk,tHHS,out+mk,tLHS,out.In the formula,mk,tL,H is the load of hydrogen station in season k at time t. mk,tH,inject is the amount of hydrogen injection in nature gas network in season k at time t.


## 3.1. Upper Model

The upper model is aimed at minimizing the total cost such as construction cost and operation cost within the service life of IEHS, as shown in Equation (31). Among them, the decision variable is the configuration capacity of each equipment.
(31)minCtot=Cope+Cinv.In the formula,Ctot is the total cost within the service life of IEHS. Cope is the operation cost of IEHS. Cinv is the total investment and construction cost of equipment in IEHS.
(32)Cinv=CCHPSinstCHP+CGBSinstGB+CP2HSinstP2H+CHFCSinstHFC+CHESSinstHES+CLHSSinstLHS+CHHSSinstHHS.In the formula,CCHP, CGB, CP2H, CHFC, CHES, CLHS, and CHHS are the unit construction cost of CHP, GB, P2H, HFC, HES, and hydrogen storage device, respectively. SinstCHP, SinstGB, SinstP2H, SinstHFC, SinstHES, SinstLHS, and SinstHHS are the installed capacity of CHP, GB, P2H, HFC, HES, and hydrogen storage device, respectively.


## 3.2. Lower Model

The lower model is aimed at minimizing the annual operation cost of the IEHS, including annual energy purchase costCbuy, carbon emission cost Cem, and operation and maintenance cost Com, as shown in Equation (33). Among them, the decision variable is each equipment output.
(33)minCope=365∑m=1N0m1+rm×∑k∈sum,tra,winNk∑t=1NtCbuy+Cem+Com,(34)Cbuy=Pk,tgridck,tgrid+Gk,tNGck,tNG+mk,tH,bck,tH,b,(35)Cem=cemτgridPk,tgrid+τNGHHVNGGk,tNG+τHmk,tH,b,(36)Com=comCHPPk,tCHP,e+comHFCPk,tHFC,e+comGBPk,tGB,h+comP2HPk,tP2H+comWTPk,tWT+comPVPk,tPV+comHESPk,gHES,out−Pk,gHES,in+comLHSHHVHmk,tLHS,out−mk,tLHS,in+comHHSHHVHmk,tHHS,out−mk,tHHS,in.In the formula,N0 is the planning period. r is the discount rate. Nk is the proportion of season kin the whole year. This paper mainly considers three typical seasons, such as summer, transitional season, and winter. k∈sum,tra,win. Nt is the total number of hours in a day. ck,tgrid, ck,tNG, and ck,tH,b are the unit prices of electricity purchased from the power grid, nature gas purchased from the gas network, and hydrogen purchased from the hydrogen energy market in season k at time t. Pk,tgrid is the exchange power between IEHS and power grid in season k at time t. When Pk,tgrid>0, IEHS purchases power from power grid. When Pk,tgrid<0, IEHS sells power to power grid. Gk,tNG is the nature gas volume purchased by IEHS from the gas network in season k at time t. mk,tH,b is the volume of hydrogen purchased by IEHS from the hydrogen energy market in season k at time t. cem is the unit price of carbon emission. τgrid, τNG, and τH is the unit carbon emissions of electricity purchase, nature gas purchase, and hydrogen purchase, respectively. comCHP, comHFC, comGB, comP2H, comWT, comPV, comHES, comLHS, and comHHS are unit operation and maintenance costs of CHP, HFC, GB, P2H, WT, PV, HES, and hydrogen storage device. Pk,tCHP,e, Pk,tHFC,e, Pk,tGB,h, Pk,tP2H, Pk,tWT, and Pk,tPV are the output power of CHP, HFC, GB, P2H, WT, and PV, respectively. For the absolute value of nonlinear term in Equation (34), auxiliary variable P¯k,tgrid is introduced to linearize it, and constraint such as P¯k,tgrid≥±Pk,tgrid is added.


## 3.3. Constraint Condition

The bilevel optimization of IEHS shall meet the energy balance constraints of electricity, gas, heat, and hydrogen, as detailed below.(1)Electrical Power Balance ConstraintThe power system model in this paper adopts the linear power flow model in Ref. [33], and the nodal power balance constraint is shown in the following equation.
(37)Pk,tL,e=Pk,tgrid+Pk,tWT+Pk,tPV+Pk,tCHP,e+Pk,tHFC,e−Pk,tP2H.In the formula,Pk,tL,e is the electrical load in season k at time t.(2)Natural Gas Flow Balance Constraint(38)Gk,tmix=Gk,tCHP+Gk,tGBGminNG≤Gk,tmix≤GmaxNG.In the formula,GmaxNG and GminNG are the upper and lower limit of natural gas pipeline flow rate.(3)Thermal Energy Balance Constraint(39)Pk,tL,h=Pk,tCHP,h+Pk,tHFC,h+Pk,tGB,h.In the formula,Pk,tL,h is the heating load in season k at time t.(4)Hydrogen Energy Balance Constraint(40)mk,tL,H+mk,tHFC+mk,tH,inject−mk,tHHS,in−mk,tLHS,in=mk,tP2H+mk,tH,b+mk,tHHS,out+mk,tLHS,out.In the formula,mk,tL,H is the load of hydrogen station in season k at time t. mk,tH,inject is the amount of hydrogen injection in nature gas network in season k at time t.


## 4. Solution of Optimal Model of Integrated Electrothermal and Hydrogen Coupling System

### 4.1. Interval Linear Programming Method

Firstly, the uncertainties of wind power and photovoltaic output, load, electricity price, and gas price are expressed as interval-number according to the method in Section 1.2 of Ref. [34], and wind power and photovoltaic output, load, electricity price, and gas price are defined as interval variables. Then, the lower model Equation (33) is rewritten with interval variables and the general form of interval linear programming model (Equation (41)), according to the method in Section 1.3 of Ref. [34]. Finally, this paper adopts the method of enhanced-interval linear programming proposed in Ref. [35]. After strict mathematical proof, Ref. [35] proposed a two-stage decomposition algorithm, which can directly introduce multiple uncertainties into the interval linear programming, such as Corollary 4 and Corollary 5. Therefore, the rewritten Equation (33) can be decomposed into two sub problems, and the rewritten Equation (33) can be replaced by the optimal submodel and worst submodel. It can provide the optimal value of the objective function and the interval of the optimal value and ensure the absolute feasibility of the final solution space. The CPLEX software needs to be called in the process of solving the model. Due to limited space, the process of specific speculation is detailed in Refs. [34, 35].


#### 4.1.1. General Format of Interval Linear Programming Model

The general format of interval linear programming model in Ref. [34] is shown in the following equation.
(41)minXf=CXs.t.AX≤BX≥0.In the formula,C=cij1×n is the coefficient matrix of the objective function, cij=cij−,cij+. X=xijn×1 is the decision variable of the objective function, xij=xij−,xij+. A=aijm×n is the coefficient matrix of inequality constraint, aij=aij−,aij+. B=bijm×n is the interval value of the constraint, bij=bij−,bij+.The interval linear programming model can be decomposed into optimal submodel and worst submodel by the two-stage decomposition method in Ref. [35], and the uncertainty problem can be transformed into two deterministic problems for solution.(1)Optimal Submodel(42)minf+=∑j=1k10.5cj++cj−xj++∑j=k1+1n0.5cj++cj−xj−s.t.∑j=1k1aij±−Signaij±xj++∑j=k1+1naij±+Signaij±xj−≤bi+,∀ixj±≥0,∀jxj+≥xj,opt−,j=1,2,⋯,k1xj−≤xj,opt+,j=k1+1,k1+2,⋯,n.In the formula, whenj=1,2,⋯,k1, cj±≥0, xj± is the interval variable with positive coefficient in the objective function. When j=k1+1,k1+2,⋯,n, cj±≤0, xj± is the interval variable with negative coefficient in the objective function.(2)Worst Submodel(43)minf−=∑j=1k10.5cj++cj−xj−+∑j=k1+1n0.5cj++cj−xj+s.t.∑j=1k1aij±+Signaij±xj−+∑j=k1+1naij±−Signaij±xj+≤bi−,∀ixj±≥0,∀jxj−≤xj,opt+,j=1,2,⋯,k1xj+≥xj,opt−,j=k1+1,k1+2,⋯,n.In order to ensure that the optimal solutionxj,opt is completely feasible, the constraint Equation (44) is added, as shown below.
(44)∑j=k1−p−1k1−aδj±+xj−−aδj±−xj,opt++∑j=n−q+1naδj±−xj+−aδj±+xj,opt−≤0,∀δ∑j=1k1aδj±−Signaδj±xj,opt++∑j=k1+1naδj±+Signaδj±xj,opt−=bδ+aδj±≤0,j=k1−p−1,⋯,k1aδj±≥0,j=n−q+1,⋯,n.From Equation (42), the upper bound xj,opt+ (j=1,2,⋯,k1) of the optimal value of xj± (j=1,2,⋯,k1), the lower bound xj,opt− (j=k1+1,k1+2,⋯,n) of the optimal value of xj± (j=k1+1,k1+2,⋯,n), and the upper bound f+ of the optimal value of the objective function can be obtained. From Equations (43) and (44), and the lower bound xj,opt− (j=1,2,⋯,k1) of the optimal value of xj± (j=1,2,⋯,k1), the upper bound xj± (j=k1+1,k1+2,⋯,n) of the optimal value of xj± (j=k1+1,k1+2,⋯,n), and the lower bound f− of the optimal value of the objective function can be obtained.Combined with the solution results of the two submodels, the interval value of decision variables, such asx1=xj,opt−,xj,opt+ (j=1,2,⋯,k1), x2=xj,opt−,xj,opt+ (j=k1+1,k1+2,⋯,n), and the interval value of objective function f=f−,f+ can be obtained.


#### 4.1.2. Bilevel Optimization Model Based on Interval Linear Programming Method

According to the models in Part 3 and Section4.1.1, a bilevel optimization model considering multiuncertainty can be established. The model adopts the interval linear programming method to deal with multiuncertainty in the lower model. The objective function interval value Cope=Cope−,Cope+ of lower model can be obtained by optimal and worst submodels constructed, and then the influence of the multiuncertainty on the operation cost of IEHS can be analyzed by the following equation.
(45)Copeaver=0.5Cope−+Cope+Copewidth=0.5Cope−−Cope+.In formula,Copeaver is the expected value of operation cost in IEHS under the influence of uncertain factors. Copewidth is the fluctuation width of the operating cost in IEHS.The lower model will returnCopeaverX to the upper model, and the upper model Equation (31) can be expressed as follows.
(46)maxCtot′=CopeaverX+Cinv.For the bilevel optimization model established in this paper, the upper model is solved by a multiverse optimizer. The lower model uses the interval linear programming method to transform the uncertain equation into two deterministic submodels and then calls the mature commercial solver CPLEX, which can solve the interval linear programming model efficiently and ensure the accuracy of the solution.


### 4.2. Multiverse Optimizer

For the optimal configuration and scheduling of IES, existing researches mostly use traditional optimization algorithms such as particle swarm optimization (PSO) in Ref. [36] and genetic algorithm (GA) in Ref. [37]. In this paper, the multiverse optimizer (MVO) in Ref. [38] is used to solve the configuration capacity of IEHS. This algorithm has simple structure, relatively few parameters and good optimization effect within low-dimension. Therefore, MVO is suitable for solving the planning problem in the upper model. MVO simulates the natural phenomenon that matter in the universe transfers through a wormhole. The main principles are summarized as follows.Firstly, objects in the universe with high expansion rate always tend to objects in the universe with low expansion rate.Secondly, the adjacent universes transfer objects through the white hole/black hole mechanism, and the cosmic position is updated, as shown in the following equation.(47)xij=xkjr1<NIUixijr1≥NIUi.In the formula,xij is the variable j in the universe i. NIUi is the standard expansion rate of the universe i. r1 is a random number between 0 and 1. xkj is the variable j of the universe k selected according to the roulette selection mechanism.Thirdly, objects are transferred between the current optimal universe and the other universes through wormholes, and the position of the universe is updated, as shown in the following equation.(48)xij=Xj+TDRubj−lbjr4+lbjr3<0.5Xj−TDRubj−lbjr4+lbjr3≥0.5r2<WEPxijr2≥WEP.In the formula,xij is the variable j of the universe i. Xj is the variable j of the current optimal universe. ubj and lbj are the maximum and minimum values of the variable j, respectively. r2, r3, and r4 are random numbers between 0 and 1. WEP and TDR are the probability of wormhole existence and travel distance rate, respectively.


### 4.3. Solution Process of Bilevel Optimization Model

The upper model adopts MVO, and the lower model adopts the interval linear programming method to deal with the uncertain problems in IEHS and calls the CPLEX solver to solve the bilevel optimization model. The solving steps are shown in Figure3, which are detailed as follows.Figure 3
Flow chart of solution for the bilevel optimization model.Step 1, input the wind turbine, photovoltaic and load forecasting data of typical days in each season, time of use electricity price, natural gas price, and relevant equipment efficiency and cost parameters. Set the parameters of MVO, including the total number of universes, the maximum number of iterations, the maximum and minimum probability of wormhole existence, and the upper and lower limits of the search space for decision variables. According to the capacity constrain of each equipment, the random initialization of the cosmic population is completed.Step 2, judge whether the convergence condition is met. If not, jump to Step 3. If satisfied, the optimal capacity configuration of each equipment and the minimum total cost of IEHS can be obtained.Step 3, the randomly generated capacity configuration of each equipment is input into the lower model, and take it as the operation constraint of each equipment.Step 4, the interval linear programming method is used to decompose the lower uncertainty model into the optimal and worst deterministic submodels.Step 5, the CPLEX solver is called to solve the optimal submodel. According to the optimal valueCope+, its corresponding solution, and the additional constraint Equation (44), the CPLEX solver is called to solve the worst submodel to obtain Cope−. And then, the operating cost interval of the lower model Cope can be obtained.Step 6, the average valueCopeaver and interval value Copewidth of operation cost can be calculated according to Equation (45).Step 7, the total construction cost of each equipment can be calculated according to the capacity configuration of each equipment randomly generated by the upper model, and the cosmic expansion rate can be calculated with Equation (46), i.e., the objective function value of the upper model. Then, the current optimal universe can be obtained.Step 8, white holes and black holes are generated by roulette selection mechanism, objects are transferred between different universes through white holes/black holes, and the cosmic position is updated according to Equation (47).Step 9, objects are transferred from the current optimal universe to other universes through the wormhole, the position of the universe is updated according to Equation (48), and jump to Step 4.


## 4.1. Interval Linear Programming Method

Firstly, the uncertainties of wind power and photovoltaic output, load, electricity price, and gas price are expressed as interval-number according to the method in Section 1.2 of Ref. [34], and wind power and photovoltaic output, load, electricity price, and gas price are defined as interval variables. Then, the lower model Equation (33) is rewritten with interval variables and the general form of interval linear programming model (Equation (41)), according to the method in Section 1.3 of Ref. [34]. Finally, this paper adopts the method of enhanced-interval linear programming proposed in Ref. [35]. After strict mathematical proof, Ref. [35] proposed a two-stage decomposition algorithm, which can directly introduce multiple uncertainties into the interval linear programming, such as Corollary 4 and Corollary 5. Therefore, the rewritten Equation (33) can be decomposed into two sub problems, and the rewritten Equation (33) can be replaced by the optimal submodel and worst submodel. It can provide the optimal value of the objective function and the interval of the optimal value and ensure the absolute feasibility of the final solution space. The CPLEX software needs to be called in the process of solving the model. Due to limited space, the process of specific speculation is detailed in Refs. [34, 35].


### 4.1.1. General Format of Interval Linear Programming Model

The general format of interval linear programming model in Ref. [34] is shown in the following equation.
(41)minXf=CXs.t.AX≤BX≥0.In the formula,C=cij1×n is the coefficient matrix of the objective function, cij=cij−,cij+. X=xijn×1 is the decision variable of the objective function, xij=xij−,xij+. A=aijm×n is the coefficient matrix of inequality constraint, aij=aij−,aij+. B=bijm×n is the interval value of the constraint, bij=bij−,bij+.The interval linear programming model can be decomposed into optimal submodel and worst submodel by the two-stage decomposition method in Ref. [35], and the uncertainty problem can be transformed into two deterministic problems for solution.(1)Optimal Submodel(42)minf+=∑j=1k10.5cj++cj−xj++∑j=k1+1n0.5cj++cj−xj−s.t.∑j=1k1aij±−Signaij±xj++∑j=k1+1naij±+Signaij±xj−≤bi+,∀ixj±≥0,∀jxj+≥xj,opt−,j=1,2,⋯,k1xj−≤xj,opt+,j=k1+1,k1+2,⋯,n.In the formula, whenj=1,2,⋯,k1, cj±≥0, xj± is the interval variable with positive coefficient in the objective function. When j=k1+1,k1+2,⋯,n, cj±≤0, xj± is the interval variable with negative coefficient in the objective function.(2)Worst Submodel(43)minf−=∑j=1k10.5cj++cj−xj−+∑j=k1+1n0.5cj++cj−xj+s.t.∑j=1k1aij±+Signaij±xj−+∑j=k1+1naij±−Signaij±xj+≤bi−,∀ixj±≥0,∀jxj−≤xj,opt+,j=1,2,⋯,k1xj+≥xj,opt−,j=k1+1,k1+2,⋯,n.In order to ensure that the optimal solutionxj,opt is completely feasible, the constraint Equation (44) is added, as shown below.
(44)∑j=k1−p−1k1−aδj±+xj−−aδj±−xj,opt++∑j=n−q+1naδj±−xj+−aδj±+xj,opt−≤0,∀δ∑j=1k1aδj±−Signaδj±xj,opt++∑j=k1+1naδj±+Signaδj±xj,opt−=bδ+aδj±≤0,j=k1−p−1,⋯,k1aδj±≥0,j=n−q+1,⋯,n.From Equation (42), the upper bound xj,opt+ (j=1,2,⋯,k1) of the optimal value of xj± (j=1,2,⋯,k1), the lower bound xj,opt− (j=k1+1,k1+2,⋯,n) of the optimal value of xj± (j=k1+1,k1+2,⋯,n), and the upper bound f+ of the optimal value of the objective function can be obtained. From Equations (43) and (44), and the lower bound xj,opt− (j=1,2,⋯,k1) of the optimal value of xj± (j=1,2,⋯,k1), the upper bound xj± (j=k1+1,k1+2,⋯,n) of the optimal value of xj± (j=k1+1,k1+2,⋯,n), and the lower bound f− of the optimal value of the objective function can be obtained.Combined with the solution results of the two submodels, the interval value of decision variables, such asx1=xj,opt−,xj,opt+ (j=1,2,⋯,k1), x2=xj,opt−,xj,opt+ (j=k1+1,k1+2,⋯,n), and the interval value of objective function f=f−,f+ can be obtained.


### 4.1.2. Bilevel Optimization Model Based on Interval Linear Programming Method

According to the models in Part 3 and Section4.1.1, a bilevel optimization model considering multiuncertainty can be established. The model adopts the interval linear programming method to deal with multiuncertainty in the lower model. The objective function interval value Cope=Cope−,Cope+ of lower model can be obtained by optimal and worst submodels constructed, and then the influence of the multiuncertainty on the operation cost of IEHS can be analyzed by the following equation.
(45)Copeaver=0.5Cope−+Cope+Copewidth=0.5Cope−−Cope+.In formula,Copeaver is the expected value of operation cost in IEHS under the influence of uncertain factors. Copewidth is the fluctuation width of the operating cost in IEHS.The lower model will returnCopeaverX to the upper model, and the upper model Equation (31) can be expressed as follows.
(46)maxCtot′=CopeaverX+Cinv.For the bilevel optimization model established in this paper, the upper model is solved by a multiverse optimizer. The lower model uses the interval linear programming method to transform the uncertain equation into two deterministic submodels and then calls the mature commercial solver CPLEX, which can solve the interval linear programming model efficiently and ensure the accuracy of the solution.


## 4.1.1. General Format of Interval Linear Programming Model

The general format of interval linear programming model in Ref. [34] is shown in the following equation.
(41)minXf=CXs.t.AX≤BX≥0.In the formula,C=cij1×n is the coefficient matrix of the objective function, cij=cij−,cij+. X=xijn×1 is the decision variable of the objective function, xij=xij−,xij+. A=aijm×n is the coefficient matrix of inequality constraint, aij=aij−,aij+. B=bijm×n is the interval value of the constraint, bij=bij−,bij+.The interval linear programming model can be decomposed into optimal submodel and worst submodel by the two-stage decomposition method in Ref. [35], and the uncertainty problem can be transformed into two deterministic problems for solution.(1)Optimal Submodel(42)minf+=∑j=1k10.5cj++cj−xj++∑j=k1+1n0.5cj++cj−xj−s.t.∑j=1k1aij±−Signaij±xj++∑j=k1+1naij±+Signaij±xj−≤bi+,∀ixj±≥0,∀jxj+≥xj,opt−,j=1,2,⋯,k1xj−≤xj,opt+,j=k1+1,k1+2,⋯,n.In the formula, whenj=1,2,⋯,k1, cj±≥0, xj± is the interval variable with positive coefficient in the objective function. When j=k1+1,k1+2,⋯,n, cj±≤0, xj± is the interval variable with negative coefficient in the objective function.(2)Worst Submodel(43)minf−=∑j=1k10.5cj++cj−xj−+∑j=k1+1n0.5cj++cj−xj+s.t.∑j=1k1aij±+Signaij±xj−+∑j=k1+1naij±−Signaij±xj+≤bi−,∀ixj±≥0,∀jxj−≤xj,opt+,j=1,2,⋯,k1xj+≥xj,opt−,j=k1+1,k1+2,⋯,n.In order to ensure that the optimal solutionxj,opt is completely feasible, the constraint Equation (44) is added, as shown below.
(44)∑j=k1−p−1k1−aδj±+xj−−aδj±−xj,opt++∑j=n−q+1naδj±−xj+−aδj±+xj,opt−≤0,∀δ∑j=1k1aδj±−Signaδj±xj,opt++∑j=k1+1naδj±+Signaδj±xj,opt−=bδ+aδj±≤0,j=k1−p−1,⋯,k1aδj±≥0,j=n−q+1,⋯,n.From Equation (42), the upper bound xj,opt+ (j=1,2,⋯,k1) of the optimal value of xj± (j=1,2,⋯,k1), the lower bound xj,opt− (j=k1+1,k1+2,⋯,n) of the optimal value of xj± (j=k1+1,k1+2,⋯,n), and the upper bound f+ of the optimal value of the objective function can be obtained. From Equations (43) and (44), and the lower bound xj,opt− (j=1,2,⋯,k1) of the optimal value of xj± (j=1,2,⋯,k1), the upper bound xj± (j=k1+1,k1+2,⋯,n) of the optimal value of xj± (j=k1+1,k1+2,⋯,n), and the lower bound f− of the optimal value of the objective function can be obtained.Combined with the solution results of the two submodels, the interval value of decision variables, such asx1=xj,opt−,xj,opt+ (j=1,2,⋯,k1), x2=xj,opt−,xj,opt+ (j=k1+1,k1+2,⋯,n), and the interval value of objective function f=f−,f+ can be obtained.


## 4.1.2. Bilevel Optimization Model Based on Interval Linear Programming Method

According to the models in Part 3 and Section4.1.1, a bilevel optimization model considering multiuncertainty can be established. The model adopts the interval linear programming method to deal with multiuncertainty in the lower model. The objective function interval value Cope=Cope−,Cope+ of lower model can be obtained by optimal and worst submodels constructed, and then the influence of the multiuncertainty on the operation cost of IEHS can be analyzed by the following equation.
(45)Copeaver=0.5Cope−+Cope+Copewidth=0.5Cope−−Cope+.In formula,Copeaver is the expected value of operation cost in IEHS under the influence of uncertain factors. Copewidth is the fluctuation width of the operating cost in IEHS.The lower model will returnCopeaverX to the upper model, and the upper model Equation (31) can be expressed as follows.
(46)maxCtot′=CopeaverX+Cinv.For the bilevel optimization model established in this paper, the upper model is solved by a multiverse optimizer. The lower model uses the interval linear programming method to transform the uncertain equation into two deterministic submodels and then calls the mature commercial solver CPLEX, which can solve the interval linear programming model efficiently and ensure the accuracy of the solution.


## 4.2. Multiverse Optimizer

For the optimal configuration and scheduling of IES, existing researches mostly use traditional optimization algorithms such as particle swarm optimization (PSO) in Ref. [36] and genetic algorithm (GA) in Ref. [37]. In this paper, the multiverse optimizer (MVO) in Ref. [38] is used to solve the configuration capacity of IEHS. This algorithm has simple structure, relatively few parameters and good optimization effect within low-dimension. Therefore, MVO is suitable for solving the planning problem in the upper model. MVO simulates the natural phenomenon that matter in the universe transfers through a wormhole. The main principles are summarized as follows.Firstly, objects in the universe with high expansion rate always tend to objects in the universe with low expansion rate.Secondly, the adjacent universes transfer objects through the white hole/black hole mechanism, and the cosmic position is updated, as shown in the following equation.(47)xij=xkjr1<NIUixijr1≥NIUi.In the formula,xij is the variable j in the universe i. NIUi is the standard expansion rate of the universe i. r1 is a random number between 0 and 1. xkj is the variable j of the universe k selected according to the roulette selection mechanism.Thirdly, objects are transferred between the current optimal universe and the other universes through wormholes, and the position of the universe is updated, as shown in the following equation.(48)xij=Xj+TDRubj−lbjr4+lbjr3<0.5Xj−TDRubj−lbjr4+lbjr3≥0.5r2<WEPxijr2≥WEP.In the formula,xij is the variable j of the universe i. Xj is the variable j of the current optimal universe. ubj and lbj are the maximum and minimum values of the variable j, respectively. r2, r3, and r4 are random numbers between 0 and 1. WEP and TDR are the probability of wormhole existence and travel distance rate, respectively.


## 4.3. Solution Process of Bilevel Optimization Model

The upper model adopts MVO, and the lower model adopts the interval linear programming method to deal with the uncertain problems in IEHS and calls the CPLEX solver to solve the bilevel optimization model. The solving steps are shown in Figure3, which are detailed as follows.Figure 3
Flow chart of solution for the bilevel optimization model.Step 1, input the wind turbine, photovoltaic and load forecasting data of typical days in each season, time of use electricity price, natural gas price, and relevant equipment efficiency and cost parameters. Set the parameters of MVO, including the total number of universes, the maximum number of iterations, the maximum and minimum probability of wormhole existence, and the upper and lower limits of the search space for decision variables. According to the capacity constrain of each equipment, the random initialization of the cosmic population is completed.Step 2, judge whether the convergence condition is met. If not, jump to Step 3. If satisfied, the optimal capacity configuration of each equipment and the minimum total cost of IEHS can be obtained.Step 3, the randomly generated capacity configuration of each equipment is input into the lower model, and take it as the operation constraint of each equipment.Step 4, the interval linear programming method is used to decompose the lower uncertainty model into the optimal and worst deterministic submodels.Step 5, the CPLEX solver is called to solve the optimal submodel. According to the optimal valueCope+, its corresponding solution, and the additional constraint Equation (44), the CPLEX solver is called to solve the worst submodel to obtain Cope−. And then, the operating cost interval of the lower model Cope can be obtained.Step 6, the average valueCopeaver and interval value Copewidth of operation cost can be calculated according to Equation (45).Step 7, the total construction cost of each equipment can be calculated according to the capacity configuration of each equipment randomly generated by the upper model, and the cosmic expansion rate can be calculated with Equation (46), i.e., the objective function value of the upper model. Then, the current optimal universe can be obtained.Step 8, white holes and black holes are generated by roulette selection mechanism, objects are transferred between different universes through white holes/black holes, and the cosmic position is updated according to Equation (47).Step 9, objects are transferred from the current optimal universe to other universes through the wormhole, the position of the universe is updated according to Equation (48), and jump to Step 4.


## 5. Example Analysis

### 5.1. Basic Data

Based on the frame structure of integrated electrothermal and hydrogen coupling system shown in Figure2, the example analysis is carried out in Part 5 considering three load demands, such as electrical load (Eload), heating load (Hload), and hydrogen load (H2load). Three typical days in summer, transitional season, and winter are selected, with 24 hours as the scheduling period and 1 hour as the scheduling interval. The curves of wind turbine (WT), photovoltaic (PV), and load in each typical day are shown in Figure 4. The data of hydrogen load in hydrogen station comes from Ref. [33], and the rest data comes from an actual power system in North China. The data of energy price is shown in Table 1. The planning period is 10 years, and the discount rate is 5%. The hydrogen injection ratio in summer, transitional season, and winter is set to 10%, 14%, and 18%, respectively, such as in Ref. [13]. Other parameters of the example are shown in Table 2, data sources can refer to Ref. [2], Ref. [4], Ref. [13], and Ref. [33].Figure 4
Prediction curves of renewable energy power generation and multiple loads on typical days in each season. (a) Summer. (b) Transitional season. (c) Winter.
(a)(b)(c)Table 1
The numbers of energy prices.
ElectricityTime periodPricePeak load period11:00~15:00, 19:00~22:000.191 $/kWhMedium load period8:00~10:00, 16:00~18:000.131 $/kWhValley load period1:00~7:00, 23:00~24:000.055 $/kWhNature gas0.362 $/m3Hydrogen energy6.11 thousand dollar/tonTable 2
Parameters used in the example analysis.
ParametersValueParametersValueηP2H/ηcomp/ηeHFC/ηhHFC0.6/0.7/0.6/0.2ηeCHP/ηhCHP/ηhGB0.3/0.5/0.75HHVH/HHVNG (MJ/m3)12.75/41.04τH (kg/kg)25b (m3/mol)0.027 × 10-3a (pa·m6/Mol)0.025R (J/mol/K)8.314pminHHS/pmaxHHS (MPa)3/20comCHP/comGB ($/kW)0.035/0.008pminLHS/pmaxLHS (MPa)0.3/3comP2H/comHES ($/kW)0.023/0.008cem (thousand dollar/ton)0.018comHHS/comLHS ($/kW)0.011/0.011CCHP/CGB/CP2H/CHES ($/kW)1370/130/840/25comHFC/comWT/comPV ($/kW)0.011/0.011/0.011CHFC ($/kWh)460τgrid/τNG (kg/kWh)0.137/0.237CLHS/CHHS ($/m3)0.458/0.763


### 5.2. Analysis of Optimization Results on Typical Days

Because the operation scheduling of each equipment in different seasons is similar, Section5.2 selects a typical day in summer to analyze the operation scheduling of each equipment in IEHS, as shown in Figure 5.Figure 5
Optimal operation scheduling results of typical day in summer. (a) Electrical energy. (b) Thermal energy.
(a)(b)It can be seen from Figure5(a) that there was low power demand during the period of 1~8 and 23~24, and the electrical energy is mainly supplied by wind turbine (WT). The surplus power of wind turbine can be converted into hydrogen energy by P2H, which can effectively reduce the wind power curtailment in IEHS. During the period of 7~19, the power demand shall be met by wind turbine and photovoltaic (PV) first, and the residual power demand shall be met by CHP and power purchase from utility grid (UG). As can be seen from Figure 5(b), due to the low thermal energy demand in summer, CHP is shutdown during the period of 7~8, and the gas-fired boiler (GB) is used to meet the thermal energy demand. In other periods, CHP and heat energy storage device (HES) jointly meet the thermal energy demand. HES releases thermal energy during the period of 11~13 and 17~20 and stores thermal energy in other periods. This operation mode is conducive to realize peak load shifting, reduce the output of CHP, and improve the operation economy of IEHS. Through the above analysis, it can be seen that IEHS can integrate a variety of energy and energy consuming equipment, give full play to the time and space complementary advantages of the operation characteristics of energy consuming equipment, and realize the coordination and optimization of multienergy.In Section5.2, the PSO algorithm and MVO algorithm are used to solve the upper model, respectively. Through comparison, the advantages of the MVO algorithm are explained, and the iteration process is shown in Figure 6. It can be seen from Figure 6 that the MVO algorithm is better than the PSO algorithm in convergence speed and convergence ability.Figure 6
Iterative process of IEHS total cost optimization.The results of the two algorithms are shown in Table3. The difference between the configuration schemes of the two algorithms is mainly reflected in the capacity of energy storage devices. The total cost of IEHS in MVO algorithm is 0.2 million dollars less than that of the PSO algorithm, and the total system cost is reduced by 3.58%.Table 3
Comparison of configuration schemes in IEHS between PSO and MVO.
Algorithm nameCHP (MW)GB (MW)P2H (MW)HFC (MWh)HES (MW)LHS (m3)HHS (m3)Total cost (million dollar)PSO2.061.961.6700.225610305.578MVO2.081.84200.5300155.378


### 5.3. Analysis of Optimization Results considering Uncertain Factors

Section5.3 continues to select typical day in summer as the research object, sets a certain fluctuation proportion for multiple uncertain factors such as wind power, photovoltaic, load, and energy price, and analyzes the impact of different uncertain factors on the configuration capacity of equipment and the total cost of IEHS. The influence degree of different uncertain factors on the optimization results is calculated to assist the decision-maker to make reasonable judgment and decision.


#### 5.3.1. Impact of Single Uncertain Factor on IEHS Capacity Configuration and Total Cost

In order to clarify the impact of a single uncertain factor on the capacity configuration and total cost of IEHS, the fluctuation proportions of wind turbine, photovoltaic, load, electricity price, and nature gas price are set separately as ±3%, ±6%, ±9%, and ±12%, and the capacity configuration of each equipment and the total cost of IEHS under different fluctuation proportions are analyzed, as shown in Figure7.Figure 7
The influence of each single uncertain factor on the capacity configuration of IEHS. (a) Fluctuation of renewable energy power generation. (b) Fluctuation of load. (c) Fluctuation of nature gas price. (d) Fluctuation of electricity price.
(a)(b)(c)(d)It can be seen from Figures7(a)–7(c) that with the increase of fluctuation proportion, the installed capacity of CHP, GB, and hydrogen storage device fluctuates the most, and the capacity of other units does not fluctuate greatly. This is mainly because IEHS considers the complementation, substitution and transformation of multienergy, and the main energy consumption in IEHS is natural gas. Therefore, the fluctuation of wind turbine, photovoltaic, load, and nature gas price will have a great impact on the capacity configuration of gas equipment in IEHS. As can be seen from Figure 7(d), the fluctuation of electricity price basically has no impact on the capacity configuration of equipment in IEHS.It can be seen from Table4 that with the increase of the fluctuation ratio of wind power, photovoltaic, and load, the total cost of IEHS shows an upward trend. With the increase of the fluctuation ratio of nature gas price and electricity price, the total cost of IEHS does not change significantly. This shows that the impact of wind turbine, photovoltaic, and load fluctuations on the total cost of IEHS is greater than that of electricity price and nature gas price fluctuations. When the fluctuation ratio is the same, the total cost of IEHS under the fluctuation of electricity price is the smallest.Table 4
The impact of a single uncertain factor on the total cost of IEHS. Unit: million dollar.
Fluctuation ratio, %0±3±6±9±12Renewable energy5.3785.4985.5415.5965.639Load5.3785.4535.4875.5495.618Nature gas price5.3785.3825.3865.395.39Electricity price5.3785.3795.3845.3845.384


#### 5.3.2. Impact of Multiple Uncertain Factors on IEHS Capacity Configuration and Total Cost

Since multiple uncertain factors in IEHS may exist at the same time, in order to make the configuration and operation scheme of IEHS conform to the actual situation, the combined action of multiple uncertain factors is considered in Section5.3.2. The fluctuation proportion of that is set as ±3%, ±6%, ±9%, and ±12%, and the capacity configuration and total cost of IEHS is analyzed under different fluctuation ratios. It can be seen from Figure 8 that under the combined action of multiple uncertain factors, the installed capacity of CHP, GB, and hydrogen storage device is still the most affected, and its fluctuation range is greater than that of any single uncertain factor.Figure 8
The influence of multiple uncertain factors on the capacity configuration of IEHS.It can be seen from Table5, when the combined action of multiple uncertain factors is considered, the total cost of IEHS increases gradually with the increase of the fluctuation ratio of uncertain factors. The greater the fluctuation ratio, the greater the increase of the total cost of IEHS. Compared with any single uncertain factor, the range of total cost fluctuant under multiple uncertain factors is less than the cumulative fluctuation range of total cost under each single uncertain factor. This shows that the effects of multiple uncertain factors on the total cost may cancel each other out, reducing the range of the IEHS total cost fluctuant. Therefore, considering the combined action of multiple uncertain factors in IEHS is more in line with the actual situation, which can make IEHS stronger robustness and better economy.Table 5
The impact of multiple uncertain factors on the total cost of IEHS. Unit: million dollar.
Fluctuation ratio/%0±3±6±9±12Multiple uncertain factors5.3785.395.3445.6315.851


### 5.4. Impact of Multiple Hydrogen Energy Utilization Modes on the Optimization Results

Section5.4 focuses on the impact of multiple hydrogen energy utilization modes on the operation of IEHS through the following three scenarios. Among them, Sections 5.2 and 5.3 are analyzed based on Scenario 3.Scenario 1, the optimal operation of IEHS is realized, without considering multiple utilization modes of hydrogen energy. That is, Scenario 1 only buys hydrogen from the market to meet the hydrogen load demand of hydrogen station, without considering other hydrogen utilization modes.Scenario 2, only HFC mode is considered to realize the optimal operation of IEHS.Scenario 3, consider the multiple utilization modes of hydrogen energy proposed in this paper to realize the optimal operation of IEHS.The change curve of residual capacity in low-pressure hydrogen storage device (LHS) under three scenarios is shown in Figure9. As can be seen from Figure 9, LHS stores hydrogen in transitional season and releases hydrogen in summer and winter under Scenario 2 and Scenario 3, which has obvious seasonal characteristics.Figure 9
Change curve of residual capacity in low-pressure hydrogen storage device under Scenario 2 and Scenario 3.The hydrogen energy utilization in each mode under Scenario 2 and Scenario 3 is shown in Figure10. HSD is the abbreviation of hydrogen storage device. Among them, a positive value represents hydrogen production and a negative value represents hydrogen consumption. The production of hydrogen by alkaline electrolytic water in Scenario 2 and Scenario 3 is 2.78 ton and 9.7 ton, respectively. In order to meet the hydrogen load demand of the hydrogen station, the production of hydrogen in Scenario 3 is larger. In Scenario 2, the hydrogen consumption of the HFC mode occurs in the peak load period of summer and winter. In Scenario 3, the modes of hydrogen injection in nature gas network and hydrogen selling to hydrogen station are adopted, not the HFC mode. The main reason is that the HFC mode has large loss in energy conversion link and poor economy.Figure 10
Results of multiple hydrogen energy utilization modes. (a) Scenario 2. (b) Scenario 3.
(a)(b)The annual operation cost of IEHS under the three scenarios is shown in Table6. It can be seen from Table 6 that the IEHS operation costs of the three scenarios are 0.222 million dollar, 0.219 million dollar, and 0.185 million dollar, respectively, and the power curtailment of renewable energy is 38.41 MWh, 21.78 MWh, and 0 MWh, respectively. Compared with Scenario 1 and Scenario 2, the IEHS operation cost of Scenario 3 is reduced by 16.6% and 15.5%, respectively, and the power curtailment of renewable energy is reduced by 100%, which has significant advantages. It can be seen that multiple utilization modes of hydrogen energy can effectively reduce the IEHS operation cost and improve the capacity of renewable energy utilization. Scenario 1 and Scenario 2 mainly rely on the hydrogen energy market to meet the hydrogen energy demand. While the LHS in Scenario 3 can store a large amount of hydrogen produced by abandon wind and light in transitional season, and LHS can be used in summer and winter to effectively reduce the dependence on the hydrogen energy market and support multiple utilization modes of hydrogen energy with low hydrogen purchase cost. Compared with Scenario 1 and Scenario 2, Scenario 3 has the lowest hydrogen purchase cost, which is only 1281 $. Because the production of hydrogen by alkaline electrolytic water in Scenario 3 is large and needs more electrical energy, its power purchase cost is the highest. Scenario 2 only uses the HFC mode to consume hydrogen for cogeneration, which reduces the power and thermal energy demand of IEHS, resulting in the minimum power purchase cost and gas purchase cost in Scenario 2.Table 6
Operating results under different scenarios.
Number of scenariosPower purchase cost (104 $)Nature gas purchase cost (104 $)Hydrogen purchase cost (104 $)Carbon emission cost (104 $)Operation and maintenance cost (104 $)Total operation cost (104 $)Power curtailment of renewable energy (MWh)15.7388.5595.9741.3670.59822.23638.4125.4148.2855.9741.3350.92221.9321.7837.3748.550.1281.0011.48518.5380


### 5.5. Impact of Hydrogen Purchase Price on the Optimization Results

Section5.5 analyzes the impact of hydrogen purchase price on multiple utilization modes of hydrogen energy, as shown in Figure 11. Among them, a positive value represents hydrogen production, and a negative value represents hydrogen consumption. It can be seen from Figure 11 that with the decrease of hydrogen purchase price, the hydrogen purchases of IEHS increase, the production of hydrogen by alkaline electrolytic water decreases, and the operation cost of IEHS decreases. When the hydrogen purchase price is 1.53 thousand dollar/ton, the operation cost is reduced by 7.6% compared with Scenario 3. Therefore, the hydrogen purchase price has an important impact on the operation cost of IEHS. When the hydrogen purchase price is 4.58 thousand dollar/ton or above, the hydrogen consumption of HFC mode is 0. When the hydrogen purchase price is less than 3.05 thousand dollar/ton, the hydrogen consumption of HFC mode increases, and the economy of HFC model is gradually reflected.Figure 11
The influence of hydrogen purchase price on multiple utilization modes and operation cost of IEHS.


### 5.6. Analysis of Hydrogen Physical Characteristics and Compressor Power Consumption

Section5.6 firstly analyzes the impact of hydrogen physical characteristics on the operation of IEHS through the Scenario 3 and Scenario 4.Scenario 4, the model of high-pressure hydrogen storage device (HHS) is established based on the perfect gas equation, such as in Ref. [2] and Ref. [18].The pressure change curve of HHS and simulation results in each typical day are shown in Figure12 and Table 7. It can be seen from Figure 12 that there has a certain deviation in HHS pressure between the above two scenarios, especially on the April to June in summer with large hydrogen load. Scenario 4 does not consider the hydrogen physical characteristics, and the gas pressure level of HHS is underestimated, which may affect the accurate judgment of the safe operation state of HHS. Compared with Scenario 3, Scenario 4 has a serious misjudgment on the maximum capacity with an error of 13.1%, resulting in overestimation of the IEHS flexibility. Therefore, people mistakenly believe that the hydrogen purchase cost of Scenario 4 has decreased by 40%, and the total cost has decreased by 0.46%. In terms of calculation, the hydrogen physical characteristics are considered in Scenario 3, which increases the complexity of the model. However, the calculation scale of Scenario 3 is similar to that of Scenario 4 through the linearization method, and the calculation time is only increased by 5.9%. Therefore, the model of HHS considering the hydrogen physical characteristics cannot only quickly and accurately simulate the real operation cost of IEHS but also effectively ensure the safe operation of HHS.Figure 12
Pressure change curve of HHS in typical days.Table 7
Comparison of results under Scenario 3 and Scenario 4.
Number of scenariosMaximum capacity of HHS (m3)Hydrogen purchase cost (104 $)Total operation cost (104 $)Computing time (second)3150.12818.53810.14170.07718.4549.5Secondly, Section5.6 quantitatively analyzes the impact of compressor power consumption on the operation of IEHS through the Scenario 3, Scenario 5, and Scenario 6.Scenario 5, the compressor power consumption is calculated as 10% of the hydrogen energy system, such as in Ref. [33].Scenario 6, compressor power consumption is calculated as fixed compression ratio, such as in Ref. [31].The simulation results of three scenarios are shown in Table8, where Ec represents the proportion of compressor power consumption in the power consumption of the hydrogen energy system. It can be seen from Table 8 that the dynamic compression ratio has less impact on the operation cost of IEHS but has a greater impact on the compressor power consumption. Compared with Scenario 5 and Scenario 6, the compressor power consumption in Scenario 3 is reduced by 43.1% and 25.1%, respectively. It can be seen that without considering the dynamic compression ratio, the compressor power consumption will be overestimated.Table 8
Comparison of results of compressor power consumption under different scenarios.
Number of scenariosTotal operation cost (104 $)Compressor power consumption, MWhEc, %318.53830.445.9518.66353.5110618.58540.637.9


## 5.1. Basic Data

Based on the frame structure of integrated electrothermal and hydrogen coupling system shown in Figure2, the example analysis is carried out in Part 5 considering three load demands, such as electrical load (Eload), heating load (Hload), and hydrogen load (H2load). Three typical days in summer, transitional season, and winter are selected, with 24 hours as the scheduling period and 1 hour as the scheduling interval. The curves of wind turbine (WT), photovoltaic (PV), and load in each typical day are shown in Figure 4. The data of hydrogen load in hydrogen station comes from Ref. [33], and the rest data comes from an actual power system in North China. The data of energy price is shown in Table 1. The planning period is 10 years, and the discount rate is 5%. The hydrogen injection ratio in summer, transitional season, and winter is set to 10%, 14%, and 18%, respectively, such as in Ref. [13]. Other parameters of the example are shown in Table 2, data sources can refer to Ref. [2], Ref. [4], Ref. [13], and Ref. [33].Figure 4
Prediction curves of renewable energy power generation and multiple loads on typical days in each season. (a) Summer. (b) Transitional season. (c) Winter.
(a)(b)(c)Table 1
The numbers of energy prices.
ElectricityTime periodPricePeak load period11:00~15:00, 19:00~22:000.191 $/kWhMedium load period8:00~10:00, 16:00~18:000.131 $/kWhValley load period1:00~7:00, 23:00~24:000.055 $/kWhNature gas0.362 $/m3Hydrogen energy6.11 thousand dollar/tonTable 2
Parameters used in the example analysis.
ParametersValueParametersValueηP2H/ηcomp/ηeHFC/ηhHFC0.6/0.7/0.6/0.2ηeCHP/ηhCHP/ηhGB0.3/0.5/0.75HHVH/HHVNG (MJ/m3)12.75/41.04τH (kg/kg)25b (m3/mol)0.027 × 10-3a (pa·m6/Mol)0.025R (J/mol/K)8.314pminHHS/pmaxHHS (MPa)3/20comCHP/comGB ($/kW)0.035/0.008pminLHS/pmaxLHS (MPa)0.3/3comP2H/comHES ($/kW)0.023/0.008cem (thousand dollar/ton)0.018comHHS/comLHS ($/kW)0.011/0.011CCHP/CGB/CP2H/CHES ($/kW)1370/130/840/25comHFC/comWT/comPV ($/kW)0.011/0.011/0.011CHFC ($/kWh)460τgrid/τNG (kg/kWh)0.137/0.237CLHS/CHHS ($/m3)0.458/0.763


## 5.2. Analysis of Optimization Results on Typical Days

Because the operation scheduling of each equipment in different seasons is similar, Section5.2 selects a typical day in summer to analyze the operation scheduling of each equipment in IEHS, as shown in Figure 5.Figure 5
Optimal operation scheduling results of typical day in summer. (a) Electrical energy. (b) Thermal energy.
(a)(b)It can be seen from Figure5(a) that there was low power demand during the period of 1~8 and 23~24, and the electrical energy is mainly supplied by wind turbine (WT). The surplus power of wind turbine can be converted into hydrogen energy by P2H, which can effectively reduce the wind power curtailment in IEHS. During the period of 7~19, the power demand shall be met by wind turbine and photovoltaic (PV) first, and the residual power demand shall be met by CHP and power purchase from utility grid (UG). As can be seen from Figure 5(b), due to the low thermal energy demand in summer, CHP is shutdown during the period of 7~8, and the gas-fired boiler (GB) is used to meet the thermal energy demand. In other periods, CHP and heat energy storage device (HES) jointly meet the thermal energy demand. HES releases thermal energy during the period of 11~13 and 17~20 and stores thermal energy in other periods. This operation mode is conducive to realize peak load shifting, reduce the output of CHP, and improve the operation economy of IEHS. Through the above analysis, it can be seen that IEHS can integrate a variety of energy and energy consuming equipment, give full play to the time and space complementary advantages of the operation characteristics of energy consuming equipment, and realize the coordination and optimization of multienergy.In Section5.2, the PSO algorithm and MVO algorithm are used to solve the upper model, respectively. Through comparison, the advantages of the MVO algorithm are explained, and the iteration process is shown in Figure 6. It can be seen from Figure 6 that the MVO algorithm is better than the PSO algorithm in convergence speed and convergence ability.Figure 6
Iterative process of IEHS total cost optimization.The results of the two algorithms are shown in Table3. The difference between the configuration schemes of the two algorithms is mainly reflected in the capacity of energy storage devices. The total cost of IEHS in MVO algorithm is 0.2 million dollars less than that of the PSO algorithm, and the total system cost is reduced by 3.58%.Table 3
Comparison of configuration schemes in IEHS between PSO and MVO.
Algorithm nameCHP (MW)GB (MW)P2H (MW)HFC (MWh)HES (MW)LHS (m3)HHS (m3)Total cost (million dollar)PSO2.061.961.6700.225610305.578MVO2.081.84200.5300155.378


## 5.3. Analysis of Optimization Results considering Uncertain Factors

Section5.3 continues to select typical day in summer as the research object, sets a certain fluctuation proportion for multiple uncertain factors such as wind power, photovoltaic, load, and energy price, and analyzes the impact of different uncertain factors on the configuration capacity of equipment and the total cost of IEHS. The influence degree of different uncertain factors on the optimization results is calculated to assist the decision-maker to make reasonable judgment and decision.


### 5.3.1. Impact of Single Uncertain Factor on IEHS Capacity Configuration and Total Cost

In order to clarify the impact of a single uncertain factor on the capacity configuration and total cost of IEHS, the fluctuation proportions of wind turbine, photovoltaic, load, electricity price, and nature gas price are set separately as ±3%, ±6%, ±9%, and ±12%, and the capacity configuration of each equipment and the total cost of IEHS under different fluctuation proportions are analyzed, as shown in Figure7.Figure 7
The influence of each single uncertain factor on the capacity configuration of IEHS. (a) Fluctuation of renewable energy power generation. (b) Fluctuation of load. (c) Fluctuation of nature gas price. (d) Fluctuation of electricity price.
(a)(b)(c)(d)It can be seen from Figures7(a)–7(c) that with the increase of fluctuation proportion, the installed capacity of CHP, GB, and hydrogen storage device fluctuates the most, and the capacity of other units does not fluctuate greatly. This is mainly because IEHS considers the complementation, substitution and transformation of multienergy, and the main energy consumption in IEHS is natural gas. Therefore, the fluctuation of wind turbine, photovoltaic, load, and nature gas price will have a great impact on the capacity configuration of gas equipment in IEHS. As can be seen from Figure 7(d), the fluctuation of electricity price basically has no impact on the capacity configuration of equipment in IEHS.It can be seen from Table4 that with the increase of the fluctuation ratio of wind power, photovoltaic, and load, the total cost of IEHS shows an upward trend. With the increase of the fluctuation ratio of nature gas price and electricity price, the total cost of IEHS does not change significantly. This shows that the impact of wind turbine, photovoltaic, and load fluctuations on the total cost of IEHS is greater than that of electricity price and nature gas price fluctuations. When the fluctuation ratio is the same, the total cost of IEHS under the fluctuation of electricity price is the smallest.Table 4
The impact of a single uncertain factor on the total cost of IEHS. Unit: million dollar.
Fluctuation ratio, %0±3±6±9±12Renewable energy5.3785.4985.5415.5965.639Load5.3785.4535.4875.5495.618Nature gas price5.3785.3825.3865.395.39Electricity price5.3785.3795.3845.3845.384


### 5.3.2. Impact of Multiple Uncertain Factors on IEHS Capacity Configuration and Total Cost

Since multiple uncertain factors in IEHS may exist at the same time, in order to make the configuration and operation scheme of IEHS conform to the actual situation, the combined action of multiple uncertain factors is considered in Section5.3.2. The fluctuation proportion of that is set as ±3%, ±6%, ±9%, and ±12%, and the capacity configuration and total cost of IEHS is analyzed under different fluctuation ratios. It can be seen from Figure 8 that under the combined action of multiple uncertain factors, the installed capacity of CHP, GB, and hydrogen storage device is still the most affected, and its fluctuation range is greater than that of any single uncertain factor.Figure 8
The influence of multiple uncertain factors on the capacity configuration of IEHS.It can be seen from Table5, when the combined action of multiple uncertain factors is considered, the total cost of IEHS increases gradually with the increase of the fluctuation ratio of uncertain factors. The greater the fluctuation ratio, the greater the increase of the total cost of IEHS. Compared with any single uncertain factor, the range of total cost fluctuant under multiple uncertain factors is less than the cumulative fluctuation range of total cost under each single uncertain factor. This shows that the effects of multiple uncertain factors on the total cost may cancel each other out, reducing the range of the IEHS total cost fluctuant. Therefore, considering the combined action of multiple uncertain factors in IEHS is more in line with the actual situation, which can make IEHS stronger robustness and better economy.Table 5
The impact of multiple uncertain factors on the total cost of IEHS. Unit: million dollar.
Fluctuation ratio/%0±3±6±9±12Multiple uncertain factors5.3785.395.3445.6315.851


## 5.3.1. Impact of Single Uncertain Factor on IEHS Capacity Configuration and Total Cost

In order to clarify the impact of a single uncertain factor on the capacity configuration and total cost of IEHS, the fluctuation proportions of wind turbine, photovoltaic, load, electricity price, and nature gas price are set separately as ±3%, ±6%, ±9%, and ±12%, and the capacity configuration of each equipment and the total cost of IEHS under different fluctuation proportions are analyzed, as shown in Figure7.Figure 7
The influence of each single uncertain factor on the capacity configuration of IEHS. (a) Fluctuation of renewable energy power generation. (b) Fluctuation of load. (c) Fluctuation of nature gas price. (d) Fluctuation of electricity price.
(a)(b)(c)(d)It can be seen from Figures7(a)–7(c) that with the increase of fluctuation proportion, the installed capacity of CHP, GB, and hydrogen storage device fluctuates the most, and the capacity of other units does not fluctuate greatly. This is mainly because IEHS considers the complementation, substitution and transformation of multienergy, and the main energy consumption in IEHS is natural gas. Therefore, the fluctuation of wind turbine, photovoltaic, load, and nature gas price will have a great impact on the capacity configuration of gas equipment in IEHS. As can be seen from Figure 7(d), the fluctuation of electricity price basically has no impact on the capacity configuration of equipment in IEHS.It can be seen from Table4 that with the increase of the fluctuation ratio of wind power, photovoltaic, and load, the total cost of IEHS shows an upward trend. With the increase of the fluctuation ratio of nature gas price and electricity price, the total cost of IEHS does not change significantly. This shows that the impact of wind turbine, photovoltaic, and load fluctuations on the total cost of IEHS is greater than that of electricity price and nature gas price fluctuations. When the fluctuation ratio is the same, the total cost of IEHS under the fluctuation of electricity price is the smallest.Table 4
The impact of a single uncertain factor on the total cost of IEHS. Unit: million dollar.
Fluctuation ratio, %0±3±6±9±12Renewable energy5.3785.4985.5415.5965.639Load5.3785.4535.4875.5495.618Nature gas price5.3785.3825.3865.395.39Electricity price5.3785.3795.3845.3845.384


## 5.3.2. Impact of Multiple Uncertain Factors on IEHS Capacity Configuration and Total Cost

Since multiple uncertain factors in IEHS may exist at the same time, in order to make the configuration and operation scheme of IEHS conform to the actual situation, the combined action of multiple uncertain factors is considered in Section5.3.2. The fluctuation proportion of that is set as ±3%, ±6%, ±9%, and ±12%, and the capacity configuration and total cost of IEHS is analyzed under different fluctuation ratios. It can be seen from Figure 8 that under the combined action of multiple uncertain factors, the installed capacity of CHP, GB, and hydrogen storage device is still the most affected, and its fluctuation range is greater than that of any single uncertain factor.Figure 8
The influence of multiple uncertain factors on the capacity configuration of IEHS.It can be seen from Table5, when the combined action of multiple uncertain factors is considered, the total cost of IEHS increases gradually with the increase of the fluctuation ratio of uncertain factors. The greater the fluctuation ratio, the greater the increase of the total cost of IEHS. Compared with any single uncertain factor, the range of total cost fluctuant under multiple uncertain factors is less than the cumulative fluctuation range of total cost under each single uncertain factor. This shows that the effects of multiple uncertain factors on the total cost may cancel each other out, reducing the range of the IEHS total cost fluctuant. Therefore, considering the combined action of multiple uncertain factors in IEHS is more in line with the actual situation, which can make IEHS stronger robustness and better economy.Table 5
The impact of multiple uncertain factors on the total cost of IEHS. Unit: million dollar.
Fluctuation ratio/%0±3±6±9±12Multiple uncertain factors5.3785.395.3445.6315.851


## 5.4. Impact of Multiple Hydrogen Energy Utilization Modes on the Optimization Results

Section5.4 focuses on the impact of multiple hydrogen energy utilization modes on the operation of IEHS through the following three scenarios. Among them, Sections 5.2 and 5.3 are analyzed based on Scenario 3.Scenario 1, the optimal operation of IEHS is realized, without considering multiple utilization modes of hydrogen energy. That is, Scenario 1 only buys hydrogen from the market to meet the hydrogen load demand of hydrogen station, without considering other hydrogen utilization modes.Scenario 2, only HFC mode is considered to realize the optimal operation of IEHS.Scenario 3, consider the multiple utilization modes of hydrogen energy proposed in this paper to realize the optimal operation of IEHS.The change curve of residual capacity in low-pressure hydrogen storage device (LHS) under three scenarios is shown in Figure9. As can be seen from Figure 9, LHS stores hydrogen in transitional season and releases hydrogen in summer and winter under Scenario 2 and Scenario 3, which has obvious seasonal characteristics.Figure 9
Change curve of residual capacity in low-pressure hydrogen storage device under Scenario 2 and Scenario 3.The hydrogen energy utilization in each mode under Scenario 2 and Scenario 3 is shown in Figure10. HSD is the abbreviation of hydrogen storage device. Among them, a positive value represents hydrogen production and a negative value represents hydrogen consumption. The production of hydrogen by alkaline electrolytic water in Scenario 2 and Scenario 3 is 2.78 ton and 9.7 ton, respectively. In order to meet the hydrogen load demand of the hydrogen station, the production of hydrogen in Scenario 3 is larger. In Scenario 2, the hydrogen consumption of the HFC mode occurs in the peak load period of summer and winter. In Scenario 3, the modes of hydrogen injection in nature gas network and hydrogen selling to hydrogen station are adopted, not the HFC mode. The main reason is that the HFC mode has large loss in energy conversion link and poor economy.Figure 10
Results of multiple hydrogen energy utilization modes. (a) Scenario 2. (b) Scenario 3.
(a)(b)The annual operation cost of IEHS under the three scenarios is shown in Table6. It can be seen from Table 6 that the IEHS operation costs of the three scenarios are 0.222 million dollar, 0.219 million dollar, and 0.185 million dollar, respectively, and the power curtailment of renewable energy is 38.41 MWh, 21.78 MWh, and 0 MWh, respectively. Compared with Scenario 1 and Scenario 2, the IEHS operation cost of Scenario 3 is reduced by 16.6% and 15.5%, respectively, and the power curtailment of renewable energy is reduced by 100%, which has significant advantages. It can be seen that multiple utilization modes of hydrogen energy can effectively reduce the IEHS operation cost and improve the capacity of renewable energy utilization. Scenario 1 and Scenario 2 mainly rely on the hydrogen energy market to meet the hydrogen energy demand. While the LHS in Scenario 3 can store a large amount of hydrogen produced by abandon wind and light in transitional season, and LHS can be used in summer and winter to effectively reduce the dependence on the hydrogen energy market and support multiple utilization modes of hydrogen energy with low hydrogen purchase cost. Compared with Scenario 1 and Scenario 2, Scenario 3 has the lowest hydrogen purchase cost, which is only 1281 $. Because the production of hydrogen by alkaline electrolytic water in Scenario 3 is large and needs more electrical energy, its power purchase cost is the highest. Scenario 2 only uses the HFC mode to consume hydrogen for cogeneration, which reduces the power and thermal energy demand of IEHS, resulting in the minimum power purchase cost and gas purchase cost in Scenario 2.Table 6
Operating results under different scenarios.
Number of scenariosPower purchase cost (104 $)Nature gas purchase cost (104 $)Hydrogen purchase cost (104 $)Carbon emission cost (104 $)Operation and maintenance cost (104 $)Total operation cost (104 $)Power curtailment of renewable energy (MWh)15.7388.5595.9741.3670.59822.23638.4125.4148.2855.9741.3350.92221.9321.7837.3748.550.1281.0011.48518.5380


## 5.5. Impact of Hydrogen Purchase Price on the Optimization Results

Section5.5 analyzes the impact of hydrogen purchase price on multiple utilization modes of hydrogen energy, as shown in Figure 11. Among them, a positive value represents hydrogen production, and a negative value represents hydrogen consumption. It can be seen from Figure 11 that with the decrease of hydrogen purchase price, the hydrogen purchases of IEHS increase, the production of hydrogen by alkaline electrolytic water decreases, and the operation cost of IEHS decreases. When the hydrogen purchase price is 1.53 thousand dollar/ton, the operation cost is reduced by 7.6% compared with Scenario 3. Therefore, the hydrogen purchase price has an important impact on the operation cost of IEHS. When the hydrogen purchase price is 4.58 thousand dollar/ton or above, the hydrogen consumption of HFC mode is 0. When the hydrogen purchase price is less than 3.05 thousand dollar/ton, the hydrogen consumption of HFC mode increases, and the economy of HFC model is gradually reflected.Figure 11
The influence of hydrogen purchase price on multiple utilization modes and operation cost of IEHS.


## 5.6. Analysis of Hydrogen Physical Characteristics and Compressor Power Consumption

Section5.6 firstly analyzes the impact of hydrogen physical characteristics on the operation of IEHS through the Scenario 3 and Scenario 4.Scenario 4, the model of high-pressure hydrogen storage device (HHS) is established based on the perfect gas equation, such as in Ref. [2] and Ref. [18].The pressure change curve of HHS and simulation results in each typical day are shown in Figure12 and Table 7. It can be seen from Figure 12 that there has a certain deviation in HHS pressure between the above two scenarios, especially on the April to June in summer with large hydrogen load. Scenario 4 does not consider the hydrogen physical characteristics, and the gas pressure level of HHS is underestimated, which may affect the accurate judgment of the safe operation state of HHS. Compared with Scenario 3, Scenario 4 has a serious misjudgment on the maximum capacity with an error of 13.1%, resulting in overestimation of the IEHS flexibility. Therefore, people mistakenly believe that the hydrogen purchase cost of Scenario 4 has decreased by 40%, and the total cost has decreased by 0.46%. In terms of calculation, the hydrogen physical characteristics are considered in Scenario 3, which increases the complexity of the model. However, the calculation scale of Scenario 3 is similar to that of Scenario 4 through the linearization method, and the calculation time is only increased by 5.9%. Therefore, the model of HHS considering the hydrogen physical characteristics cannot only quickly and accurately simulate the real operation cost of IEHS but also effectively ensure the safe operation of HHS.Figure 12
Pressure change curve of HHS in typical days.Table 7
Comparison of results under Scenario 3 and Scenario 4.
Number of scenariosMaximum capacity of HHS (m3)Hydrogen purchase cost (104 $)Total operation cost (104 $)Computing time (second)3150.12818.53810.14170.07718.4549.5Secondly, Section5.6 quantitatively analyzes the impact of compressor power consumption on the operation of IEHS through the Scenario 3, Scenario 5, and Scenario 6.Scenario 5, the compressor power consumption is calculated as 10% of the hydrogen energy system, such as in Ref. [33].Scenario 6, compressor power consumption is calculated as fixed compression ratio, such as in Ref. [31].The simulation results of three scenarios are shown in Table8, where Ec represents the proportion of compressor power consumption in the power consumption of the hydrogen energy system. It can be seen from Table 8 that the dynamic compression ratio has less impact on the operation cost of IEHS but has a greater impact on the compressor power consumption. Compared with Scenario 5 and Scenario 6, the compressor power consumption in Scenario 3 is reduced by 43.1% and 25.1%, respectively. It can be seen that without considering the dynamic compression ratio, the compressor power consumption will be overestimated.Table 8
Comparison of results of compressor power consumption under different scenarios.
Number of scenariosTotal operation cost (104 $)Compressor power consumption, MWhEc, %318.53830.445.9518.66353.5110618.58540.637.9


## 6. Conclusions

Considering the hydrogen physical characteristics, the model of hydrogen storage device is firstly established based on the van Deemter equation in this paper. A linear model of compressor power consumption is proposed secondly, considering dynamic compression ratio. And then, by analyzing multiple uncertain factors of integrated electrothermal and hydrogen coupling system (IEHS) and multiple utilization modes of hydrogen energy, the joint optimization of planning and operation for IEHS are proposed based on multiverse optimizer algorithm (MVO) and interval linear programming method. Finally, the correctness and effectiveness of this method are verified by an example analysis. The main conclusions are summarized as follows.(1)Compared with particle swarm optimization algorithm, the convergence speed and convergence ability of MVO for solving the upper model are better, and the total cost of IEHS is smaller(2)The interval linear programming method can effectively describe the uncertainty of wind turbine, photovoltaic, load, and energy price. On the one hand, the multiple uncertain factors have the greatest impact on the capacity configuration of CHP, GB, and hydrogen storage device. On the other hand, the total cost of IEHS and its fluctuation range considering comprehensive factors is better than the accumulation of that of single factors. Therefore, considering the combined action of multiple uncertain factors in IEHS is more in line with the actual situation compared with any single uncertain factor, considering the combined action of multiple uncertain factors in IEHS is more in line with the actual situation and can reduce the total cost of IEHS and its fluctuation range(3)The multiple utilization modes of hydrogen energy one hand improve the capacity of renewable energy utilization, and the coordination and optimization ability of IEHS. On the other hand, it can effectively deal with the seasonal imbalance between the supply and demand of hydrogen energy, and reduce the operation cost of IEHS(4)The competitiveness of HFC model is closely related to the hydrogen purchase price. When the hydrogen purchase price is higher than 3.05 thousand dollar/ton, this mode will not be adopted(5)Compared with the model of hydrogen storage device (HSD) based on the perfect gas equation, HSD considering the hydrogen physical characteristics can accurately simulate the operating characteristics of high-pressure hydrogen storage device (HHS) and estimate the gas pressure in HHS. Significantly, ignoring the dynamic compression ratio will overestimate the compressor power consumption



---

*Source: 4559960-2023-02-03.xml*