Two-stage Robust Operation of Electricity-Gas-Heat Integrated Multi-Energy Microgrids Considering Heterogeneous Uncertainties

solving the operation problem with the column and constraint generation (C & CG) algorithm. Finally, case studies are conducted to verify our proposed method, demonstrating that it can reduce the multi-energy supply cost while the stepped carbon trading mechanism can also significantly reduce carbon emissions.


Background and motivation
With the emerging low-carbon goal around the world [1,2], the invested capacity of renewables such as wind power and photovoltaic cells is continuously expanding.However, the uncertainty and volatility of renewable outputs pose a significant threat to the reliability and economy of the energy system [3].Besides, traditionally, the power, heat, and natural gas systems operate and plan independently with renewables absorbed by the power system only.Nonetheless, with limited transmission capacity [4] and harsh stability constraints [5], there could be a large amount of waste of energy resources.As a result, it is both urgent and necessary to seek coordinated operation methods for multienergy systems with a high penetration level of renewables [6].
A multi-energy system on the distribution level, which is typically called a multi-energy microgrid (MEMG) [7,8], can enhance holistic operation flexibility and accommodate part of renewable generations [9,10].Still, excessive renewables with intrinsic stochasticity pose a huge threat to the reliability of both the long-term (like yearly) and short-term (like daily) MEMG operations with huge energy waste [11,12].Therefore, some scholars make use of the complementary characteristics of solar energy and wind to form a wind-photovoltaic complementary system [13], which can improve energy efficiency.However, the wind-photovoltaic complementary system depends on Efficiency of P2HH T s/r.min /T s/r,max Upper/lower temperature of supply/return heat pipeline ms b/ mr b Water flowing rate of supply/return pipeline y r,min/ y r,max Upper and lower limit of gas source r output π i,min /π i,max Upper and lower limit of node pressure i of the gas system C ij Constants related to temperature, length, diameter, friction, etc. in pipelines.δ e /δ g The carbon quota coefficients for electricity/natural gas.
η GT /η e-chp /η h-chp The carbon emission coefficients for GT/CHP solar resources and cannot be applied in places with less solar energy.Then, the wind-storage hybrid system is widely used in the power system with its better flexibility and independence.Luckily, with the promising advancement of hydrogen storage technology, excessive renewables could now be converted into hydrogen that can be stored in tanks for later use on long-term and short-term scales.
Hydrogen energy storage (HES) is an attractive option compared to battery energy storage (BES) due to its ability to enable large-scale and long-term energy storage.It also has a higher energy density than pumped storage [14,15] and compressed air storage [16], is more flexible in deploying and storing and can be used in a wider range of application scenarios, including transport and industry.In this sense, the MEMG operation can be more flexible, reliable, and economical [17].

