Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Bi-level Robust Clearing Framework of Integrated Electricity and Gas Market Considering Robust Bidding of Smart Energy Hubs  PDF

  • Yanqiu Hou
  • Minglei Bao
  • Yi Ding
College of Electrical Engineering, Zhejiang University, Hangzhou, China

Updated:2025-01-22

DOI:10.35833/MPCE.2024.000093

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

With the implementation of the integrated electricity and gas market (IEGM), the smart energy hubs (SEHs) tend to participate in the market clearing for the optimization of the energy purchase portfolio. Meanwhile, the renewable energy is mushrooming at different scales of energy systems, which can introduce utility-level and distribution-level uncertainties to the operation of the IEGM and SEHs, respectively. Considering the impacts of divergent uncertainties, there exist complicated interactions between the IEGM clearing and the robust bidding of SEHs. The lack of consideration of such interactions may lead to inaccurate modeling of the IEGM clearing and cause potential market inefficiency. To handle this, a bi-level robust clearing framework of the IEGM considering the robust bidding of SEHs is proposed, which simultaneously considers the impacts of utility-level and distribution-level uncertainties. The proposed framework is partitioned into two levels. The upper level is the robust clearing mechanism of the IEGM. At this level, the uncertainty locational marginal electricity and gas prices are derived considering the utility-level uncertainties and the uncertainty-based bidding of SEHs. Given the price signals deduced in the upper level, the lower-level robust bidding of the SEH seeks the optimal bidding strategies while hedging against distribution-level uncertainties. To address the proposed framework, an effective algorithm combining column-and-constraint generation (C&CG) algorithm with the best-response decomposition (BRD) algorithm is formulated. The devised algorithm can efficiently solve the individual robust optimization model and coordinate the interaction of two levels. Numerical experiments are carried out to verify the effectiveness of the proposed framework. Moreover, the impacts of uncertainties on the market clearing results along with the optimal biddings of SEHs are further demonstrated within the proposed framework.

I. Introduction

OWING to the proliferation of desirable efficient and highly flexible gas-fired power plants (GPPs), the mutual dependency of the electricity and gas system is largely strengthened [

1]. Naturally, the strong coupling of these two energy systems drives the tight coupling between the electricity and gas markets [2]. The integrated electricity and gas market (IEGM) is emerging and has attracted widespread attention. The cooperation of IEGM is beneficial to optimally allocate energy resources [3] and thereby enhance the overall economic efficiency of energy systems [4]. Some inspiring research efforts have been put into the coupled operation of the electricity and gas markets [5]-[8]. Reference [5] studies the independent operation of the electricity and gas markets and the operational equilibrium is exploited considering the coordination of the two markets. Two types of coordination schemes from the perspective of energy companies are presented in [6] to couple the electricity and gas markets. In contrast with the separate clearing of the electricity and gas markets, several interesting works are conducted to investigate the centralized operation of the two markets to improve the holistic operational and economic efficiency [7], [8]. Reference [7] demonstrates the benefits of the synchronized and co-cleared mechanism of the electricity and gas markets. The pool-based coupled clearing of electricity and gas markets considering the strategic offerings of large energy producers is presented in [8].

With the rapid progress of advanced communication and information technologies, smart energy hubs (SEHs) are becoming important resources in the IEGM [

9]. Managed by the energy hub operator (EHO), the SEH can purchase and import electricity and gas in the IEGM to optimally satisfy terminal energy demands such as electricity, heat, and gas [10]. Currently, extensive efforts have been devoted to incorporating the SEHs into the IEGM to improve the market efficiency and reduce the demand cost [11]-[13]. A distributed framework for the optimal operation of SEHs in the IEGM is proposed in [11] to economically satisfy the multi-energy demand (MED). To exploit the response behavior of MEDs to the varying energy prices, a two-level framework is developed in [12] to incorporate the multi-elasticity in the IEGM clearing. The exploration of the multi-elasticity of SEHs in this work is proven to be conducive to improving the efficiency and welfare of IEGM. Moreover, integrating the carbon price into the energy prices, a decentralized low-carbon coordination scheme of the IEGM and SEH is proposed in [13], which is beneficial to reducing the carbon emission of energy system. However, the aforementioned studies usually utilize deterministic models, where the renewable energy uncertainties are seldom considered.

As the centralized and distributed renewable energies are increasingly deployed in energy systems, the uncertainties can concurrently reside in the IEGM and SEHs, i.e., at the utility and distribution levels [

14]. To investigate the influence of utility-level renewable energy uncertainties on the IEGM, existing research works have adopted different optimization approaches. Reference [15] proposes a scenario-based clearing mechanism of IEGM based on stochastic optimization (SO). In this work, the energy prices in different scenarios are respectively derived and further averaged to define the clearing energy prices. Featuring numerous scenarios to ensure the feasible solution, the SO would yet lead the heavy computation burden [16]. By contrast, the robust optimization (RO) can well ensure the desirable solution hedges against all the realizations of uncertainties within pre-defined uncertainty sets [16]. Along this line, [17] devises an RO-based market clearing mechanism to price the uncertainties. On this basis, a robust pricing framework for energy and ancillary services in the IEGM is proposed in [18]. Thereinto, the price components in the identified worst-case scenarios are added to the base-case electricity and gas prices to quantify uncertainties. Accordingly, the uncertainty locational marginal electricity prices (ULMEPs) and uncertainty locational marginal gas prices (ULMGPs) are formulated. However, these studies mainly investigate the impacts of utility-level uncertainties on the IEGM clearing, where the integration of SEHs and their internal distribution-level uncertainties are not considered.

To better adapt to the distribution-level uncertainties in SEHs, many researchers begin to study the uncertainty-based optimal bidding of SEHs in the IEGM [

19]-[21]. Reference [19] presents a stochastic probabilistic operation model of the energy hub considering the MED response. An optimal bidding strategy is proposed in [20], where the SO is utilized to model the uncertainties of energy market prices and renewable energies. An RO framework for the optimal operation management of SEH is presented in [21], where distributed renewable energy uncertainties in the SEH are modeled with the robust chance-constrained model. However, the previous studies mainly determine the optimal bidding of SEHs by specifying scenarios given the probabilistic distribution of the energy prices or prescribing the uncertain price intervals. Nevertheless, in practice, the market prices are closely related to the IEGM clearing and challenging to predict precisely [13].

Considering the influences of different levels of uncertainties, the IEGM clearing and the bidding of SEHs are complicatedly intertwined [

14]. On the one hand, the utility-level uncertainties will affect the market clearing prices [18], which will have a great impact on the bidding of SEHs [13]. On the other hand, based on the distribution-level uncertainties and market prices, the SEHs will adjust the bidding behavior. This bidding will in turn alter the energy demands and further the clearing result of IEGM [13]. Unfortunately, as far as the authors are aware, existing studies have not yet been carried out to study the mutual influence between the IEGM and SEHs considering different levels of uncertainties. The incomplete consideration of this bi-directional interplay may lead to the inaccurate reflection of the energy supply-demand situation, which will further lead to the loss of market efficiency [14]. Moreover, this limitation will also potentially cause the sub-optimal decision for SEHs [14].

To fill the research gap, a robust clearing framework of the IEGM considering the robust biddings of SEHs is innovatively proposed, which considers both the utility-level and distribution-level uncertainties. The proposed framework includes two levels, where the upper level is the robust clearing of the IEGM and the lower level is the robust bidding of SEHs. For the IEGM clearing at the upper level, the robust pricing mechanism is applied to effectively reflect the influence of uncertainties on the clearing results given the biddings of SEHs. To deal with the distribution-level uncertainties, a robust bidding model is formulated at the lower-level stage to procure the optimal strategies for the SEHs, where the uncertainty market clearing prices are embedded. With the precise modeling of the bi-directional influence of two levels of uncertainties, the proposed framework could more accurately price and allocate resources in the IEGM. Moreover, the biddings of SEHs are also ensured to be optimally cost-effective since the market clearing is properly embedded within the framework. Compared with existing research incorporating the alternating direction method of multipliers (ADMM) to solve the distributed problems [

22], [23], the proposed framework could elaborately embed the clearing prices into the interaction of the IEGM and SEHs. Also, different from solving the bi-level problems with mathematical programming with equilibrium constraints (MPECs) [2], [24], the best-response decomposition (BRD) algorithm in this paper could efficiently coordinate the bi-level nonconvex robust models.

