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

,


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 Figure 1.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.

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 Figure 2. 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 2 International Journal of Energy Research (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).
In the frame structure of IEHS shown in Figure 2, 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), highpressure 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   k,t = 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.

× 〠
k∈ sum,tra,win f g In the formula, N k is the proportion of season k in the whole year.N t is the total number of hours in a day.

Hydrogen Storage Link
(1) Hydrogen Storage Device.According to the pressure of hydrogen storage device, which can be divided into highpressure 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 pres-sure 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.
In the formula, p LHS t is the gas pressure of LHS at time t.n LHS t is the amount of substance of hydrogen in LHS at time t.m LHS t is the mass of hydrogen in LHS at time t.R is the ideal gas constant.T H is the gas temperature.M H is the relative molecular weight of hydrogen.V LHS is the LHS volume.According to Equations ( 5) and ( 6), when temperature was constant, p LHS t was linear to m LHS t 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.
In the formula, p HHS t is the gas pressure of HHS at time t.V HHS is the HHS volume.n HHS t 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 n HHS t and mass m HHS t in HHS is shown in Equation (3).
According to Equation (8), p HHS t and m HHS t 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].
International Journal of Energy Research In the formula, N p is the total number of linearized segments.m p,i is the segment point.θ i reflects the position of m HHS t 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).
In the formula, m HS t is the residual capacity of hydrogen storage device at time t, HS ∈ fLHS, HHSg.p HS t is the gas pressure of hydrogen storage device at time t.m HS,in t and m HS,out t are the hydrogen storage rate and hydrogen release rate of hydrogen storage device at time t.m HS,in max and m HS,out max are the maximum hydrogen storage rate and hydrogen release rate.p HS min and p HS max are the minimum and maximum gas pressure of hydrogen storage device.B HS,in t and B HS,out t are the state variables of hydrogen storage and hydrogen release of hydrogen storage device at time t.When B HS,in t = 1, the hydrogen storage device is in hydrogen storage state.When B HS,out t = 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.
In the formula, P CP t is the power consumption of compressor at time t.W CP t is the power consumption of the compressor to compress the unit molar mass of hydrogen at time t.m in t 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.p in and p out are the inlet and outlet pressures of the compressor, respectively, and p out /p in 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 mW x , 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.
In the formula, m and W x are continuous variables.δ x n is a binary variable.τ x n is an auxiliary variable.Equation ( 18) discretizes the continuous variable W x , which is expressed by δ x n .y is the resolution.Equations ( 20) and ( 22) linearize the nonlinear term mW x based on the M method.

Hydrogen Utilization Link
(1) Hydrogen Fuel Cell.The gas pressure of hydrogen supplied to hydrogen fuel cell (HFC) is usually 0.03~0.06MPa 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).(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.
In the formula, ω H k 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, ω H k is generally not higher than 20%, such as in Ref. [13].G mix k,t is the flow of hydrogen, nature gas mixed gas in the natural gas network in season k at time t.G is the hydrogen injection flow in season k at time t .G NG k,t 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.
In formula, HHV mix is the higher heating value of hydrogen, nature gas mixed gas.

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.
In the formula, P

Bilevel Optimization Model of Integrated Electrothermal and Hydrogen Coupling System
According to the frame structure of IEHS shown in Figure 2, 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.

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. min In the formula, C tot is the total cost within the service life of IEHS.C ope is the operation cost of IEHS.C inv is the total  33).Among them, the decision variable is each equipment output.In the formula, P L,e k,t is the electrical load in season k at time t.
(2)Natural Gas Flow Balance Constraint In the formula, G NG max and G NG min are the upper and lower limit of natural gas pipeline flow rate.
(3)Thermal Energy Balance Constraint In the formula, P L,h k,t is the heating load in season k at time t.
(4)Hydrogen Energy Balance Constraint In the formula, m L,H k,t is the load of hydrogen station in season k at time t.m is the amount of hydrogen injection in nature gas network in season k at time t.

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.International Journal of Energy Research [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].

General Format of Interval Linear Programming
Model.The general format of interval linear programming model in Ref. [34] is shown in the following equation.
In the formula, ½C = ð½c ij Þ 1×n is the coefficient matrix of the objective function, The interval linear programming model can be decomposed into optimal submodel and worst submodel by the twostage decomposition method in Ref. [35], and the uncertainty problem can be transformed into two deterministic problems for solution.
(2)Worst Submodel In order to ensure that the optimal solution x j,opt is completely feasible, the constraint Equation ( 44) is added, as shown below.
Combined with the solution results of the two submodels, the interval value of decision variables, such as In formula, C aver ope is the expected value of operation cost in IEHS under the influence of uncertain factors.C width ope is the fluctuation width of the operating cost in IEHS.
The lower model will return C aver ope ðXÞ to the upper model, and the upper model Equation ( 31) can be expressed as follows.
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.

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.
In the formula, x j i is the variable j in the universe i. NI ðU i Þ is the standard expansion rate of the universe i. r 1 is a random number between 0 and 1. x j k 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.
In the formula, x j i is the variable j of the universe i. X j is the variable j of the current optimal universe.ub j and lb j are the maximum and minimum values of the variable j, respectively.r 2 , r 3 , and r 4 are random numbers between 0 and 1. WEP and TDR are the probability of wormhole existence and travel distance rate, respectively.

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 Figure 3, which are detailed as follows.
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 value C + ope , its corresponding solution, and the additional constraint Equation (44), the CPLEX solver is called to solve the worst submodel to obtain C − ope .And then, the operating cost interval of the lower model C ope can be obtained.
Step 6, the average value C aver ope and interval value C width ope 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.International Journal of Energy Research Step 8, white holes and black holes are generated by roulette selection mechanism, 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.