Literature review and research gap
Numerous studies have been done on the operation of green hydrogen-based MEMGs.Authors in [18] investigated the operation of the electricity-gas-based MEMG, highlighting the efficacy of power-togas (P2G) in enhancing wind power consumption and suppression load fluctuations in terms of peak shaving and valley filling.Similarly, reference [19] explored the integration of power-to-hydrogen (P2H) in the electricity-gas integrated MEMGs, emphasizing its impact on network operations and potential re-dispatch requirements with a steady-state gas flow model.Reference [20] formulated an integrated model for gas and electrical networks to evaluate the technical feasibility and advantages of P2G over shorter timeframes.It also considers the changes in gas characteristics resulting from the P2G scheme.Though the coordinated power and gas or hydrogen dispatch are fulfilled, the above research ignores the heat (heat and/or cooling) energy generated by the P2H process, resulting in a significant heat energy loss in the range of 20%-30% at the temperatures between 60 • C and 80 • C without proper heat recovery mechanisms.
Recently, the utilization of heat recovery in MEMGs has gained attention, especially in harnessing the potential heat storage capacity of pipelines, boilers, and storage tanks within the heat network.Several scholars have introduced integrated solutions that combine hydrogen production with heat recovery, known as power-to-hydrogen-and-heat (P2HH) conversion.Compared to conventional P2H, the P2HH scheme exhibits higher efficiency, particularly in MEMGs with comprehensive networks.In [21], a P2HH model was introduced, capable of generating both electrical and heat energy by controlling the temperature of the alkaline electrolyzer.Reference [22] proposed a control model of a P2HH-based MEMG considering the utilization of waste heat and verified that the P2HH can contribute to the system efficiency enhancement.
Though heat recovery is considered in the P2HH process, they all assume that the predictions of renewable energy generation and loads are accurate enough.This violates the fact that they are stochastic given uncertain weather conditions and random human energy usage behaviors [23].Unlike solar energy whose output is highly predictable because solar radiation and weather patterns are predictable to some extent [24][25][26].In contrast, the wind power output mainly depends on wind speed, wind direction, and turbine efficiency, which are stochastic and are influenced by a variety of factors such as terrain, seasons, and weather.Although there are many forecasting models and techniques that can predict wind speed and direction, the prediction accuracy is usually lower than that of solar energy.To deal with the uncertainties from wind power and load demands, authors in [27] introduced a lowcarbon operation scheme solved by the stochastic programming method, which incorporates P2G and combined cooling, heat, and power (CCHP) technology for increased flexibility and energy efficiency.In [28], a model predictive control integrated robust dispatch method for a regional MEMG with a district heat network was presented, considering multiple flexibilities from sources, networks, and loads to maintain the indoor temperature within a comfortable range and prevent fast fluctuations.In [29], the quasi-dynamic properties of pipeline temperatures The TSRO method was used to optimize reactive power, aiming to coordinate reactive (continued on next page) R. Zhang et al. in the heat network were studied in the modeling of MEMG operation, and a single-stage robust optimization approach was applied to address wind power uncertainty.However, the stochastic programming method in [27] is too optimistic and time-consuming if the number of scenarios increases, while practical challenges may arise in obtaining the probability distributions of all uncertainties to generate required scenarios.In [28], the model prediction method has too strict and high requirements for prediction accuracy, which limits its application in real-world cases; Finally, the single-stage robust optimization employed in [29] is often too conservative, which would not reflect the true operation condition.
To overcome the deficiencies identified above in dealing with uncertainty sources, the two-stage robust optimization (TSRO) method, which not only ensures the operation reliability in extreme cases by only knowing the boundary information of uncertain factors but also reduces the solution conservatism through dynamic iteration, can be applied.In [30], a TSRO approach considering the gas properties and wind power uncertainty for the economy and safety of MEMG operation was proposed.A comprehensive sensitivity analysis was also conducted in terms of its operating parameters.The TSRO method was used in [31] to optimize reactive power, aiming to coordinate reactive power compensations and manage uncertain wind power outputs.The study in [32] proposed a network hardening model for a transmission-gird level electricity-gas coupling system to seek optimal network hardening plans.It aims to minimize the worst-case weighted gas and power load shedding through network hardening plans.The TSRO method is employed with box sets describing the uncertainty of wind power and line interruptions.Furthermore, the role of P2G and storage equipment was also analyzed.The TRSO method is effectively utilized in the above research to reduce the total energy supply cost and enhance the reliability of MEMGs, but they all neglect the critical role that the hydrogen technology can play for a low carbon future, that is, it could be super effective for the reduction of gas emissions as the only byproduct of bring hydrogen is water.
With the proposal of carbon neutrality goals all over the world, apart from reducing costs or maximizing the reliability index in [30]- [32], minimization of carbon emissions, which is commonly in the form of the carbon trading mechanism should also be considered.In [33], a method for optimal dispatching of low-carbon MEMG economy is proposed, which effectively reduces operation costs while limiting the carbon emissions of MEMGs.Furthermore, [34] presented an optimal scheduling method that simultaneously considers low-carbon and economic factors based on the CCHP plant and carbon capture.The significant application value of demands response and ladder-type carbon trading mechanism was also verified.Nonetheless, the research in [33]- [34] failed to handle system uncertainties arising from the renewables or demands.
At last, different energy networks tend to own distinct properties.For example, there is instant electricity transmission from its source to its load side given the fast transmission speed.However, since transmission speeds of the gas and water (transmission media for heat energy) speed in gas and heat networks are relatively low, there can be huge inertia, which makes both networks potential energy storage sources [23].Nevertheless, currently, in the literature, almost no such comprehensive research has been conducted to make full use dynamic characteristics of both networks to enhance the overall operation flexibility and immunize against uncertainties.
Given the above, the summary of the current literature can be shown in Table 1.

Contribution
Given all the above insights, this article presents a TSRO operation method for the MEMG operation considering both the dynamic characteristics of the heat and gas networks.Besides, uncertainties from wind power and loads are also well addressed for both economic and reliable solutions.The key contributions of this study can be given below.
(1) A coordinated operation model is introduced for a MEMG integrated with electricity-gas-heat networks considering the P2HH devices whose byproduct, i.e., heat energy is also fully utilized via a closed-loop water circulation system.Compared with P2H and P2G, The P2HH scheme increases the system operation flexibility and energy efficiency by converting excessive renewable energies.In addition, the dynamic properties of the gas and heat systems are also considered systematically, which serve as the potential storage sources for more flexible and economic operations.
(2) A two-stage robust coordinated dispatch method for the MEMG is formulated, which alleviates all negative effects of diverse uncertainties from wind power and loads.All internal units can be coordinated on different timescales to reduce the cost and avoid too conservative solutions by conventional single-stage robust optimization methods.(3) The ladder-type carbon emission cost mechanism is proposed to reduce gas emissions.The concept of "virtual carbon emission" is used to measure emissions in the trading behavior of power and natural gas.Meanwhile, given the carbon capture capability of methane reactors, subsidies are provided as economic benefits in the operation model, promoting green and low-carbon development.Minimize the worst-case weighted gas and power load shedding through network hardening plans by the TSRO method [33] Low-carbon economy Propose a method for optimal dispatch of lowcarbon MEMG economy Fail to handle system uncertainties rising from the wind power or demands [34] Present an optimal scheduling method that considers low-carbon and economic factors based on the CCHP plant and carbon capture