In light of the discussion, the contributions of this paper are two-fold.

1) A bi-level robust clearing framework of the IEGM is innovatively designed to concurrently consider the impacts of utility-level uncertainties and the robust bidding of SEHs with the distribution-level uncertainties. In the proposed framework, the RO-based market clearing model is developed to derive the uncertainty energy prices, which quantify the influence of the utility-level uncertainties and considers the participation of SEHs. At the lower level, a robust bidding model of the SEH is devised to characterize the distribution-level uncertainties in the premise of updated market clearing prices. The framework could more systematically simulate the mutual influence of different levels of uncertainties and more accurately reflect the energy supply-demand situation. Considering that, the effectiveness of the clearing results can be improved.

2) An efficient algorithm combining the BRD algorithm and the column-and-constraint generation (C&CG) algorithm is developed to elaborately resolve the robust clearing model. The upper-level and lower-level models are iteratively solved while the collaboration of these two levels is achieved in a distributed manner.

II. Problem Framework

The robust clearing framework of the IEGM considering the robust bidding of SEHs is presented in Fig. 1. The holistic framework is hierarchically divided into two levels. The IEGM is cleared at the upper level while SEHs optimally bid their energy inputs in the IEGM governed by individual EHO at the lower level.

Fig. 1  Robust clearing framework of IEGM considering robust bidding of SEHs.

At the upper level, the energy system operator (ESO) manages the available resources including deterministic generation resources, e.g., coal-fired power plants (CPPs), GPPs, and gas wells. Besides, as one of the representative types of utility-level uncertainties, the stochastic power generation of wind farms is also considered. The ESO is responsible for the IEGM clearing. Taking into consideration the utility-level uncertainties in the IEGM, a robust electricity market clearing mechanism [

18] is adopted. Given the energy bids of individual SEH, the robust clearing mechanism determines the ULMEPs and ULMGPs [18]. The energy prices can reflect the impact of the utility-level uncertainties on the market clearing result.

At the lower level, the individual EHOs bid their energy inputs in the IEGM based on the energy price signals derived from the upper level. Similarly, a robust bidding model of the SEH is constructed to obtain the optimal energy dispatch portfolio incorporating the distribution-level uncertainties in the SEH. Each SEH includes the combined heat and power (CHP) unit, electric boiler (EB), gas boiler (GB), and electricity storage (ES) system. The distribution-level uncertainties are herein considered to stem from distributed wind farms (DWFs), as indicated by the red dotted lines. Energy inputs in response to the prices, i.e., electricity and gas inputs of each SEH in the IEGM, are accordingly determined and effectively in turn transferred to the upper level. This model well incorporates the clearing result of IEGM into the decision of the optimal bidding.

The whole robust clearing framework will iterate between the two levels until an equilibrium is reached. Thereupon, the proposed framework can well address utility-level uncertainties in the IEGM while simultaneously accounting for the influence of distribution-level uncertainties in the SEH.

III. Mathematical Formulation

The major assumptions made in the proposed framework are summarized as follows.

1) The strategic behaviors of power plants, gas wells, and SEHs are not considered. Under this assumption, the bidding prices of power plant outputs and gas well production are based on their marginal costs. The SEHs also bid their electricity and gas demands according to their endogenous optimal energy portfolios.

2) The clearing mechanism of IEGM is adopted. The electricity and gas market is operated and cleared by a centralized ESO. A similar assumption is also used in [

18] and [25].

3) The linearized lossless DC power flow model is introduced to depict the operating condition of the power network [

26], [27]. And unit commitment (UC) decisions of power plants are assumed to have been made a prior. If a unit is shut down, its output is fixed at zero [16].

4) The simplified linepack-based gas network model is introduced to feature the dynamic characteristics of the gas networks. The approximated gas model can well depict the storage capacity of pipelines while preserving the desirable computing efficiency [

16], [18], [28].

The proposed framework is divided into two sub-models depicted in Fig. 2. The upper-level model for the robust clearing mechanism of the IEGM determines the ULMEP βb,tE and ULMGP βn,tG considering the utility-level uncertainties and the robust bidding of the SEH. Regarding the clearing energy prices, the lower-level model of the EHO decides the optimal electricity input Pe,tin and gas input Ge,tin, which are subsequently updated and transmitted to the upper level.

Fig. 2  Mathematical model of proposed framework.

A. Upper-level Model: Robust Clearing Mechanism of IEGM

The IEGM clearing is based on the adaptive robust security-constrained economic dispatch (SCED) model [

18].

The operational goal of the clearing model is to minimize the total costs at the day-ahead operational stage and the worst-case costs at the real-time regulating stage [

18], as shown in (1). The first four terms pertain to the cost of gas production, the generation cost of CPPs, and the upward and downward reserve costs of power plants, respectively. The real-time regulating costs include upward/downward output adjustments of CPPs and GPPs, electric and gas load shedding penalties as well as the wind curtailment cost.

miny(I)tuCuSSu,t(I)+i𝒞(Ci,tEpi,t(I)+Ci,tUri,tU,(I)+Ci,tDri,tD,(I))(I): day-ahead operational stage+maxz𝒰minx(II)ti𝒞𝒢(Ci+Δpi,t+,(II)+Ci-Δpi,t-,(II))+bCbEΔDb,tE,(II)+n𝒩CnGΔDn,tG,(II)+w𝒲Cw(Pw,tU,(II)-pw,t(II))(II): real-time regulating stage (1)

where the superscript (I) and (II) represent the day-ahead operational stage and real-time regulating stage, respectively; the subscripts i, b, n, w, u, and t are the indices of the power plant, electric bus, gas node, wind farm, gas well, and time slot, respectively; 𝒞, 𝒢, and 𝒲 are the sets of CPPs, GPPs, and wind farms, respectively; and 𝒩 are the sets of electric buses and gas nodes, respectively; Su,t(I) is the production of gas well; pi,t(I) is the output of power plant; ri,tU,(I) and ri,tD,(I) are the upward and downward reserve capacities of power plant, respectively; CuS, Ci,tE, Ci,tU, and Ci,tD are the unit costs of gas production, generation of power plant, upward reserve, and downward reserve, respectively; Δpi,t+,(II) and Δpi,t-,(II) are the upward and downward adjustments of power plant, respectively; ΔDb,tE,(II) and ΔDn,tG,(II) are the electric and gas load shedding, respectively; Pw,tU,(II) is the real-time wind power output; pw,t(II) is the dispatched wind power; Ci+ and Ci- are the unit costs of upward and downward adjustments for the power plant, respectively; CbE and CnG are the unit costs of electric and gas load shedding, respectively; Cw pertains to the unit wind power curtailment penalty; y(I)=[Su,t(I),pi,t(I),ri,tU,(I),ri,tD,(I)] is the first-stage decision variable vector; x(II)=[Δpi,t+,(II),Δpi,t-,(II), ΔDb,tE,(II),ΔDn,tG,(II),pw,t(II)] is the second-stage decision variable vector; and z𝒰=[Pw,tU,(II)] is the vector denoting utility-level uncertainties.

Formula (2) admits a bunch of representative day-ahead operational constraints [

16]. These constraints comprise power capacity constraints considering reserve limits ((2a) and (2b)), ramping constraints ((2c) and (2d)), and power balance equation ((2e)). Additionally, the system reserve requirements ((2f) and (2g)), transmission flow capacity limit ((2h)), and the nodal active power injection constraint ((2i)) are also included.

pi,t(I)+ri,tU,(I)Pimaxi𝒞𝒢 (2a)
pi,t(I)-ri,tD,(I)Pimini𝒞𝒢 (2b)
pi,t(I)-pi,t-1(I)RUii𝒞𝒢 (2c)
pi,t-1(I)-pi,t(I)RDii𝒞𝒢 (2d)
i𝒞𝒢pi,t(I)+w𝒲pw,tF=b/Db,tE+ePe,tin,(I) (2e)
i𝒞𝒢ri,tU,(I)RtU (2f)
i𝒞𝒢ri,tD,(I)RtD (2g)
bHf-bpbinj,(I)Pfmax (2h)
pbinj,(I)=i𝒞(b)𝒢(b)pi,t(I)+w𝒲(b)pw,tF+e(b)Pe,tin,(I)-Db,tE (2i)