Example Analysis
5.1.Basic Data.Based on the frame structure of integrated electrothermal and hydrogen coupling system shown in Figure 2, 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].

Analysis of Optimization Results on Typical Days.
Because the operation scheduling of each equipment in different seasons is similar, Section 5.2 selects a typical day in summer to analyze the operation scheduling of each equipment in IEHS, as shown in Figure 5.
It can be seen from Figure 5(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 Section 5.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.
The results of the two algorithms are shown in Table 3.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

Generate the initial cosmos of each equipment capacity configuration
Convergence condition is met?
The decision variables of the upper model Eq.39 are input into the lower model Eq.41 and used as constraints According to the upper objective function Eq.54, the expansion rate of each universe is calculated and sorted, and the current optimal universe is further selected.Update the position of each universe using the white hole/black hole transmission mechanism as Eq.55 The wormhole transmission mechanism is used to update the position of each universe as Eq.56 Configuration capacity of each equipment and the total cost of IEHS can be obtained.

Interval linear programming model considering multi-uncertainty
Constructing optimal submodel Constructing worst submodel Constrain Eq.52 The optimal value C ope and corresponding solution are obtained  International Journal of Energy Research is 0.2 million dollars less than that of the PSO algorithm, and the total system cost is reduced by 3.58%.
Analysis of Optimization Results considering Uncertain Factors.Section 5.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.

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   7.
It can be seen from Figures 7(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 Table 4 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    International Journal of Energy Research 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.

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 Section 5.3.2.The fluctuation proportion of that is set as ±3%, ±6%, ±9%, and ±12%,    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.
It can be seen from Table 5, 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.

Impact of Multiple Hydrogen Energy Utilization Modes
on the Optimization Results.Section 5.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 Figure 9.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.
The hydrogen energy utilization in each mode under Scenario 2 and Scenario 3 is shown in Figure 10.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.
The annual operation cost of IEHS under the three scenarios is shown in Table 6.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 14 International Journal of Energy Research 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.

Impact of Hydrogen Purchase Price on the Optimization
Results.Section 5.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    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 Figure 12 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 underes-timated, 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.
Secondly, Section 5.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 Table 8, where E c 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.

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 International Journal of Energy Research 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

Figure 1 :
Figure 1: Logical structure of this thesis.

Figure 2 :
Figure 2: Frame structure of integrated electrothermal and hydrogen coupling system.
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 Section 4.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 C ope = ½C − ope , C + ope of lower model can be obtained by optimal and worst submodels constructed, and then the influence of the multiuncertainty on the operation cost of 8 International Journal of Energy Research IEHS can be analyzed by the following value C ope and interval Value C ope of operation cost based on the Eq.53 aver width

-Figure 3 :
Figure 3: Flow chart of solution for the bilevel optimization model.

Figure 4 :
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.

PowerFigure 5 :
Figure 5: Optimal operation scheduling results of typical day in summer.(a) Electrical energy.(b) Thermal energy.

Figure 6 :
Figure 6: Iterative process of IEHS total cost optimization.

Figure 7 :
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.

Figure 11 : 5 . 6 .
Figure 11: The influence of hydrogen purchase price on multiple utilization modes and operation cost of IEHS.

Figure 12 :
Figure 12: Pressure change curve of HHS in typical days.
is the thermal power of GB in season k at time t.η GB h is the heat generation efficiency of GB.G GB k,t is the gas consumption of GB in season k at time t.PGB,hmax is the maximum thermal power of GB.
max is the maximum electric power of CHP.2.2.2.Gas-Fired Boiler.P GB,h k,t = η GB h HHV mix G GB k,t , 0 ≤ P GB,h k,t ≤ P GB,h max :ð29ÞIn the formula, P GB,h k,t t = 1, HES is in the heat release state.
CHP , C GB , C P2H , C HFC , C HES , C LHS , and C HHS are the unit construction cost of CHP, GB, P2H, HFC, HES, and hydrogen storage device, respectively.The lower model is aimed at minimizing the annual operation cost of the IEHS, including annual energy purchase cost C buy , carbon emission cost C em , and operation and maintenance cost C om , as shown in Equation ( [33]he formula, N 0 is the planning period.r is the discount rate.N k is the proportion of season k in the whole year.This paper mainly considers three typical seasons, such as summer, transitional season, and winter.k∈fsum,tra,wing.N t is the total number of hours in a day.c 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.P IEHS sells power to power grid.G NG k,t is the nature gas volume purchased by IEHS from the gas network in season k at time t.m H,b k,t is the volume of hydrogen purchased by IEHS from the hydrogen energy market in season k at time t.c em 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.cCHP,andP PV k,t 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 Electrical Power Balance Constraint The 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.
om , c HFC om , c GB om , c P2H om , c WT om , c PV om , c HES om , c LHS om , and c HHS om are unit operation and maintenance costs of CHP, HFC, GB, P2H, WT, PV, HES,

Table 1 :
The numbers of energy prices.

Table 2 :
Parameters used in the example analysis.

Table 3 :
Comparison of configuration schemes in IEHS between PSO and MVO.

Table 4 :
The impact of a single uncertain factor on the total cost of IEHS.Unit: million dollar.the capacity configuration and total cost of IEHS is analyzed under different fluctuation ratios.It can be seen from Figure and

Table 6 :
Operating results under different scenarios.

Table 7 :
Comparison of results under Scenario 3 and Scenario 4.

Table 8 :
Comparison of results of compressor power consumption under different scenarios.