Paper organization
The rest of this paper is organized as follows.In Section 2, the framework of EH-IES and the district heating system, natural gas system, and hydrogen system mode are presented.In Section 3, the ladder-type carbon trading mechanism combined with carbon sequestration benefits.Section 4 formulates a two-stage coordination dispatch model for EH-IES considering a ladder-type carbon trading mechanism.Section 5 gives the solution method of the two-stage model.In section 6, case studies are analyzed.Our conclusions are shown in Section 7.

System description
The basic structure of the electricity-heat-gas integrated MEMG is demonstrated in Fig. 1.As can be seen, the system is mainly composed of wind turbines (WTs), heat power units (TP), gas turbine units (GT), combined heat and power units (CHP), electrolyzer (EL), methanation reactor (MR) and various energy loads for power, heat, natural gas, hydrogen.The WT, TP, GT, and CHP provide electricity to the customers, while CHP units can generate heat energy simultaneously to meet the heat demands [28].The P2HH system, consisting primarily of EL, hydrogen storage tank (HST), and heat exchanger (HE), converts excessive wind power into hydrogen and stores it in the HS.Meanwhile, it collects the heat energy generated during the electrolysis process to improve holistic energy efficiency.Some of the hydrogen is processed into natural gas through MR and put into the gas system to meet the fuel requirements of CHP and GT units.Finally, the MEMG is also tied to the main gas network and power grid for bi-directional energy transactions [35].

Hydrogen system considering heat recovery
The configuration of the P2HH in Fig. 2, mainly includes the alkaline electrolyzer, a closed water circulation system, and a hydrogen storage tank.
The electric network supplies P DC to an alkaline electrolyzer through a transformer and rectifier with an efficiency of η 1 .During electrolysis, a portion of the power is converted into hydrogen (P H2 ), which is then separated from the remaining electrolyte on the cathode and pumped back into the cell.The remaining power is converted into heat (P heat ).The working process of P2HH involves two processes, namely, electrolysis and heat transfer.
According to the analysis conducted in reference [5], the P2HH model offers the capability to adjust the temperature of the alkaline electrolysis, thereby achieving variable outputs of hydrogen and heat power.The analysis indicates that as the operation temperature rises, the output heat power decreases, while the output hydrogen power increases.Furthermore, the P2HH model proves effective in reducing heat loss in the DHN, enhancing the GT operating flexibility, improving voltage stability, and reducing system operation costs.

1) Electrolysis Process
During the electrolysis process, electricity is converted into hydrogen as well as heat.The relationship between voltage and energy, as well as the relationship between P H2 and P Heat during temperature variations over a given time interval, are described in (1) and ( 2).These equations illustrate the proportional relationship between the distribution of energy between hydrogen and heat, which is influenced by changes in temperature.
2) Heat-transfer process A closed-loop water circulation scheme has been integrated into the P2HH process, consisting primarily of a heat pump and an exchanger designed to eliminate waste heat and enhance EL efficiency.Consequently, apart from the dissipated heat Pdiss, the residual portion H rec will be eliminated by the cold water.Given the losses incurred by the heat pump and exchanger, the heat power into the network is denoted as η 2 H rec .The heat power generated is determined by the temperature and input electricity, as shown in (5).The whole process can be denoted by the lumped capacitance heat dynamic model represented by ( 3)- (5).
τ = RC means the heat time constant; Δt indicates the time step that τ> > Δt [21].Thus, the dynamic heat model in (4) can be converted into the quasi-steady state model in (8).
Therefore, the temperature T is closely linked to H rec ; the variation in H rec results in a change in T, consequently altering the apportion of P H2 and P Heat .Eq. ( 8) presents the electrolyzer working temperature variation that relies on the input electric power and heat power exchanged with the heat exchanger [1].
According to (1)-( 8), the P2HH scheme could potentially generate more heat power than its production, leading to a decrease in temperature.The heat exchanger can absorb the heat power emitted by the electrolyzer and adjust its working temperature.Without considering the heat exchanger losses, we can obtain the following heat exchanger operation states corresponding to the relationships with the working temperature of the electrolyzer.3) Equivalent Circuit for the P2HH scheme:

P t
With the above, by combining the electrolysis process with the heat transfer process, the P2HH model can be derived, represented as a circuit illustrated in Fig. 2.
First and foremost, the efficiency of the P2HH scheme based on the variables associated with AND, DHN, and HST can be defined in (12):