where Pimax and Pimin are the maximum and minimum power outputs, respectively; RUi and RDi are the maximum ramping-up and ramping-down power of power plant, respectively; pw,tF is the forecasting power of wind farm; Db,tE is the electric load; the subscript e is the index of SEH; is the set of SEHs; RtU and RtD are the system-wide upward and downward reserve requirements, respectively; pbinj,(I) is the net injection power at electric bus; Hf-b is the power transfer distribution factor; Pfmax is the maximum active power on the line f; and 𝒞(b), 𝒢(b), 𝒲(b), and (b) are the sets of CPPs, GPPs, wind farms, and SEHs connected to bus b, respectively.

The operating condition of the natural gas system is constrained by (3) [

16]. The nodal gas flow balance is ensured by (3a). The production range and nodal pressure of the gas system are constrained by (3b) and (3c), respectively. The gas flow in the pipeline p is related to the gas pressures of its head node (p) and tail node 𝒯(p) by the Weymouth equation (3d). The linepack is calculated by (3e), whose continuity can be guaranteed by (3f). The linepack at the end of the dispatching horizon T is bounded by (3g). Additionally, the terminal pressures and gas flow of compressor q are depicted by (3h) and (3i), respectively. Formula (3j) correlates the gas consumption of GPP with its power output.

u𝒰(n)Su,t(I)+dis(q)=nFq,t(I)+𝒯(p)=nFTp,t(I)=g𝒢(n)Gg,t(I)+Dn,tG+e(n)Ge,tin,(I)+suc(q)=nFq,t(I)+(p)=nFFp,t(I) (3a)
SuminSu,t(I)Sumax (3b)
ΠnminΠn,t(I)Πnmax (3c)
FFp,t(I)+FTp,t(I)22=(1-2xp,t(I))Cp2((Π(p),t(I))2-(Π𝒯(p),t(I))2) (3d)
Lp,t(I)=Kp(Π(p),t(I)+Π𝒯(p),t(I))/2 (3e)
Lp,t(I)=Lp,t-1(I)+(FFp,t(I)-FTp,t(I)) (3f)
Lp,T(I)Lp,0 (3g)
RqminΠdis(q),t(I)Πsuc(q),tIRqmax (3h)
FqminFq,t(I)Fqmax (3i)
Gg,t(I)=Φpg,t(I)/ηg (3j)

where the subscripts g, q, and p are the indices of GPP, compressor, and pipeline, respectively; Fq,t(I) is the gas flow of compressor; Dn,tG and Gg,t(I) are the gas load and the gas consumption of GPP, respectively; FFp,t(I) and FTp,t(I) are the gas flows at nodes (p) and 𝒯(p), respectively; xp,t(I) is the binary variable denoting the direction of gas flow; Πn,t(I) is the gas pressure of gas node; Lp,t(I) is the linepack of pipeline; Cp and Kp are the Weymouth and transient constants of pipeline, respectively; suc(q) and dis(q) represent the suction and discharge nodes of the compressor, respectively; 𝒰(n) is the gas source connected to gas node; Fqmax and Fqmin are the maximum and minimum gas flows through the compressor, respectively; Rqmin and Rqmax are the minimum and maximum compressing ratios, respectively; pg,t(I) is the power generation of GPP; Φ is the electricity-to-gas conversion factor; ηg is the efficiency of GPP; Sumax and Sumin are the maximum and minimun productions of gas well, respectively; and Πnmax and Πnmin are the maximum and minimum gas pressures of gas node, respectively.

The real-time operational constraints regarding the power and natural gas system are shown in (4) and (5) [

16], respectively. Formula (4a) constrains power output adjustments of both CPPs and GPPs. The system power balance considering real-time electric load shedding and wind power curtailment is described by (4b). The real-time transmission power limits and nodal power injection are defined by (4c) and (4d), respectively. Formulas (4e) and (4f) bound the wind power curtailment and electric load shedding, respectively.

-r^i,tD,(I)Δpi,t(II)r^i,tU,(I)i𝒞𝒢 (4a)
i𝒞𝒢Δpi,t(II)+w𝒲(pw,t(II)-pw,tF)+bΔDb,tE,II=0 (4b)
bHf-b(pbinj,(I)+Δpbinj,(II))Pfmax (4c)
Δpbinj,II=i𝒞𝒢Δpi,tII+w𝒲(b)(pw,tII-pw,tF)+ΔDb,tE,II (4d)
0pw,t(II)Pw,tU,(II) (4e)
0ΔDb,tE,(II)Db,tE (4f)

where Δpi,t(II) and Δpbinj,II are the regulating power output of power plant and the power injection at electric bus, respectively; and the symbol  ^ represents the optimized value.

At the natural gas system side, the real-time gas flow balance is depicted by (5a). Formulas (5b)-(5i) are real-time equivalent formulas related to day-ahead ones. Formula (5k) calculates the real-time gas consumption of GPP.

u𝒰(n)Su,t(II)+dis(q)=nFq,t(II)+𝒯(p)=nFTp,t(II)=g𝒢(n)Gg,t(II)+(Dn,tG-ΔDn,tG,(II))+e(n)Ge,tin,(II)+suc(q)=nFq,t(II)+(p)=nFFp,t(II) (5a)
SuminSu,t(II)Sumax (5b)
ΠnminΠn,t(II)Πnmax (5c)
FFp,t(II)+FTp,t(II)22=(1-2xp,t(I))Cp2((Π(p),t(II))2-(Π𝒯(p),t(II))2) (5d)
Lp,t(II)=Kp(Π(p),t(II)+Π𝒯(p),t(II))/2 (5e)
Lp,t(II)=Lp,t-1(II)+(F(p),t(II)-F𝒯(p),t(II)) (5f)
Lp,T(II)Lp,0 (5g)
RqminΠdis(q),t(II)Πsuc(q),t(II)Rqmax (5h)
FqminFq,t(II)Fqmax (5i)
0ΔDn,tG,(II)Dn,tG (5j)
Gg,t(II)=ηg(pg,t(I)+Δpg,t(II)) (5k)

where Δpg,t(II) is the regulating power output of GPP.

After solving the adaptive robust model, the ULMEP and ULMGPs can be derived.

B. Lower-level Model: Robust Bidding of SEHs

Tightly interacting with the energy system at the demand side, each SEH is composed of multiple energy devices [

10] to satisfy the MED. Besides, the booming development of distributed renewable generation units brings distribution-level uncertainties for the SEH [14]. Due to the distribution-level uncertainties, the optimal bidding strategy of the SEH is thus also formulated as an adaptive RO problem. The optimization objective in (6) is to minimize the operational cost of SEH, i.e., the summing day-ahead energy purchasing costs from the IEGM and the degradation cost of the ES while considering the least cost at the worst-case real-time regulating stage.

minye,(I)t[βe(b),tEPe,tin,(I)+βe(n),tGGe,tin,(I)+Ce,tES(de,tES,(I)+ce,tES,(I))]+(I): day-ahead operational stagemaxz𝒟minxe,(II)t[CeELΔELe,t(II)+CeGLΔGLe,t(II)+CeHLΔHLe,t(II)+CeW(Pe,tU,(II)-Pe,tD,(II))](II): real-time regulating stage (6)

where e(b) is the electric bus connected to SEH e; e(n) is the gas node linking the SEH and the natural gas system; ce,tES,(I) and de,tES,(I) are the charging and discharging power of ES at the day-ahead operational stage, respectively; Ce,tES is the unit degradation cost; ΔELe,t(II), ΔGLe,t(II), and ΔHLe,t(II) are the unserved terminal electricity, gas, and heat demands at the real-time regulating stage, respectively, and CeEL, CeGL, and CeHL are their corresponding unserved costs; Pe,tU,(II) and Pe,tD,(II) are the available and dispatched wind power of DWF, respectively; CeW is the unit wind curtailment cost; ye,(I)=[Pe,tin,(I),Ge,tin,(I),de,tES,(I), ge,tES,(I),Pe,tCHP,(I),He,tEB,(I),He,tGB,(I)] is the first-stage decision variable vector, Pe,tCHP,(I) is the electric power output of CHP, and He,tGB,(I) and He,tEB,(I) are the heat energy outputs of GB and EB, respectively;xe,(II)=[ΔELe,t(II),ΔGLe,t(II),ΔHLe,t(II),de,tES,(II),Pe,tCHP,(II), He,tEB,(II),He,tGB,(II),Pe,tD,(II)] is the second-stage decision variable vector; and z𝒟=[Pe,tU,(II)] is the vector denoting distribution-level uncertainties.