1) Heat Source
In this study, the heat stations mean the combined heat and power units (CHP), which deliver heat energy to the network as well as water pumps.Unlike the TP units, CHP units can produce both electricity and heat.The electricity produced by CHP is fed into the power network, while the heat energy is harnessed to fulfill the corresponding loads.The electricity and heat power generated by CHP units can be expressed as a linear combination of extreme points within a feasible region [36], as shown in Fig. 3.There exists a coupling relationship between the electricity and heat output produced by CHP units denoted by ( 8)- (10).Any point with a feasible region can be represented by the coordinates of the extreme point.
The CHP heat generation is utilized to heat the water flow, which is expressed as (11): The supply temperature should be raised over a threshold while staying below the upper bound to avoid steam formation and guarantee load-serving quality.
2) Temperature Mixing The water flowing into one node would be mixed.With energy conservation law [9], the temperature of the mixed flow can be computed as ( 14)- (15).
It is assumed that the mixing temperature at each node is equal to that at the starting point of the pipelines tied to that node, demonstrated as

3) Heat Interia and Loss
The temperature in water pipes changes gradually as the water flows from the inlet to the outlet.Additionally, heat exchange occurs during the transmission process due to the temperature difference between the water in the network and the environment, which leads to heat loss and a decrease in the temperature of hot water flows.To address this problem, a nodal method is utilized, which takes into account the time delay and heat loss.The approach can be summarized in two steps: firstly, estimating the outlet temperature based on the historical inlet water temperature with the total transmission time; secondly, adjusting the estimated output temperature to account for the impact of heat losses.
To explain the heat inertia, i.e., the above time delay, a vertical section of a pipeline in Fig. 4 is taken as an instance.
In Fig. 4, the blue and green areas represent the volumes of water mass flowing into and out of the pipeline, while the block on the right depicts the water mass that has flowed out by the end of time t.In other words, water will flow out after time χs t, and all water mass completely flows out after ζs t.Consequently, the values of χs t and ζs t could be determined according to (18) Finally, the temperature of the node at the end of pipelines can be calculated as the average temperature of water in the green area depicted in Fig. 4, as illustrated in (21). ( The heat loss by the exchange between the hot water in the pipelines and the external temperature during the transmission is related to the pipeline length, heat transfer coefficient, and water mass flow as (22).

1) Gas Supply
The gas source in the natural gas network primarily includes gas storage stations and wells.The supply of gas from gas sources should not exceed its upper and lower limits as in (23).

2) Node Pressure and Line Pack Constraints
To ensure safe and stable operation, the node pressure of the gas system must be within a reasonable operation range as: When the gas flows in pipelines, due to its slow transmission speed, there is a line pack effect as shown in ( 25)-( 27) which indicates the dynamic characteristic of the gas network

3) Pipeline Flow Rate
For a specific pipeline, the relationship between its flow rate and node pressure could be represented by the Weymouth eq.[29] as (28)

Ladder-type carbon trading mechanism
To limit MEMG carbon emissions and pave the way for a clean energy future, a ladder-type carbon trading mechanism [26] combined with carbon sequestration benefits could then be proposed.The carbon emissions in the MEMG are mainly from two aspects: (1) the purchase of electricity and gas; (2) generators and coupling units.
Since the of electricity does not directly lead to carbon emissions, the "virtual carbon emission" concept [27] can be introduced to measure the amount of carbon emissions in electricity consumption with a reasonable tax, i.e., the excess of carbon emission quotas.Meanwhile, subsidies are provided for the methane reactor's carbon sequestration benefits, which are included as benefits in the operation.

Carbon emission quotas and actual model
For a low-carbon future, a carbon trading mechanism has been formulated that carbon emission allowances are treated as tradable commodities [26].Therefore, governments allocate initial carbon emission quotas to different sources based on certain principles.During the operation, it is essential to purchase emission quotas if the original emissions are beyond the allocated ones.The total initial carbon emission quotas as shown in (29), and the electricity and gas quotas can be formulated as (30).
E 0,e = δ e ∑ T t=1 ( P G,t + P grid,t ) The total actual carbon emission can be expressed as (31), and carbon dioxide generated from electricity and natural gas is represented by (32).
For the carbon trading, the actual carbon emission in the market is shown in (33):

Ladder-type carbon trading mechanism
To further enhance the restriction of carbon emissions, a ladder-type carbon trading mechanism is adopted based on the approach proposed in [26].Different from the widely-used single carbon emission rate, the ladder-type mechanism as shown in (34) divides multiple carbon emission intervals, then, they are priced differently according to the intervals based on the consumption within a certain period.The higher emissions incur higher taxes.

Carbon sequestration
The methane reactor can produce CH4 through a synthesis reaction of hydrogen and CO2.It serves as carbon sequestration and can be subsidized based on the amount as (35).

Objective function
Traditionally, the MEMG optimization is typically carried out based on a single deterministic snapshot.To deal with the uncertainties associated with wind generation and load demand, a two-stage optimization model will be employed within the unpredictability of wind power generation and loads, the robust optimal solution obtained from this study protects potential realizations.Moreover, the formulation of the two-stage robust optimization model is as follows:

1) First-stage Objective Function
The first stage, known as the day-ahead dispatch operation, operates on a relatively long dispatch timescale.Given the inputs, the first stage is optimized using the RO approach to determine the dispatch of the MEMGs(power output of each unit), the carbon emission generated by fossil fuel and natural gas, the transacted electricity quantity declared to the power grid, and the demand response of electricity, heat, and gas load.The first-stage decisions are decided before the revelation of wind generation uncertainties and are indicated as: f buy (t) = c t e P grid,t + c gas F gas,t

2) Second-stage Objective Function
Due to the errors of wind generation in the worst-case scenario (wind curtailment), it is necessary to adjust the flexible sources (mainly gasfired power units, gas turbine units, and CHP units) in the MEMGs to compensate for the mismatched wind power and load demand, which will result in corresponding adjustment costs.Therefore, the adjustment costs primarily include wind curtailment penalty costs, as well as the costs associated with adjusting the flexible resources in the system to account for uncertain variables.

Uncertain sets modeling
The widespread utilization of flexible loads and distributed wind power has introduced significant uncertainties, which have adverse effects on reliable operation.Therefore, it is necessary to deal with the uncertainties.
Since the probability distribution of all uncertainty sources is hard to find, the interval method is employed to characterize the uncertainties as described ( 56) and (57). (57)