The typical day-ahead operational constraints of SEH are described in (7). The multi-energy flow balances are ensured with (7a)-(7c). Formula (7d) limits the power output of CHP as well as the heat outputs of EB and GB. The operating states of ES are depicted with (7e)-(7n) [

29]. Specifically, the net power output of ES is presented by (7e). The exclusiveness of the charging and discharging states is ensured by (7f). In (7g) and (7h), the charging power considering the down- and up-reserve is characterized [29], respectively. Similarly, (7i) and (7j) feature the discharging power considering down- and up-reserve, respectively. The continuity of the state of charge (SOC) of the ES is assumed with (7k). Moreover, the bound of SOC and its consistency requirement are introduced with (7l) and (7m), respectively. The non-negativity of reserves is guaranteed by (7n).

Pe,tin,(I)+Pe,tCHP,(I)+Pe,tf+de,tES,(I)=Pe,tEB,(I)+ce,tES,(I)+ELe,t (7a)
Ge,tin,(I)=1ηtuPe,tCHP,(I)+1ηGBHe,tGB,(I)+GLe,t (7b)
(1-ηtu-ηlo)ηheηtuPe,tCHP,(I)+He,tGB,(I)+He,tEB,(I)=HLe,t (7c)
P̲eCHPPe,tCHP,(I)P¯eCHPH̲eEBHe,tEB,(I)H¯eEBH̲eGBHe,tGB,(I)H¯eGB (7d)
Pe,tES,(I)=de,tES,(I)-ce,tES,(I) (7e)
xe,tc,(I)+xe,td,(I)1    xe,tc,(I),xe,td,(I){0,1} (7f)
0de,tES,(I)-rde,tDN,(I)xed,(I)D¯ES (7g)
0ce,tES,(I)-rce,tUP,(I)xec,(I)C¯ES (7h)
0ce,tES,(I)+rce,tDN,(I)xec,(I)C¯ES (7i)
0de,tES,(I)+rde,tUP,(I)xed,(I)D¯ES (7j)
SOCe,t(I)-SOCe,t-1(I)=ηecce,tES,(I)-de,tES,(I)/ηed (7k)
SOC̲eSOCe,t(I)SOC¯e (7l)
SOCe,T(I)SOCe,0 (7m)
rce,tUP,(I)0rce,tDN,(I)0rde,tUP,(I)0rde,tDN,(I)0 (7n)

where Pe,tf and Pe,tEB,(I) are the forecasting wind power of DWF and the power consumption of EB, respectively; ηtu, ηhe, and ηlo are the efficiencies of turbine, heat, and loss characterizing the operation of CHP [

30], respectively; xe,tc,(I) and xe,td,(I) are the binary variables denoting the charging and discharging states, respectively; rce,tUP,(I) and rce,tDN,(I) are the charging up-reserve and down-reserve, respectively; rde,tUP,(I) and rde,tDN,(I) are the discharging up-reserve and down-reserve, respectively; SOCe,t(I) is the SOC of ES; ηec and ηed are the charging and discharging efficiencies, respectively; ηGB is the heat conversion efficiency of GB; ELe,t, GLe,t, and HLe,t are the terminal electric, gas, and heat loads, respectively; and ()¯ and ()̲ denote the upper limit and lower limit, respectively.

The real-time regulating stage regulates the dispatch of each energy device inside the SEH to hedge against the uncertainty of wind power output. Formulas (8a)-(8c) re-dispatch multi-energy flows. The charging and discharging power considering the down- and up-reserve is depicted by (8f) and (8g), respectively [

29]. Other formulas ((8h)-(8j)) correspond to those at the day-ahead operational stage.

Pe,tin,(I)+Pe,tD,(II)+Pe,tCHP,(II)+de,tES,(II)=Pe,tEB,(II)+ce,tES,(II)+ELe,t-ΔELe,t(II) (8a)
Ge,tin,(I)=1ηtuPe,tCHP,(II)+1ηGBHe,tGB,(II)+GLe,t-ΔGLe,t(II) (8b)
(1-ηtu-ηlo)ηheηtuPe,tCHP,(II)+He,tGB,(II)+He,tEB,(II)=HLe,t-ΔHLe,t(II) (8c)
P̲eCHPPe,tCHP,(II)P¯eCHPH̲eEBHe,tEB,(II)H¯eEBH̲eGBHeGB,(II)H¯eGB0Pe,tD,(II)Pe,tU,(II) (8d)
Pe,tES,(II)=de,tES,(II)-ce,tES,(II) (8e)
de,tES,(I)-rde,tDN,(I)de,tES,(II)de,tES,(I)+rde,tUP,(I) (8f)
ce,tES,(I)-rce,tUP,(I)ce,tES,(II)ce,tES,(I)+rce,tDN,(I) (8g)
SOCe,t(II)-SOCe,t-1(II)=ηecce,tES,(II)-de,tES,(II)/ηed (8h)
SOC̲eSOCe,t(II)SOC¯e (8i)
SOCe,T(II)SOCe,0 (8j)

After solving the robust bidding of SEHs, the electricity and gas inputs, i.e., Pe,tin,(I) and Ge,tin,(I), can be determined and subsequently employed as the input of the upper-level market clearing.

C. Modeling Utility-level and Distribution-level Uncertainties of Wind Power Output

A box-like uncertainty set [

17] is uniformly introduced to depict the stochastic outputs of wind generation both with utility and distribution scales.

z𝒰/𝒟ϒ𝒰/𝒟:=(z˜w/e,t𝒰/𝒟,+,z˜w/e,t𝒰/𝒟,-)w/e(z˜w/e,t𝒰/𝒟,++z˜w/e,t𝒰/𝒟,-)Γ1𝒰/𝒟,t(z˜w/e,t𝒰/𝒟,++z˜w/e,t𝒰/𝒟,-)Γ2𝒰/𝒟,z˜w/e,t𝒰/𝒟,++z˜w/e,t𝒰/𝒟,-1,z˜w/e,t𝒰,𝒟,+,z˜w/e,t𝒰/𝒟,-{0,1} (9a)
Pw,tU,(II)=pw,tF+z˜w,t𝒰,+P¯w,t-z˜w,t𝒰,-P̲w,t (9b)
Pe,tU,(II)=Pe,tf++z˜e,t𝒟,+P¯e,t-z˜e,t𝒟,-P̲e,t (9c)

where the superscripts 𝒰 and 𝒟 denote the corresponding variables regarding the utility level and distribution level, respectively; z˜w/e,t𝒰/𝒟,+ and z˜w/e,t𝒰/𝒟,- are the uncertainty variables regarding the output of wind farm w and SEH e, while the spatial and temporal uncertainty budgets Γ1𝒰/𝒟 and Γ2𝒰/𝒟 constrain the conservativeness; P¯w,t and P̲w,t are the upper and lower bounds for the wind power of wind farm w, respectively; and P¯e,t and P̲e,t are the upper and lower bounds of SEH e, respectively. The actual power generation of wind farms is accordingly represented with (9b).

IV. Solution Methodology

A. Solution Methodology: C&CG Algorithm

The upper-level and lower-level models are both individually constructed with the two-stage adaptive RO model, which is hard to be solved directly [

16]. Moreover, the Weymouth equations depicted in (3d) and (5d) bring the nonconvexity of the upper-level IEGM clearing. To bear these, the second-order cone (SOC) relaxation method is first introduced to recover the convexity of the upper-level problem [31]. By the SOC relaxation, the mixed integer nonconvex constraints (3d) and (5d) are compactly transformed as a bunch of constraints denoted as [16]:

[(FFp,t+FTp,t)/2]2Cp2((Πp,ta)2-(Πp,tb)2)xp,t(Π(p),tmin-Π𝒯(p),tmax)Π(p),t-Πp,taxp,t(Π(p),tmax-Π𝒯(p),tmin)xp,t(Π𝒯(p),tmin-Π(p),tmax)Π𝒯(p),t-Πp,tbxp,t(Π𝒯(p),tmax-Π(p),tmin)(1-xp,t)(Π(p),tmin-Π𝒯(p),tmax)Π(p),t-Πp,tb          (1-xp,t)(Π(p),tmax-Π𝒯(p),tmin)(1-xp,t)(Π𝒯(p),tmin-Π(p),tmax)Π𝒯(p),t-Πp,ta          (1-xp,t)(Π𝒯(p),tmax-Π(p),tmin) (10)