Gas flow linearization
The Weymouth equation in ( 53) is a set of quadratic non-convex constraints, which are difficult to solve.First, nonlinear equations are usually more complicated which will highly increase the computational burden and solution time.Second, the nonlinearities in the system operation model may cause the local optimal solutions and difficulties to converge.
Given the natural gas system's slow dynamic operation, the flow direction of natural gas typically remains constant throughout the day and can be forecasted in advance [24], the start and end nodes of the pipeline can be fixed based on the flow direction (assuming the symbols at the start and end nodes of pipeline lg align with a positive direction of gas flow).
Introducing a new variable fpl, the pipeline flow for the natural gas system is obtained as ( 59)-(60): where Π lf,t and Π lf,t represent the square of the pressure at the start and end of nodes of pipeline l g respectively.The quadratic constraint in (60) can be relaxed into an inequality constraint in rotating second-order cone (SOC) form: To transform the rotate SOC constraint into the standard SOC form, the specific steps are given in (62): To guarantee the secure and steady operation of the gas system, it is essential to impose limits on gas node pressure: ( (63)

Column and constraint generation (C&CG) algorithm
To solve the TSRO problem, the C&CG algorithm [22] can be employed.This algorithm decomposes the original problem into a master one and a slave one.Initially, the master problem constraints are relaxed to a certain extent, creating a relaxed master problem [23].Then, decision variable values obtained by the master problem are transferred to the slave one.The slave problem generates scenarios and incorporates additional variables and constraints into the master problem to generate cutting planes.Finally, the above process is repeated iteratively until the objection function values solved by the master and slave problems are equal.
• Master problem: where Where ℤ represents the objective function values of the slave problem; r denotes the current C&CG iteration and l indicates the index of the iteration.The worst scenarios dk,l * within the uncertainty set k, and k is determined from the slave problem in iteration l.The master problem is a MILP model that aims at determining the first-stage decision xr *.
• Slave problem: the slave problem will generate the worst-case scenarios within uncertainty sets and calculate an upper bound by calculating the expectation of Zr c(x r *).
During the alternating solution process, the relaxation of some constraints in the master problem makes the objective function value lead to a lower bound on the total cost of the MEMG.Conversely, the slave problem provides an upper bound on the total operation cost.The upper and lower bounds are iteratively updated in the alternating solution process.
In summary, the C&CG algorithm flowchart is shown in Fig. 5.

Karush-Kuhn-Tucker (KKT) with big-M method
In the TSRO model, the second stage is a bilayer problem which involves both maximization and minimization problems, which is complex to solve.To simplify the solution process, KKT conditions [19] and the big-M method can be used.KKT conditions allow the conversion of an optimization problem with an objective function into multiple equations, and the big-M [9] method linearizes the constraints.Then, the bilayer problem in the second stage can be reformulated as a single-layer maximization one, take eq.( 69)-(70) as an example, which is shown as follows: 6. Case study

System description
To verify the above model, this paper uses a test MEMG which is based on the IEEE 33-bus distribution power system, 7-node DHN, and 6-node natural gas distribution system, as shown in Fig. 6.In the distribution power network, node E1 is connected to the main grid, and nodes E22, E25, and E33 are connected to 700 kW wind turbines.In the natural gas network, nodes N1 and N2 are connected to the gas retailer and natural gas source.Four coupling devices are there: GT1 connects N5 and E6, GT2 connects N6 and E25, CHP connects N6, H1, and E18, and P2HH and methane conversion device connect E2, N1, and H3.
The internal power, heat, gas, and hydrogen loads are shown in Fig. 7.The wind power forecasting errors are assumed to be 20%, and electricity load errors are 10%.The operation cost for wind generation is set to 0.35 ¥/kWh, and the penalty cost for wind curtailment is 0.3 ¥/kWh.The electricity prices are the Time-of-Use ones, which are indicated in Table 2.
The natural gas price is set as 3.5 ¥/Nm3.The carbon trading base price is 0.26 ¥/kg, and the carbon sequestration subsidy of MR is 0.1 ¥/kg.The length of the ladder-type tax interval is 500 kg.The carbon emission factor of the power and gas networks is set at 1 kg/kWh and

Table 2
Electricity price in case studies [23,24].0.6 kg/kWh, respectively.The upper limit for exchanging power is set as 400 kW.The operating parameters of each equipment are given in Table 3.
All subsequent simulations are implemented on a computer with an i7-4510U CPU with 8GB of memory.The proposed model is developed using CPLEX with Matlab R2016a and solved by CPLEX v12.6.0.The total operation horizon is 24 h with a 1-h unit interval.The computation time for solving the proposed optimization model is 208.9 s.

Analysis of dispatch result for the proposed model
The results of the dispatch for three systems are shown in Fig. 8 (a)-(d).
Based on Fig. 8 (a)-(d), the following can be observed: a) During periods 0:00-6:00 and 13:00-16:00, the MEMG has high wind power generation, while the electricity load is relatively low.The main electricity supply comes from renewables, with the excess wind power being absorbed by the EL for hydrogen production.The hydrogen produced is then used to meet the hydrogen loads and converted to natural gas via methanation for injection into the gas network for gas customers.b) During periods 7:00-9:00 and 19:00-21:00, the electricity loads are high while wind generation is relatively low.As a result, the output of TP, GT, and CHP units increases.However, due to the high purchase price of natural gas, some wind power is still converted to natural gas through the EL and MR to meet the gas loads.Besides, the wind power operates close to its forecasted output upper limit, resulting in a significant reduction in curtailment.c) During the entire process, the heat generated by the EL will be utilized for supplying the heat loads in the heat system through a heat exchanger.This utilization of heat reduces the heat output of CHP units, thereby increasing the holistic energy efficiency, system economy, and flexibility of CHP units for responding to wind power fluctuations.
Based on the above results, to assess the impact of the configuration of hydrogen energy storage (HES) considered heat recovery mechanism on the optimal scheduling of the MEMGs, a comparative analysis between the battery energy storage (BES) and HES is carried out based on the proposed TSRO model.The basic parameters of the BES are detailed in [29].
Fig. 9 shows the comparison of the dispatched wind power and heat output of CHP units between the BES and HES, and Table 4 presents the dispatch results of these two modes.It can be observed that compared to BES, the HES improves comprehensive energy efficiency by 4.4% since it combines the capability of multi-energy storage and supply.The utilized wind power is 2.20 × 10 4 kW, which is increased by 0.28 × 10 4 kW when compared with that of BES.The reason for this is that HES converts the excess wind power into hydrogen for storage through the EL and provides a large amount of heat energy to supply the heat load during the electrolysis process.This decreases the heat output and electrical power output of the CHP units, giving more space to enhance the ability to accommodate wind power.The heat loads are completely supplied by CHP units, which increases the electrical power output of the CHP units.The BES maintains the balance of charging and discharging in general, but due to the lack of multi-energy conversion ability, there is still a certain amount of wind power curtailment.In addition, the HES has more advantages in reducing carbon emissions and energy purchasing costs.Since the hydrogen generated by EL can be used to produce natural gas by MR, it reduces the purchase cost of the system.Besides, the MR absorbs part of the CO 2 released by each unit, thus it reduces the carbon emission to a certain extent.(See Fig. 9.) From the aforementioned finding, it is evident that the EL and MR play a crucial role in converting surplus electricity into natural gas and then injecting it into the gas pipelines.The generated heat is transferred to the heat system, coordinating the operation of the three systems.This improves the consumption of renewables in the system and consequently reduces the MEMG operation costs.

Analysis of MR carbon sequestration
To analyze the impact of MR carbon sequestration on the MEMG operation and the consumption of wind power, two contrasting scenarios are set in this study.The simulation results are given in Table 5.
Case 1.The MR process is not considered.
Case 2. The MR process and its benefits are considered.
From Table 4, it can be seen that in Case 2, the hydrogen energy consumption stage is subdivided into hydrogen load supply and MR.In other words, a portion of the hydrogen energy is allocated to supplying the hydrogen load, while the remaining is utilized for the MR process.
The surplus wind energy is harnessed through MR to generate natural gas, which is then integrated into the gas network.This results in a 30% reduction in overall gas purchases for the MEMG.Conversely, in Case 1, since the methanation process is not considered, the hydrogen produced through EL is limited to supply the hydrogen load, which compromises its flexibility.As a result, the overall wind curtailment rate remains high.Although the operation of MR incurs certain maintenance costs, the benefits brought by MR still lead to a decrease in the system's operating costs by 2% compared to Case 1.
The comparison of wind curtailment and natural gas system dispatch results between Case 1 and Case 2 is illustrated in Fig. 10.In Case 2, the gas load is predominantly supplied by the gas source and MR, with a relatively low level of gas purchased from the gas network.However, in Case 1, the gas load is mainly supplied by the gas source and gas purchased from the higher-level gas network, resulting in a significantly higher gas purchase amount.Although the MR incurs certain operational costs in Case 2, the high gas purchase cost in Case 1 leads to a cost reduction of 5293.13 ¥ in Case 2. Comparing the wind curtailment under both cases, it can be observed that, under the same wind power output, the operation of MR in Case 2 significantly reduces the wind curtailment, that is the wind curtailment only occurs during the 1:00-6:00 and 23:00-24:00 periods.This indicates that the MR could greatly enhance the wind power accommodation rate in the MEMG.This indicates that the MR could greatly enhance the wind power accommodation rate in the MEMG.The reason is that when methanation to natural gas is considered, the hydrogen produced by excess wind power, in addition to satisfying the hydrogen load, will be used for the MR process.When the system is further given a certain amount of subsidy, the system will encourage the MR process, which not only reduces the carbon emissions but also realizes the full use of wind power.

Analysis of the impact of carbon trading mechanism
To assess the effectiveness of the ladder-type carbon trading mechanism and further explore the potential for emission reduction in the system, the dispatch results for the following three scenarios are analyzed.In a conventional carbon trading mechanism, we assume that the carbon trading price is fixed at 0.26 ¥/kg.Table 6 presents the comparison of dispatching results in different cases and Fig. 11 shows the output of heat power units, gas turbines, and CHP units.
From Table 6 and Fig. 11, it can be observed that in Case 3, to optimize traditional economic operations, the system prioritizes electricity generation from TP due to their lower operational costs compared to gas turbines and CHP units.However, coal combustion emits more CO 2 compared to natural gas, resulting in a significant increase in CO 2 emissions in Case 3. Additionally, since Case 3 does not consider carbon trading, the system needs to purchase a large amount of carbon emission quotas from the carbon trading market, leading to the highest cost.In Case 4, carbon emissions decreased by 111.62 kg compared to Case 3, and the operational cost decreased by 0.79%.However, the traditional carbon trading mechanism in Case 4 treats the carbon trading price as a constant and only penalizes emissions based on the carbon emission factor, which is not sufficient to limit high-carbon TP.In Case 5, considering the ladder-type carbon emission mechanism, the price of carbon emission quotas increases in steps.This further restricts carbon emissions while ensuring economic efficiency.Therefore, compared to Case 4, Case 5 has an increase in operational cost by 824.47 ¥, while the holistic carbon emissions decrease by 21.79%.

Method comparison
The TSRO dispatch model in this paper focuses on the planning of generation unit output and power purchase in the day-ahead phase to counteract the risk due to the fluctuation of wind power and load demand in the next day's operation phase.To verify that the model is both robust and economical, two sets of robust optimization cases with uncertainty are set up and compared with the deterministic dispatch method and the stochastic dispatch method, respectively.For the deterministic dispatch method, we assume that the wind power predictions and load demand are accurate enough for real-time operation, and constraints under the worst-case scenarios are not considered.The stochastic dispatch method uses 10,000 Monte Carlo sampling scenarios to describe the normal probability distribution of wind power.For all the uncertainty cases, the lower and upper limits of the uncertain sources are shown in Fig. 12.The solution time and operation time are listed in Table 7.
Table 7 presents the comparison of different dispatch methods.It can be observed that the operation cost of the TSRO method is slightly higher than those of the deterministic and stochastic methods in the dayahead stage.As the set of uncertainties grows, the costs will increase further.The day-ahead operation costs for the TSRO method with the highest uncertainty increased by 35.9% and 26.2% compared to the deterministic and stochastic methods.However, the discrepancy between the planned power generation and the actual power generation of the next day due to the forecasting error of the wind power output and load demand needs to be compensated by purchasing power in real time.Although the day-ahead operation cost under the deterministic method is the lowest, the lack of consideration of the uncertainty of the wind power and load demand results in more costs to compensate for mismatched power.In contrast, the TSRO method does not need to purchase large amounts of power for compensation in the real-time stage, and the cost is 13.7% and 10.2% less.The stochastic method, although slightly less expensive than the deterministic method in the real-time phase, is still not as economical as the TSRO method.
The solution time for each of the three methods is 41.2 s, 1125.8 s, and 208.9 s respectively.As 10,000 samples of uncertain variables and the average of optimization solution with 10,000 samples are done, the  computation time of the stochastic method is the longest.The deterministic method has the shortest solution time since it has the least decision variables.The TSRO method takes longer to compute than the deterministic method by requiring multiple iterations to find the optimal solution, with each iteration requiring an update of the objective function.
To further analyze the operation of the system when the uncertainty parameter changes, we set up a comparative analysis under different uncertainty parameters.Fig.13 demonstrates the fluctuation of load and wind power generation under the condition of Γ l/w values of 0, 12, and 24, and the output of heat power units and gas turbines.In this paper, the fluctuation deviation of wind power output is 20%, and the load is 10% of the predicted value.
From Fig. 13, it can be observed that the worst scenario occurs at the lower bound of the predicted wind power, and the upper bound of the load.In the MEMG, the operation cost will be the highest when the wind output reaches the minimum and the load response reaches the maximum, resulting in the occurrence of the worst scenarios.As Γ l and Γ w increase, there will be higher operation costs.When Γ l = Γ w = 12, the worst-case scenario for wind power occurs during the 6-8 and 21-23 time periods, while the worst-case scenario for load occurs during the 7-14 and 15-23 time periods.In this situation, the output of the TP and GT needs to be adjusted to cope with the fluctuations in wind power and load.Therefore, compared to Γ l = Γ w = 0, the output of TP and GT increases at this time to compensate for the power deficit.With an increase in the number of uncertain scenarios, at Γ l = Γ w = 24, the worst-case scenario for wind power extends to the 2-6, 12-17, and 22-24 time periods, while the load remains in the worst-case scenario throughout the dispatch cycle.Therefore, the heat power units and gas turbines need to increase their power output based on these conditions from Γ l = Γ w = 12.