where Πp,ta and Πp,tb are the end nodes of the pipeline with higher and lower gas pressures, respectively. When the gas pressure of head node Π(p),t is larger than that of the tail node Π𝒯(p),t, i.e., xp,t=0, the constraint (10) ensures Πp,ta=Π(p),t and Πp,tb=Π𝒯(p),t. Otherwise, when Π(p),t is lower than Π𝒯(p),t, i.e., xp,t=1, the constraint (10) limits Πp,tb=Π(p),t and Πp,ta=Π𝒯(p),t. By this means, the binary variable xp,t is successfully separated from the original Weymouth equation, which can be further transformed into the convex SOC.

On the other hand, the constraints of the lower-level model are all linear. The linear constraints can be regarded as the special case of SOC optimization.

To effectively solve the adaptive robust models both for upper-level and lower-level problems, a modified SOC-dual-based C&CG algorithm is incorporated.

For ease of interpretation, the robust clearing model of the IEGM and the robust bidding of SEH can both be formulated as the tri-level “min-max-min” problem with a compact form:

minyttatTyt(I): first-level+maxuminxttgtTxt(II): second-  and  third-level (11)

where at and gt are the coefficients regarding the objective functions (1) and (6), respectively; and yt and xt are the decision variable vectors regarding the day-ahead operational stage and real-time regulating stage, respectively.

The constraints of compositive first-level “min” problem are denoted by (12a) and (12b). The linear constraint (12a) integrates all the linear constraints at the first stage, i.e., (2), (3a)-(3c), (3e)-(3j) for the upper-level problem and (7) for the lower-level problem. Formula (12b) encompasses the first-stage SOC constraints, i.e., SOC reformulations of constraint (3d) with (10).

A1ytbt (12a)
Bp,t'yt(Ep,t1)Tyt (12b)

where A1, bt, Bp,t', and Ep,t1 are the coefficients regarding the constraints at the day-ahead operational stage.

The constraints of the second- and third-level “max-min” problem are integrated in (12c)-(12e). Formula (12c) includes all the linear constraints in the second and third levels (i.e., (4), (5a)-(5c), and (5e)-(5k) for the upper-level problem and (8) for the lower-level problem). Formula (12d) is composed of the second-stage SOC constraints, i.e., SOC reformulations of constraint (5d) with (10). Formula (12e) denotes the uncertainty interval ϒ of uncertainty variables related to (9).

A2xtFtyt+Gtut+dt (12c)
Bp,t''xt2(Ep,t2)Txt (12d)
uϒ (12e)

where A2, Ft, Gt, dt, Bp,t, and Ep,t2 are the coefficients regarding constraints at the real-time regulating stage.

The model in (11) and (12) respects a standard adaptive RO model with SOC constraints, which could be effectively solved by the C&CG algorithm [

32].

Specifically, the master problem (MP) is defined as:

minyttatTyt+φ (13a)

s.t.

φtgtTxt(k)k (13b)
A1ytbt (13c)
Bp,t'yt(Ep,t1)Tyt (13d)
A2xt(k)Ftyt+Gtut(k),*+dtk (13e)
Bp,t''xt(k)2(Ep,t2)Txt(k)k (13f)

where the superscript k denotes the C&CG iteration index.

The MP determines the optimal first-stage solution yt* in the identified worst-case scenarios ut(k),*. Accordingly, with the updated yt*, the primal subproblem (Primal-SP) is constructed as:

maxuminxttgtTxt (14a)

s.t.

A2xtFtyt*+Gtut+dt:πt (14b)
Bp,t''xt2(Ep,t2)Txt: (νp,t,μp,t) (14c)

where πt is the dual vector for inequality constraint (14b); and νp,t and μp,t are the dual vectors for SOC constraint (14c).

The two-level “max-min” Primal-SP can be re-formulated to the single-level dual subproblem (Dual-SP) with the aid of the strong duality theory [

31] as:

maxu,π,w,λt[(Ftyt*+bt)Tπt+utTGtπt+dtTπt] (15a)
A2Tπt+p(Bp,t''vp,t+Ep,t2μp,t)=gt (15b)
vp,t2μp,t (15c)

The bilinear term utTGtπt in (15a) causes the non-convexity of the Dual-SP, which could be easily linearized by the big-M method [

31]. By this means, the non-convex Dual-SP is hence transformed into a mixed-integer second-order cone programming (MISOCP) problem and could be resolved by the existing commercial solvers.

The optimal solution could be explored by iteratively solving the MP and subproblem using C&CG algorithm, as shown in Algorithm 1. It should be noted that since no SOC constraint is involved in the lower-level model, the SOC constraints are omitted when solving the lower-level model.

Algorithm 1  : C&CG algorithm

1:

Initialize the lower bound LB=-, the upper bound UB=+, the C&CG iteration index k=1, and C&CG convergence tolerance ψ1

2:

Solve the MP problem in (13) with given ut(k-1),*, and obtain yt(k-1),* and φ(k),*. Update LB=maxLB,tatTyt(k),*+φ(k),*

3:

Solve the Dual-SP in (15) with given yt(k),*, and obtain ut(k),* and πt(k),*. Update UB=minUB,t[(Ftyt(k),*+bt)Tπt(k),*+((ut(k),*)TGt+dtT)πt(k),*]

4:

If |(UB-LB)/LB|ψ1, terminate; otherwise, k=k+1, and go to Step 2

The upper-level IEGM will be cleared when the SOC-dual-based C&CG algorithm converges. Recalling the MP at the last iteration, inspired by [

16] and [18], the ULMEP and ULMGP can be derived from the Lagrangian function of (13a) L with (16a) and (16b), respectively.

βb,tE=LDb,tE=λb,tE+f(λf,tl--λf,tl+)Hf-b+          k𝒦λb,tΔ,(k)+f(λf,tl-,(k)-λf,tl+,(k))Hf-b (16a)
βn,tG=LDn,tG=μn,tG+k𝒦(μn,tΔ,(k)+μn,tG,(k)) (16b)

where λb,tE, λf,tl+/-, λf,tl+/-,(k), and λb,tΔ,(k) are the dual variables of constraints (2e), (2h), (4c), and (4f), respectively; and μn,tG, μn,tG,(k), and μn,tΔ,(k) are the dual variables of (3a), (5a), and (5j), respectively.

From (16a) and (16b), we can observe that unlike the traditional energy prices, the robust marginal prices can reflect the impact of uncertainties on the prices.

B. Coordinating Upper- and Lower-level Models Based on BRD Algorithm

The upper-level model determines the energy prices considering uncertainties, i.e., ULEMPs and ULGMPs. Given the energy prices, the lower-level model schedules the energy input of SEH for the upper-level model. To effectively reach an equilibrium of the upper-level and lower-level models and improve the convergence of the solving algorithm, a BRD algorithm [

16] is developed, as shown in Algorithm 2.

Algorithm 2  : BRD algorithm

1:

Initialize the energy inputs of SEHs Pe,tin,0 and Ge,tin,0 by setting the ULMEP and ULMGP as the prices of the cheapest CPP and gas well, respectively, and solve (6)-(8). Set the BRD iteration index o=1 and BRD convergence tolerance ψ2

2:

Solve the upper-level model for the IEGM clearing with (1)-(5) given Pe,tin,(o-1) and Ge,tin,(o-1), and obtain the ULMEP and ULMGP βb,tE,(o) and βn,tG,(o)

3:

Solve the lower-level model for the robust bidding of SEH with (6)-(8) given βb,tE,(o) and βn,tG,(o), and obtain the energy inputs of SEHs Pe,tin,(o) and Ge,tin,(o)

4:

If the maximal residual satisfies (17), terminate; otherwise, o=o+1, and go to Step 2

maxPe,tin,(o)-Pe,tin,(o-1)Pe,tin,(o),Ge,tin,(o)-Ge,tin,(o-1)Ge,tin,(o)ψ2    (17)

It should be emphasized that the BRD algorithm is heuristic. Accordingly, its convergence is sensitive to the initial point [

16]. In addithon, the iterative-based BRD algorithm is susceptible to falling into the oscillation mode, which could cause repeated iterations [33].

To address these limitations and enhance the convergence performance of the BRD algorithm, two effective measures are introduced in the proposed solving procedures.