Conclusion
This article focuses on the operation of the MEMGs via the TSRO method with network dynamic characteristics and introduces the ladder-type carbon trading mechanism.The C&CG algorithm and KKT condition are used to solve the model with a case study conducted on a MEMG based on the IEEE 33-bus distribution power system, 7-node heat network, and 6-node gas network.The main conclusions are as follows: 1) The proposed method effectively manages uncertainty across multiple timescales through coordinated units and multiple energy flows.The numerical study indicates that although the TSRO method has the highest cost in the day-ahead dispatch period, the real-time stage cost can be reduced by 0.07 × 10 5 ¥ and 0.05 × 10 5 ¥ compared with the deterministic and stochastic dispatch methods, respectively.It enhances system robustness and facilitates the seamless integration of wind power.2) Compared with BES, the HES can supply part of the heat load and reduce the heat output and electrical power output of the CHP units due to the heat recovery mechanism, resulting in an increase of 4.4% in the comprehensive energy efficiency.
3) The ladder-type carbon trading mechanism is more conducive to limiting carbon emissions.Under this incentive, it not only provides more online space for wind power but also reduces the purchase of natural gas and further lowers the operation cost of the system.On the load side, by increasing the electricity demand of the EL, the consumption of wind power is increased.
In the future, firstly, the impact of the dynamic response of transport network loads such as hydrogen vehicles and electric vehicles on the IES will be considered.Secondly, the seasonal storage characteristics of hydrogen energy storage systems will also be studied.Finally, it is also possible to incorporate the dynamic responses related to heat generation and distribution with time constants ranging from minutes to an hour.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 4 .
Fig. 4. The vertical section of a pipeline.

Case 3 .
Without carbon trading mechanism.
Sets of mixing node S v1 /S v2 inlet set/outlet set of the thermal pipeline.Pmin G,i/Pmax G,i Upper and lower active power limits of generator i Qmin G,i/Qmax G,i Upper and lower reactive power limits of generator i U W /U D Uncertain sets of wind power and electricity load CHP,t / Q CHP,t Power/heat generation of CHP unit at time t δ n,CHP,t Power-heat coefficient of nth corner point Hl v,t The heat load of node v at time t Ts v,t/Tr v,t Supply/return temperature of node v at time t Ts/r v1,t/ Ts/r v2,t inlet and outlet temperatures of the supply/ return water pipe during the time period t Ts b,t/ Tr b,t The supply water issues and return water temperature at the inlet/outlet of pipe b at time t Ts y,t/Tr y,t The inlet/outlet temperature of hot water at the mixing t Gas source r output at time t π i,t Node pressure of node i at time t Δπ i,tVariation of node pressure of node i at time t q ij(t)Line pack between node i and node j at time t qin ij(t ) /qout ij(t) Inlet/outlet gas of node i and node j at time t R.Zhang et al.

Table 1
Summary of literature.

Table 3
Operating parameters of each coupling equipment.
R.Zhang et al.

Table 4
Comparison of dispatching results.

Table 5
Comparison of dispatching results.

Table 6
Results comparisons in different cases.

Table 7
Costs and solution time under different dispatch methods.