Inspired by [

16], the modified energy price updating strategy is introduced to reduce the influence of the selection on the convergence process. In detail, between Step 2 and Step 3 of Algorithm 2, a weighted combination price of newly obtained prices at the current iteration and the previous ones are transferred to the lower-level model. Specifically, in this paper, the energy prices obtained at the current iteration (βb,tE,(o) and βn,tG,(o)) and those from the last two iterations (βb,tE,(o-1), βn,tG,(o-1), βb,tE,(o-2), and βn,tG,(o-2)) are individually weighted with 1/3 and summed to determine the energy price delivered to the lower-level model. This strategy can well lower the reliance of the convergence characteristic of the BRD algorithm on the selection of the initial points.

Secondly, to avoid the potential oscillation mode, enlightened by [

13], a bi-section strategy is incorporated. Supposing the electricity input Pe,tin,(o) oscillates at the oth iteration in Algorithm 2, set Pe,tin,max=max{Pe,tin,(o),Pe,tin,(o-1)} and Pe,tin,min=min{Pe,tin,(o),Pe,tin,(o-1)}.

Then, the following sub-steps are conducted to obtain the operational equilibrium.

Sub-step 1: embed the constraint Pe,tin,(o)=(Pe,tin,max+Pe,tin,min)/2 into the lower-level model for the robust bidding of SEH, and then solve it. If the convergence condition (17) is satisfied, terminate Algorithm 2. Otherwise, increase o by 1 and go to Sub-step 2.

Sub-step 2: insert constraint Pe,tin,minPe,tin,(o)Pe,tin,max into the lower-level model for the robust bidding of SEH and then solve it. If the convergence condition (17) is satisfied, terminate Algorithm 2. Otherwise, go to Sub-step 3.

Sub-step 3: if Pe,tin,(o)=Pe,tin,max, update the lower bound Pe,tin,min=(Pe,tin,max+Pe,tin,min)/2. Otherwise, if Pe,tin,(o)=Pe,tin,min, update the upper bound Pe,tin,max=(Pe,tin,max+Pe,tin,min)/2. Go to Sub-step 1.

By iteratively conducting Sub-steps 1-3, the operation interval of SEH will be repeatedly narrowed by half. Following this, the SEH can well converge to its optimal equilibrium operating point.

C. Wholistic Solution Process

The flow chart in Fig. 3 illustrates the solution process of the proposed model. A step-by-step solution procedure is given as follows.

Fig. 3  Flow chart for solution process of proposed model.

Step I  : initialization of the parameters and energy input of SEHs. Input the parameters of the IEGM and SEHs, and initialize the BRD iteration index and BRD convergence tolerance. Furthermore, calculate the initial electricity and gas inputs of each SEH by solving (6)-(9) with Algorithm 1.

Step II  : robust clearing mechanism of the IEGM. Given the energy inputs of SEHs, obtain the ULMEP and ULMGP by solving (1)-(5) with Algorithm 1 and (16), respectively.

Step III  : robust bidding of SEHs. Given the updated ULMEP and ULMGP, obtain and update the electricity and gas inputs of each SEH by solving (6)-(9) with Algorithm 1. Check the BRD convergence criterion (17). If it is satisfied, terminate; otherwise, go to Step II.

V. Numerical Experiments

A. Test System and Parameters

The proposed framework and its solving algorithm are tested on an integrated electricity and natural gas system coupled with SEHs. The testing system comprises a modified IEEE 39-bus electric system [

34] and a Belgian 20-node gas system [35]. The electric system includes 46 branches and 10 power plants (4 GPPs and 6 CPPs). The gas system consists of 19 pipelines, 6 gas wells, and 3 gas compressors. Ten gas loads are also included, of which 4 gas loads are connected to the GPPs. Four SEHs are linked to electric buses 9, 16, 20, and 29 while connected to the gas system through gas nodes 3, 10, 16, and 19, respectively. Additionally, the electric and gas loads of the testing system and the MED profile of SEHs can refer to [18] and [10], respectively. The uncertainty budgets Γ1 and Γ2 are preset at 4 and 1, respectively, to obtain the solutions to approximately cover 95% of the actual output of wind farms [36]. The uncertainty output interval of wind farms is ±10%. Numerical experiments are all conducted using the coding platform of MATLAB R2022 and the optimization models are solved by Gurobi 9.0.1.

B. Clearing Results of IEGM Considering Different Levels of Uncertainties

The uncertainty energy prices, i.e., the ULMEP and ULMGP, are presented in Fig. 4. The ULMEP differences among different electric buses are induced by transmission network limits. Both the transmission congestion in the base-case and worst-case scenarios will affect the ULMEP. From Hour 4 to Hour 6, the electric load is not heavy, which leads to a nearly consistent ULMEP in the electric power system. During other periods, the line congestion causes the ULMEPs at buses 34 and 36 to be relatively higher than those at other buses. Much differently, in the gas system, the ULMGP of each gas node almost remains unchanged over the entire dispatching horizon, which reflects the efficacy of considering the linepack in the gas system. However, the gas price at gas node 8 is higher, which is the discharge node of the gas compressor. Since the pipeline capacity limit is binding at the compressor branch, the gas load at node 8 can only be supplied by the expensive gas well 2.

Fig. 4  Uncertainty energy prices. (a) ULMEP. (b) ULMGP.

The energy demand difference and energy price difference of SEH 2 are depicted in Fig. 5. To facilitate the comparability, we have standardized the units to MMBtu for both electricity and gas inputs. Similarly, the ULMEP and ULMGP are both expressed in $/MMBtu. The energy demand difference is computed by subtracting the electricity demand from the gas demand in the SEH. Furthermore, the energy price difference is calculated by the electricity price minus the gas price. It can be observed that the SEH tends to consume a greater amount of gas and a reduced quantity of electricity when there is a significant price difference, e.g., during the interval from Hour 8 to Hour 12. In contrast, the energy demand gap is narrower during the period from Hour 3 to Hour 7 when the price gap is comparatively modest. This signifies that from an economic perspective, the SEH opts to utilize more electricity when the energy price gap is smaller. The result illustrates that the energy prices cleared in the IEGM could effectively guide the SEH to optimize its energy portfolio. Moreover, the mutual substitution feature of the SEH in fulfilling its terminal energy requirements is also demonstrated.

Fig. 5  Energy demand difference and energy price difference of SEH 2.

As illustrated in Section IV, the two-level robust clearing model is resolved by utilizing the BRD algorithm. The computation performance of this algorithm is conspicuously depicted in Fig. 6. The algorithm reaches convergence after 9 iterations, achieving this within a total elapsed time of 2556 s. The result thereby underscores its efficacy in addressing the clearing framework we have put forth.

Fig. 6  Computation performance of BRD algorithm.

C. Comparison Between Proposed Model and Centralized Model

The comparative results of SEH costs under different load factors are illustrated in Table I. We can observe that the obtained total costs of SEHs with the proposed model (two-level model) are approaching those with the centralized model under prescribed load factors. The results demonstrate that the proposed model will not largely compromise the optimality.

TABLE I  Comparative Results for SEH Costs Under Different Load Factors
ModelLoad factorTotal cost (105 $)Electricity cost (105 $)Gas cost (105 $)
Centralized model 0.7 8.402 6.414 1.988
1.0 13.678 11.120 2.558
1.3 19.233 16.038 3.195
Proposed model 0.7 8.419 6.425 1.994
1.0 13.694 11.129 2.565
1.3 19.250 16.047 3.203

In the base case of the load factor equal to 1, the computation time of the proposed model exhibits a certain degree of increase compared with what of the centralized model. Specifically, the proposed model consumes 2556 s to obtain the solution while the centralized model completes the computation within 1785 s. The longer computation time is primarily attributed to the iterative solving progress. Despite the iterations, the proposed model entails a lower-dimensional optimization model compared with the centralized model. In this case, the computation time of the proposed model does not significantly surpass that of the centralized one.

In conclusion, the proposed model can respect decision autonomy without significantly increasing computation time or largely compromising optimality. To conclude, the proposed model exhibits distinct advantages over the centralized one.

D. Impact of Uncertainty Degree on Result of Proposed Model

The deviations of utility-level uncertainties are varied to test the impact of uncertainty intervals on the ULMEPs and ULMGPs. The average ULMEPs and ULMGPs of system with different uncertainty degrees are depicted in Fig. 7. The results show that as uncertainty degrees intensify, the ULMEPs and ULMGPs both increase. This is because the incremental uncertaint induce the integrated electricity and gas system to deploy additional electrical reserves from GPPs and CPPs to manage the uncertainties in the worst case, consequently leading to the increase of ULMEPs. Since the rising output of GPPs contributes to the elevating gas demand, the ULMGPs are escalated as well.

Fig. 7  Average ULMEPs and ULMGPs of system with different uncertainty degrees. (a) ULMEPs. (b) ULMGPs.

E. Comparison with Model only Considering Single-level Uncertainties

To demonstrate the superiority of the proposed model for concurrently considering and modeling utility-level and distribution-level uncertainties, two cases are introduced.

1) Case 1 (C-1): the utility-level and distribution-level uncertainties are simultaneously handled with RO as in the proposed model.

2) Case 2 (C-2): only the utility-level uncertainties are considered, i.e., the deterministic optimization model is utilized for the SEH, which is similar to the model in [

18].

The comparative results regarding the costs of SEHs in C-1 and C-2 are displayed in Table II. Further, 100 random scenarios are generated with Monte Carlo simulation to mimic the actual realizations of utility-level and distribution-level uncertainties.

TABLE II  Costs of SEHs in C-1 and C-2
CaseFirst-stage cost (106 $)

Second-stage

cost (105 $)

CavgSEH (105 $)Sum of first-stage cost and CavgSEH (106 $)
C-1 1.0409 3.285 3.108 1.3517
C-2 1.0258 3.317 1.3575

In Table II, the first-stage cost and second-stage cost refer to costs of the day-ahead operational stage and the real-time regulating stage in (1), respectively. The averaged real-time regulation cost (ARRC) under actual realizations of distribution-level uncertainties CavgSEH is calculated as:

CavgSEH=esR=1SRρsR[CeELΔELsR,e,t(II)+CeGLΔGLsR,e,t(II)+CeHLΔHLsR,e,t(II)+CeW(PsR,e,tU,(II)-PsR,e,tD,(II))] (18)

where the subscript sR denotes the re-dispatch variables in the sRth scenarios of distribution uncertainties with the corresponding probability ρsR; and SR is the number of the actual realization scenarios. CavgSEH is calculated with (18) given the first-stage decisions as the boundary.

With the proposed model, the first-stage cost in C-1 is larger than that in C-2. The reason is that at the first stage, through the RO model, the SEH determines the optimal first-stage decisions with adequate reserves to prepare for the uncertainties that compromise the cost of the first stage (day-ahead stage) and the potential regulating cost with actual realizations of uncertainties. On the contrary, the deterministic model in C-2 merely gives the optimal energy use portfolio of SEH without considering the possible uncertainties. Along with this line, the ARRC in C-1 is less than that in C-2 since the SEH has dispatched adequate reserve resources to reduce the regulating cost with the realization of the uncertainties. The fifth column in Table II presents the sum of the first-stage cost and CavgSEH (hereinafter referred to as the summing cost). The lower summing cost in C-1 with respect to that in C-2 demonstrates the cost-effectiveness when considering distribution-level uncertainties with the RO model, which can lower the ARRC at the slight increase of first-stage cost, especially when the uncertainties are undesirably forecasted.

F. Validation of Proposed Model on IEEE 118-bus Electric System and 40-node Gas System

To validate the applicability of the proposed model on the larger-scale system, a large-scale system consisting of a 118-bus electric system and a 40-node gas system [

18] is further introduced. The electric system includes 186 branches, 12 GPPs, 32 CPPs, and 10 wind farms. The gas system is composed of 35 pipelines, 12 gas suppliers, and 4 gas compressors. In addition, 10 SEHs are connected to the system.

The overall computation time with the BRD and C&CG algorithms is 6348 s. The total BRD iteration reaches 13. Apparently, the relatively longer computation time is due to the growing scale of the IEGM, which leads to more time to solve the upper-level model. In addition, the increment of the number of lower-level robust bidding problems of SEHs also increases the overall computation time. Along this line, the resulting computation time is still acceptable for real-world day-ahead market implementation [

37].

Furthermore, the average electricity and gas prices of the electric and gas nodes, where SEHs are located, are depicted in Fig. 8(a). The summing electricity and gas inputs of all SEHs are displayed in Fig. 8(b). It can be observed that more electricity is input by SEHs when the electricity clearing prices are lower. As a contrary trend, the SEHs tend to consume more gas to satisfy their terminal energy use during the higher electricity price periods. The results further demonstrate the SEHs are capable of utilizing the proposed model to optimize the energy portfolio in a larger IEGM.

Fig. 8  Electricity and gas prices and inputs. (a) Average electricity and gas prices. (b) Summing electricity and gas inputs of all SEHs.

VI. Conclusion

To realize the accurate modeling of the IEGM clearing, this paper proposes a bi-level robust clearing framework of the IEGM to capture the interaction of different levels of uncertainties. In this framework, the upper-level model clears the IEGM, reflecting the influence of the utility-level uncertainties and its coupling with SEHs. The robust bidding model of the SEH at the lower level features the distribution-level uncertainties and determines the optimal energy input in response to the updated market prices. The BRD algorithm is combined with the C&CG algorithm to solve the robust clearing framework. Numerical results confirm that the proposed framework can elaborately model the influences of the intertwined uncertainties on the market clearing. Furthermore, the optimal biddings of SEHs considering the market clearing process can be also more accurately obtained.

However, in our study, the strategic bidding behaviors of the participants in the IEGM including SEHs are not considered. In reality, the market participants could potentially bid strategically to pursue the maximum economic benefits [

38]. Therefore, the robust clearing framework considering the participation of SEHs with the corresponding strategic bidding behaviors is worth exploiting in future research. Moreover, only energy trades between the IEGM and SEHs are considered. Since the operating reserve can be also traded in the market and can be supplied by the SEHs [39], [40], the robust joint energy-reserve clearing mechanism can also be exploited in future studies.

References

1

S. Wang, J. Zhai, H. Hui et al., “Operational reliability of integrated energy systems considering gas flow dynamics and demand-side flexibilities,” IEEE Transactions on Industrial Informatics, vol. 20, no. 2, pp. 1360-1373, Feb. 2024. [Baidu Scholar] 

2

S. Chen, A. J. Conejo, R. Sioshansi et al., “Equilibria in electricity and natural gas markets with strategic offers and bids,” IEEE Transactions on Power Systems, vol. 35, no. 3, pp. 1956-1966, May 2020. [Baidu Scholar] 

3

C. M. Correa-Posada and P. Sánchez-Martín, “Integrated power and natural gas model for energy adequacy in short-term operation,” IEEE Transactions on Power Systems, vol. 30, no. 6, pp. 3347-3355, Nov. 2015. [Baidu Scholar] 

4

A. Zlotnik, L. Roald, S. Backhaus et al., “Coordinated scheduling for interdependent electric power and natural gas infrastructures,” IEEE Transactions on Power Systems, vol. 32, no. 1, pp. 600-610, Jan. 2017. [Baidu Scholar] 

5

C. Wang, W. Wei, J. Wang et al., “Strategic offering and equilibrium in coupled gas and electricity markets,” IEEE Transactions on Power Systems, vol. 33, no. 1, pp. 290-306, Jan. 2018. [Baidu Scholar] 

6

M. Gil, P. Dueñas, and J. Reneses, “Electricity and natural gas interdependency: comparison of two methodologies for coupling large market models within the European regulatory framework,” IEEE Transactions on Power Systems, vol. 31, no. 1, pp. 361-369, Jan. 2016. [Baidu Scholar] 

7

C. Ordoudis, P. Pinson, and J. M. Morales, “An integrated market for electricity and natural gas systems with stochastic power producers,” European Journal of Operational Research, vol. 272, no. 2, pp. 642-654, Jan. 2019. [Baidu Scholar] 

8

C. Wang, A. R. Sayed, H. Zhang et al., “Two-stage distributionally robust strategic offering in pool-based coupled electricity and gas market,” Energy, vol. 265, p. 126288, Feb. 2023. [Baidu Scholar] 

9

S. Bahrami and A. Sheikhi, “From demand response in smart grid toward integrated demand response in smart energy hub,” IEEE Transactions on Smart Grid, vol. 7, no. 2, pp. 650-658, Mar. 2016. [Baidu Scholar] 

10

M. Bao, H. Hui, Y. Ding et al., “An efficient framework for exploiting operational flexibility of load energy hubs in risk management of integrated electricity-gas systems,” Applied Energy, vol. 338, p. 120765, May 2023. [Baidu Scholar] 

11

Y. Li, Z. Li, F. Wen et al., “Privacy-preserving optimal dispatch for an integrated power distribution and natural gas system in networked energy hubs,” IEEE Transactions on Sustainable Energy, vol. 10, no. 4, pp. 2028-2038, Oct. 2019. [Baidu Scholar] 

12

Z. Liang, M. Bao, Y. Ding et al., “A two-stage framework for exploiting the multi-elasticity of multi-energy demands in integrated electricity and gas markets,” IEEE Transactions on Power Systems, vol. 39, no. 2, pp. 3196-3210, Mar. 2024. [Baidu Scholar] 

13

Y. Cheng, N. Zhang, B. Zhang et al., “Low-carbon operation of multiple energy systems based on energy-carbon integrated prices,” IEEE Transactions on Smart Grid, vol. 11, no. 2, pp. 1307-1318, Mar. 2020. [Baidu Scholar] 

14

S. Yin and J. Wang, “Distributionally robust decentralized scheduling between the transmission market and local energy hubs,” IEEE Transactions on Power Systems, vol. 38, no. 2, pp. 1845-1856, Mar. 2023. [Baidu Scholar] 

15

R. Chen, J. Wang, and H. Sun, “Clearing and pricing for coordinated gas and electricity day-ahead markets considering wind power uncertainty,” IEEE Transactions on Power Systems, vol. 33, no. 3, pp. 2496-2508, May 2018. [Baidu Scholar] 

16

A. R. Sayed, C. Wang, W. Wei et al., “Robust operational equilibrium for electricity and gas markets considering bilateral energy and reserve contracts,” IEEE Transactions on Power Systems, vol. 36, no. 4, pp. 2891-2905, Jul. 2021. [Baidu Scholar] 

17

H. Ye, Y. Ge, M. Shahidehpour et al., “Uncertainty marginal price, transmission reserve, and day-ahead market clearing with robust unit commitment,” IEEE Transactions on Power Systems, vol. 32, no. 3, pp. 1782-1795, May 2017. [Baidu Scholar] 

18

F. Liu, X. Wang, Y. Xiao et al., “Robust pricing of energy and ancillary services in combined electricity and natural gas markets,” IEEE Transactions on Power Systems, vol. 37, no. 1, pp. 603-616, Jan. 2022. [Baidu Scholar] 

19

D. Rakipour and H. Barati, “Probabilistic optimization in operation of energy hub with participation of renewable energy resources and demand response,” Energy, vol. 173, pp. 384-399, Apr. 2019. [Baidu Scholar] 

20

V. Davatgaran, M. Saniei, and S. S. Mortazavi, “Optimal bidding strategy for an energy hub in energy market,” Energy, vol. 148, pp. 482-493, Apr. 2018. [Baidu Scholar] 

21

M. S. Javadi, M. Lotfi, A. E. Nezhad et al., “Optimal operation of energy hubs considering uncertainties and different time resolutions,” IEEE Transactions on Industry Applications, vol. 56, no. 5, pp. 5543-5552, Sept. 2020. [Baidu Scholar] 

22

Y. Wen, X. Qu, W. Li et al., “Synergistic operation of electricity and natural gas networks via ADMM,” IEEE Transactions on Smart Grid, vol. 9, no. 5, pp. 4555-4565, Sept. 2018. [Baidu Scholar] 

23

C. He, L. Wu, T. Liu et al., “Robust co-optimization scheduling of electricity and natural gas systems via ADMM,” IEEE Transactions on Sustainable Energy, vol. 8, no. 2, pp. 658-670, Apr. 2017. [Baidu Scholar] 

24

R. Li, W. Wei, S. Mei et al., “Participation of an energy hub in electricity and heat distribution markets: an MPEC approach,” IEEE Transactions on Smart Grid, vol. 10, no. 4, pp. 3641-3653, Jul. 2019. [Baidu Scholar] 

25

F. Liu, Z. Bie, and X. Wang, “Day-ahead dispatch of integrated electricity and natural gas system considering reserve scheduling and renewable uncertainties,” IEEE Transactions on Sustainable Energy, vol. 10, no. 2, pp. 646-658, Apr. 2019. [Baidu Scholar] 

26

B. Zhao, A. Zlotnik, A. J. Conejo et al., “Shadow price-based co-ordination of natural gas and electric power systems,” IEEE Transactions on Power Systems, vol. 34, no. 3, pp. 1942-1954, May 2019. [Baidu Scholar] 

27

X. Zheng, B. Zhou, X. Wang et al., “Day-ahead network-constrained unit commitment considering distributional robustness and intraday discreteness: a sparse solution approach,” Journal of Modern Power Systems and Clean Energy, vol. 11, no. 2, pp. 489-501, Mar. 2023. [Baidu Scholar] 

28

Y. Wang, K. Zhang, and K. Qu, “Segmented real-time dispatch model and stochastic robust optimization for power-gas integrated systems with wind power uncertainty,” Journal of Modern Power Systems and Clean Energy, vol. 11, no. 5, pp. 1480-1493, Sept. 2023. [Baidu Scholar] 

29

N. G. Cobos, J. M. Arroyo, N. Alguacil et al., “Robust energy and reserve scheduling considering bulk energy storage units and wind uncertainty,” IEEE Transactions on Power Systems, vol. 33, no. 5, pp. 5206-5216, Sept. 2018. [Baidu Scholar] 

30

A. T. Farsani, M. Joorabin, and S. S. Mortazavi, “Participation of energy hubs in the market clearing price-based energy market considering generation and consumption uncertainty modelling,” IET Renewable Power Generation, vol. 17, no. 8, pp. 2079-2100, Jun. 2023. [Baidu Scholar] 

31

Y. He, M. Shahidehpour, Z. Li et al., “Robust constrained operation of integrated electricity-natural gas system considering distributed natural gas storage,” IEEE Transactions on Sustainable Energy, vol. 9, no. 3, pp. 1061-1071, Jul. 2018. [Baidu Scholar] 

32

B. Zeng and L. Zhao, “Solving two-stage robust optimization problems using a column-and-constraint generation method,” Operations Research Letters, vol. 41, no. 5, pp. 457-461, Sept. 2013. [Baidu Scholar] 

33

B. Wang, N. Zhao, and F. Li, “Constraint-based interactive approach for equilibrium of interdependent gas and electricity markets,” Applied Energy, vol. 335, p. 120704, Apr. 2023. [Baidu Scholar] 

34

F. Liu, Z. Bie, S. Liu et al., “Day-ahead optimal dispatch for wind integrated power system considering zonal reserve requirements,” Applied Energy, vol. 188, pp. 399-408, Feb. 2017. [Baidu Scholar] 

35

D. de Wolf and Y. Smeers, “The gas transmission problem solved by an extension of the simplex algorithm,” Management Science, vol. 46, no. 11, pp. 1454-1465, Nov. 2000. [Baidu Scholar] 

36

A. R. Sayed, C. Wang, J. Zhao et al., “Distribution-level robust energy management of power systems considering bidirectional interactions with gas systems,” IEEE Transactions on Smart Grid, vol. 11, no. 3, pp. 2092-2105, May 2020. [Baidu Scholar] 

37

J. Almeida, F. Lezama, J. Soares et al., “Preliminary results of advanced heuristic optimization in the risk-based energy scheduling competition,” in Proceedings of the Genetic and Evolutionary Computation Conference Companion, New York, USA, pp. 1812-1816, Jul. 2022. [Baidu Scholar] 

38

S. Fan, Z. Li, J. Wang et al., “Cooperative economic scheduling for multiple energy hubs: a bargaining game theoretic perspective,” IEEE Access, vol. 6, pp. 27777-27789, May 2018. [Baidu Scholar] 

39

Y. Hou, M. Bao, M. Sang et al., “A market framework to exploit the multi-energy operating reserve of smart energy hubs in the integrated electricity-gas systems,” Applied Energy, vol. 357, p. 122279, Mar. 2024. [Baidu Scholar] 

40

B. Hu, Y. Sun, C. Shao et al., “A decentralized market framework for procurement of operating reserves from district energy systems,” IEEE Transactions on Sustainable Energy, vol. 12, no. 3, pp. 1629-1639, Jul. 2021. [Baidu Scholar]