# Two-Stage Participation of Grid-Connected Energy Hubs in Day-Ahead and Real-Time Energy Markets Based on Robust Modeling of Uncertainties

Authors: Mohammad Jafar Nafehfeshan; Bahman Bahmani-Firouzi; Mehdi Nafar
Journal: International Journal of Energy Research
 (2023)

---

## Abstract

Optimal management of sources, storage, and responsive loads in the energy network leads to improvement of economic and technical situation of the network. However, note that single-layer energy management requires to provide the network operator with more amounts of information. This results in the increase in the processing speed for the network operator. Also, coordination of sources, storage, and responsive loads in the form of an aggregator such as an energy hub can play an important role for them in promoting the economic situation and energy efficiency. So, this paper presents the participation of grid-connected energy hubs (EHs) in the day-ahead (DA) and real-time (RT) energy markets. Scheme has a two-stage optimization framework. In the first stage, participation of EHs in the hourly DA energy market is executed based on a two-layer energy management system that is modeled as bilevel problem. EH model in DA market is formulated in upper-level problem aimed at maximizing their profit subject to meeting the linearized constraints of the EH’s operation. Lower-level problem deals with the operation of energy networks to minimize cost of energy purchase from the upstream network constrained to linearized optimal power flow equations. The second stage of participation of EHs in a 5 min RT market has an objective function to minimize the unbalance between the EHs’ profits in the mentioned two markets. This stage is subject to the upper-level constraints and the lower-level formulation in the first stage. Robust optimization models the uncertainties of load, renewable power, energy price, and energy demand of mobile storage to achieve an optimal robust solution. Benders decomposition algorithm solves bilevel formulation in two-stage scheme. Numerical results confirm capabilities of the scheme in enhancing operation status of the networks and economic situation of EHs. Hence, the Benders decomposition algorithm is able to extract the optimal solution in shorter computation time. Two-layer energy management needs shorter computation time than the single-layer management. Simultaneous modeling of two real-time and day-ahead energy markets reduces the unbalance cost between these markets. The proposed scheme can enhance the operating and economic indicators of the network by up to around 18% and 26% compared to power flow studies of energy networks.


---

## Body

## 1. Introduction

### 1.1. Motivation

Nowadays, the exploitation of technologies with less emission levels has been introduced as an effort to reduce environmental pollutions caused by uncontrolled consumption of fossil fuels. Among these technologies are electric vehicles (EVs), renewable energy sources (RESs), high-efficiency energy generation systems such as combined heat and power (CHP), appropriate schemes for energy management such as energy storage systems (ESSs), and demand response programs (DRPs), which has gained the attention of energy consumers and various organizations [1]. It is predicted that the number of these sources and active loads (EVs, DRP, and ESS) will witness huge progress and utilization in the area of power systems, especially distribution networks. Since current distribution networks are difficult to control and manage due to a high number of such elements, it is recommended based on the smart distribution network theory to consider such sources and active loads as several integrated units in a microgrid (MG), virtual power plant (VPP), and energy hub (EH) designs [2]. Additionally, because of the interdependency of various energies, it has been advised to evaluate the management of several types of energies simultaneously. For instance, when modeling the CHP, it is required to consider separate models of all three electricity, gas, and heating parts of it [3]. Consequently, referring to the above statements, choosing an EH form as a coordinator of active loads and renewable and nonrenewable energy sources is advisable, which helps improve energy efficiency compared with single energy management of various energy types [4]. As EHs are connected to different energy networks such as electricity, gas, and heating networks in the form of a local system, it is anticipated that proper energy management of these networks in the presence of EHs could provide desirable capability in enhancing technical aspects of the aforementioned networks [5]. In the following, it is expected that EHs, considering their obtained capabilities from control and management of sources and active loads, can acquire optimal financial profit in the energy market for their elements.


### 1.2. Literature Review

Various pieces of research have been conducted in the field of energy management of EHs; some are elaborated below. A distributed multiperiod multienergy operational model is introduced in [6] for multicarrier energy systems, where the role of energy hubs is to make decisions and participate in the generation, distribution, and consumption of different energy systems. In an attempt to provide desirable energy conversion efficiency for several energy carriers in the energy market, Ref. [7] incorporates AC/DC MGs as EHs, as well as utilizing stochastic day-ahead (DA) bidding strategy for operators of EHs. The aim is to establish a balance between operating cost and quality of thermal service, despite being present in the DA and real-time (RT) energy markets. In another study [8], the operation of EHs is modeled using optimal probabilistic scheduling. In this research, EHs are scheduled in a way that energy carriers are purchased and stored so that energy demand is appropriately met and the overall cost of EHs is minimized. A new framework for EHs according to Lyapunov’s optimization-based energy management strategy is exploited in [9], where an energy router is also adopted together with analyzing the energy balance equations of EHs. Hybrid stochastic/information gap decision theory (IGDT) is also used in [10] to assess the scheduling of EHs with wind turbines, CHPs, and ESSs. The authors attempted to introduce a hybrid stochastic/IGDT optimization method to optimally schedule EHs consisting of wind turbines. A novel form of chance-constrained optimization approach is adopted in [11] to manage the operation of EHs that include electrical and thermal demands, where renewable energy is generated. Energy management of EHs is realized in [12] via EH collaboration in DA energy markets, where uncertainties such as energy price, power generation of RESs, and parameters of EVs are considered in the proposed formulation. EHs with electrical and thermal demands are scheduled in [13], where modeling of energy losses and costs together with feasible operation zone of CHPs are fully provided. Taking into account uncertainties of different types, smart residential EHs, known as SREHs, are optimally scheduled in [14], where the impact and consequences of energy market prices, energy demand volume, and irradiation as uncertainty parameters are considered in the scheduling of SREHs using two-stage stochastic planning constrained to risk. The authors in [15] model the management of multicarrier energy systems as an optimization problem with an objective function aiming to make the energy cost of EHs minimized.Reference [16] contributes to three different aspects. First, it proposes a novel metaheuristic approach for medical cost-efficient task schedules, where an intelligent scheduler manages the tasks, such as the rate of service schedule, and lists items utilized by users during the data processing and computation through the fog node. Second, the quality-of-service (QoS) efficient-computation algorithm, which effectively monitors performance according to the indicator (parameter) with the analysis mechanism of quality-of-experience (QoE), has been developed. Third, a framework of blockchain-distributed technology-enabled QoS (QoS-ledger) computation in healthcare applications is proposed in a permissionless public peer-to-peer (P2P) network, which stores medical processed information in a distributed ledger. In [17], it analyzes South Korea’s strategy in the global oil market. South Korean oil cooperation is characterized by the creation or termination of joint projects in the oil sector, as well as the Republic of Korea’s national project for the diversification of state-energy suppliers. In [18], it is dedicated to studying and modeling the interdependence between the oil returns and exchange-rate movements of oil-exporting and oil-importing countries. Globally, twelve countries/regions are investigated, representing more than 60% and 67% of all oil exports and imports. Reference [19] examines the impacts of Pakistani rupee volatility on monthly energy imports based on the nonlinear autoregressive distributed lag (NARDL) estimations. Reference [20] proposes the improvement of mathematical modeling apparatus at the stage of technological processes design and reconstruction. Reference [21] concerns the carbon emission and environment pollution and discussed involves considering energy, trade, society, and economy as the four bases of the topic to be addressed. Also, economic prosperity and energy development were the other topics presented in the paper that impact the emission problem. Price of energy and urbanization were assumed as the two factors impacting energy consumption [22], where the difference in incomes among different countries was taken into account. The Granger causality test approach and the impulse response function analysis were also discussed by dividing the set of countries into three categories: high-income, upper-income, and low-income countries. The data gathered by collecting the carbon emission produced in the transport sector in China during 14 years between 2005 and 2019 was adopted [23] to examine the impact of energy efficiency on carbon dioxide emission level. According to the results, as the quality of energy and its efficiency are enhanced, the level of dioxide decreases significantly. Reference [24] introduces a social index to analyze the impact of urbanization on energy generation and demand, economic situation, and the environment. Reference [25] adopts a novel decision-taking method to analyze the levelized cost imposed by some types of renewables. Another decision-taking approach is involved in [26] to discuss the benefits that can be achieved by utilizing clean energy. The authors in [27] focus on smart grids with considerable distributed energy generating units and adopts a fuzzy decision-making technique. Reference [28] provides a stochastic model for unit commitment of EHs containing various types of generating units, heat pumps, chillers, and some other elements. Another stochastic method was used in [29] to optimally operate an EH that consists of energy converters and storage units. The method acts as a biding strategy. Energy management of electrical and gas energies is handled in [30] so that the cost is reduced up to a great degree. The Eh used in the paper consists of different types of demand and supply units. A summary of the literature related to the present work is listed in Table 1.Table 1
Taxonomy of recent research works.
Ref.EMS modelEnergy market modelUncertainty modelDART[6]One layerNoNoStochastic[7]One layerYesNoStochastic[8]One layerNoNoStochastic[9]One layerNoNoStochastic[10]One layerNoNoHybrid stochastic/robust[11]One layerNoNoRobust[12]One layerYesNoStochastic[13]One layerNoNoStochastic[14]One layerYesNoStochastic[15]One layerYesNoDeterministic[16]—NoNoDeterministic[17]—NoNoDeterministic[18]—NoNoDeterministic[19]—NoNoDeterministic[20]—NoNoDeterministic[21]—NoNoDeterministic[22]—NoNoDeterministic[23]—NoNoDeterministic[24]—NoNoDeterministic[25]—NoNoDeterministic[26]—NoNoDeterministic[27]—NoNoDeterministic[28]One layerNoNoStochastic[29]One layerYesNoStochastic[30]One layerNoNoDeterministicProposed modelTwo-layerYesYesRobust


### 1.3. Research Gaps

As it is observed from the literature review and Table1, there are some major research gaps in the field of energy management of EHs as follows:
(i)
Most studies present an integrated model for energy management of energy networks and sources and active loads of the EH. To put it simply, all sources and loads available in the EH send their data to the energy network operators (ENOs) and receive commands. Thus, only ENOs are responsible for the energy management of all electricity, gas, and heating networks in the proposed study. However, this drastically increases the volume of input data of ENOs and lessens their operation speed. Hence, to compensate for this issue, sources and loads present in the EH should be coordinated with the EH operator (EHO) and the EHO is responsible for managing the elements of the EH. Furthermore, for suitable energy management of energy networks, EHOs and ENOs should enjoy bilateral coordination, helping to reduce the volume of data in ENOs compared to the previous case and thus increasing the operation speed of each operator(ii)
Most pieces of works on energy management of EHs discuss the DA energy market model, but it is worth noting that the presence of RESs and active loads, i.e., EVs, leads to unbalance between the DA and RT markets due to numerous uncertain parameters. To tackle this challenge, the EH model should be modeled in accordance with DA and RT markets in the energy management problem of EHs(iii)
A considerable number of researchers adopt stochastic modeling for uncertain parameters. Such modeling requires the identification of a suitable probability distribution function (PDF), which is a time-demanding task. For instance, the consumption load of different energy networks is based on normal PDF. To achieve this, the average and standard deviation values of the consumption load are needed. Since finding accurate figures demands statistical data of a long period (e.g., one year), achieving a proper PDF for the consumption load or any other uncertainty parameter is time-consuming. Moreover, the mentioned problem has a high number of uncertainty parameters; hence, reaching a reliable solution based on stochastic modeling needs a high number of scenarios, which leads to increase solution time of the problem. Note that scenario-based stochastic programming (SBSP) can reduce the number of scenarios. But as the number of uncertainty parameters in the energy management problem of EHs is high (due to load, energy price, RESs, EVs, etc.), the problem will have a high number of scenarios in its stochastic model. As a solution, it is recommended to use robust optimization for modeling uncertainty parameters, which has been discussed in few studies such as [11]. That requires optimization in one scenario (the worst-case scenario), so its execution speed is considerably high. Yet, it should be noted that the solution time is of notable importance in the operation of a system. For example, in RT energy market problems, market clearing is performed every 15 minutes in some cases. Consequently, solving the problem in such conditions is crucial. Accordingly, robust modeling of uncertainties plays a more effective role than stochastic modeling regarding the problem under study. Additionally, it can obtain a robust optimal solution for the proposed scheme in terms of its uncertainties


### 1.4. Contributions

Based on the research gaps presented above and in an attempt to fill the gaps, this paper introduces the optimal participation of EHs connected to electricity, gas, and heating networks in DA and RT markets. The proposed scheme is formulated in the framework of a two-stage optimization problem, where the participation of EHs in the hourly DA (5 min RT) energy market is considered in the first (second) stage. The EH in this scheme is a bilateral coordinator framework for RESs, CHP, stationary (battery) and mobile (EV) electrical energy storages, thermal energy storage (TES), boiler, and responsive loads in the DRP. Besides, to compensate for the first research gap, the first stage employs a two-layer energy management system (EMS) for the operation of the mentioned networks in the presence of EHs so that the bilateral coordination of sources and active loads with the EHO is applied in the first layer of the EMS, and the bilateral coordination between these operators and ENOs is considered in the second layer of the EMS. Accordingly, the first step will be based on a bilevel optimization. In the upper-level problem, the participation of EHs in the hourly DA energy market is formulated in accordance with maximizing their profit in this market constrained to the linearized model of operation of different sources and active loads. Furthermore, optimal operation of the mentioned energy networks in the presence of EHs is considered in the lower-level problem, where its objective function deals with minimizing the cost of energy purchase from the upstream network and its constraints include linearized optimal power flow (LOPF) equations. In the second stage, the formulation of participation of EHs in the 5 min RT energy market is presented, which considers the minimization of the unbalance between the EHs’ profit in the mentioned two markets so that the second research gap is filled. This problem is also constrained to the operation model of elements of the EH (constraints of upper-level problem) and operation formulation of energy markets (lower-level problem) at 5 min time steps, i.e., it includes a bilevel optimization. The robust optimization is employed for modeling the uncertainties of load, energy price, the active power of RESs, and consumption energy of EVs so that the third research gap is fulfilled. The proposed method relies on the bounded uncertainty-based robust optimization (BURO) technique. In the first and second stages, the problem modeling will be based on the Benders decomposition (BD) algorithm to achieve the optimal solution in the least possible time. Eventually, the contributions of the suggested scheme are summarized as follows:(i)
Formulation of participation of grid-connected EHs in DA and RT markets in the form of a two-stage optimization(ii)
Modeling the DA operation of energy networks with EHs in accordance with a two-layer EMS in the form of a bilevel optimization based on the BD algorithm(iii)
Extraction of a robust optimal solution considering uncertainties of load, energy price, RESs, and EVs for the participation of EHs in the energy marketHypothesis of the suggested scheme include the following:(i)
Simultaneous modeling of DA and RT energy markets reduces the operation unbalance of these markets(ii)
Two-layer energy management is able to extract the optimal situation of EHs and energy networks(iii)
Optimal management of EHs can improve the operating and economic situation of energy networks(iv)
Robust optimization achieves a robust solution to the suggested scheme(v)
Benders decomposition algorithm can derive an optimal solution in less computation time


### 1.5. Paper Organization

The formulation of the proposed scheme is presented in Section2. The problem solution procedure is stated in Section 3. Finally, numerical results and conclusions are given in Sections 4 and 5, respectively.


## 1.1. Motivation

Nowadays, the exploitation of technologies with less emission levels has been introduced as an effort to reduce environmental pollutions caused by uncontrolled consumption of fossil fuels. Among these technologies are electric vehicles (EVs), renewable energy sources (RESs), high-efficiency energy generation systems such as combined heat and power (CHP), appropriate schemes for energy management such as energy storage systems (ESSs), and demand response programs (DRPs), which has gained the attention of energy consumers and various organizations [1]. It is predicted that the number of these sources and active loads (EVs, DRP, and ESS) will witness huge progress and utilization in the area of power systems, especially distribution networks. Since current distribution networks are difficult to control and manage due to a high number of such elements, it is recommended based on the smart distribution network theory to consider such sources and active loads as several integrated units in a microgrid (MG), virtual power plant (VPP), and energy hub (EH) designs [2]. Additionally, because of the interdependency of various energies, it has been advised to evaluate the management of several types of energies simultaneously. For instance, when modeling the CHP, it is required to consider separate models of all three electricity, gas, and heating parts of it [3]. Consequently, referring to the above statements, choosing an EH form as a coordinator of active loads and renewable and nonrenewable energy sources is advisable, which helps improve energy efficiency compared with single energy management of various energy types [4]. As EHs are connected to different energy networks such as electricity, gas, and heating networks in the form of a local system, it is anticipated that proper energy management of these networks in the presence of EHs could provide desirable capability in enhancing technical aspects of the aforementioned networks [5]. In the following, it is expected that EHs, considering their obtained capabilities from control and management of sources and active loads, can acquire optimal financial profit in the energy market for their elements.


## 1.2. Literature Review

Various pieces of research have been conducted in the field of energy management of EHs; some are elaborated below. A distributed multiperiod multienergy operational model is introduced in [6] for multicarrier energy systems, where the role of energy hubs is to make decisions and participate in the generation, distribution, and consumption of different energy systems. In an attempt to provide desirable energy conversion efficiency for several energy carriers in the energy market, Ref. [7] incorporates AC/DC MGs as EHs, as well as utilizing stochastic day-ahead (DA) bidding strategy for operators of EHs. The aim is to establish a balance between operating cost and quality of thermal service, despite being present in the DA and real-time (RT) energy markets. In another study [8], the operation of EHs is modeled using optimal probabilistic scheduling. In this research, EHs are scheduled in a way that energy carriers are purchased and stored so that energy demand is appropriately met and the overall cost of EHs is minimized. A new framework for EHs according to Lyapunov’s optimization-based energy management strategy is exploited in [9], where an energy router is also adopted together with analyzing the energy balance equations of EHs. Hybrid stochastic/information gap decision theory (IGDT) is also used in [10] to assess the scheduling of EHs with wind turbines, CHPs, and ESSs. The authors attempted to introduce a hybrid stochastic/IGDT optimization method to optimally schedule EHs consisting of wind turbines. A novel form of chance-constrained optimization approach is adopted in [11] to manage the operation of EHs that include electrical and thermal demands, where renewable energy is generated. Energy management of EHs is realized in [12] via EH collaboration in DA energy markets, where uncertainties such as energy price, power generation of RESs, and parameters of EVs are considered in the proposed formulation. EHs with electrical and thermal demands are scheduled in [13], where modeling of energy losses and costs together with feasible operation zone of CHPs are fully provided. Taking into account uncertainties of different types, smart residential EHs, known as SREHs, are optimally scheduled in [14], where the impact and consequences of energy market prices, energy demand volume, and irradiation as uncertainty parameters are considered in the scheduling of SREHs using two-stage stochastic planning constrained to risk. The authors in [15] model the management of multicarrier energy systems as an optimization problem with an objective function aiming to make the energy cost of EHs minimized.Reference [16] contributes to three different aspects. First, it proposes a novel metaheuristic approach for medical cost-efficient task schedules, where an intelligent scheduler manages the tasks, such as the rate of service schedule, and lists items utilized by users during the data processing and computation through the fog node. Second, the quality-of-service (QoS) efficient-computation algorithm, which effectively monitors performance according to the indicator (parameter) with the analysis mechanism of quality-of-experience (QoE), has been developed. Third, a framework of blockchain-distributed technology-enabled QoS (QoS-ledger) computation in healthcare applications is proposed in a permissionless public peer-to-peer (P2P) network, which stores medical processed information in a distributed ledger. In [17], it analyzes South Korea’s strategy in the global oil market. South Korean oil cooperation is characterized by the creation or termination of joint projects in the oil sector, as well as the Republic of Korea’s national project for the diversification of state-energy suppliers. In [18], it is dedicated to studying and modeling the interdependence between the oil returns and exchange-rate movements of oil-exporting and oil-importing countries. Globally, twelve countries/regions are investigated, representing more than 60% and 67% of all oil exports and imports. Reference [19] examines the impacts of Pakistani rupee volatility on monthly energy imports based on the nonlinear autoregressive distributed lag (NARDL) estimations. Reference [20] proposes the improvement of mathematical modeling apparatus at the stage of technological processes design and reconstruction. Reference [21] concerns the carbon emission and environment pollution and discussed involves considering energy, trade, society, and economy as the four bases of the topic to be addressed. Also, economic prosperity and energy development were the other topics presented in the paper that impact the emission problem. Price of energy and urbanization were assumed as the two factors impacting energy consumption [22], where the difference in incomes among different countries was taken into account. The Granger causality test approach and the impulse response function analysis were also discussed by dividing the set of countries into three categories: high-income, upper-income, and low-income countries. The data gathered by collecting the carbon emission produced in the transport sector in China during 14 years between 2005 and 2019 was adopted [23] to examine the impact of energy efficiency on carbon dioxide emission level. According to the results, as the quality of energy and its efficiency are enhanced, the level of dioxide decreases significantly. Reference [24] introduces a social index to analyze the impact of urbanization on energy generation and demand, economic situation, and the environment. Reference [25] adopts a novel decision-taking method to analyze the levelized cost imposed by some types of renewables. Another decision-taking approach is involved in [26] to discuss the benefits that can be achieved by utilizing clean energy. The authors in [27] focus on smart grids with considerable distributed energy generating units and adopts a fuzzy decision-making technique. Reference [28] provides a stochastic model for unit commitment of EHs containing various types of generating units, heat pumps, chillers, and some other elements. Another stochastic method was used in [29] to optimally operate an EH that consists of energy converters and storage units. The method acts as a biding strategy. Energy management of electrical and gas energies is handled in [30] so that the cost is reduced up to a great degree. The Eh used in the paper consists of different types of demand and supply units. A summary of the literature related to the present work is listed in Table 1.Table 1
Taxonomy of recent research works.
Ref.EMS modelEnergy market modelUncertainty modelDART[6]One layerNoNoStochastic[7]One layerYesNoStochastic[8]One layerNoNoStochastic[9]One layerNoNoStochastic[10]One layerNoNoHybrid stochastic/robust[11]One layerNoNoRobust[12]One layerYesNoStochastic[13]One layerNoNoStochastic[14]One layerYesNoStochastic[15]One layerYesNoDeterministic[16]—NoNoDeterministic[17]—NoNoDeterministic[18]—NoNoDeterministic[19]—NoNoDeterministic[20]—NoNoDeterministic[21]—NoNoDeterministic[22]—NoNoDeterministic[23]—NoNoDeterministic[24]—NoNoDeterministic[25]—NoNoDeterministic[26]—NoNoDeterministic[27]—NoNoDeterministic[28]One layerNoNoStochastic[29]One layerYesNoStochastic[30]One layerNoNoDeterministicProposed modelTwo-layerYesYesRobust


## 1.3. Research Gaps

As it is observed from the literature review and Table1, there are some major research gaps in the field of energy management of EHs as follows:
(i)
Most studies present an integrated model for energy management of energy networks and sources and active loads of the EH. To put it simply, all sources and loads available in the EH send their data to the energy network operators (ENOs) and receive commands. Thus, only ENOs are responsible for the energy management of all electricity, gas, and heating networks in the proposed study. However, this drastically increases the volume of input data of ENOs and lessens their operation speed. Hence, to compensate for this issue, sources and loads present in the EH should be coordinated with the EH operator (EHO) and the EHO is responsible for managing the elements of the EH. Furthermore, for suitable energy management of energy networks, EHOs and ENOs should enjoy bilateral coordination, helping to reduce the volume of data in ENOs compared to the previous case and thus increasing the operation speed of each operator(ii)
Most pieces of works on energy management of EHs discuss the DA energy market model, but it is worth noting that the presence of RESs and active loads, i.e., EVs, leads to unbalance between the DA and RT markets due to numerous uncertain parameters. To tackle this challenge, the EH model should be modeled in accordance with DA and RT markets in the energy management problem of EHs(iii)
A considerable number of researchers adopt stochastic modeling for uncertain parameters. Such modeling requires the identification of a suitable probability distribution function (PDF), which is a time-demanding task. For instance, the consumption load of different energy networks is based on normal PDF. To achieve this, the average and standard deviation values of the consumption load are needed. Since finding accurate figures demands statistical data of a long period (e.g., one year), achieving a proper PDF for the consumption load or any other uncertainty parameter is time-consuming. Moreover, the mentioned problem has a high number of uncertainty parameters; hence, reaching a reliable solution based on stochastic modeling needs a high number of scenarios, which leads to increase solution time of the problem. Note that scenario-based stochastic programming (SBSP) can reduce the number of scenarios. But as the number of uncertainty parameters in the energy management problem of EHs is high (due to load, energy price, RESs, EVs, etc.), the problem will have a high number of scenarios in its stochastic model. As a solution, it is recommended to use robust optimization for modeling uncertainty parameters, which has been discussed in few studies such as [11]. That requires optimization in one scenario (the worst-case scenario), so its execution speed is considerably high. Yet, it should be noted that the solution time is of notable importance in the operation of a system. For example, in RT energy market problems, market clearing is performed every 15 minutes in some cases. Consequently, solving the problem in such conditions is crucial. Accordingly, robust modeling of uncertainties plays a more effective role than stochastic modeling regarding the problem under study. Additionally, it can obtain a robust optimal solution for the proposed scheme in terms of its uncertainties


## 1.4. Contributions

Based on the research gaps presented above and in an attempt to fill the gaps, this paper introduces the optimal participation of EHs connected to electricity, gas, and heating networks in DA and RT markets. The proposed scheme is formulated in the framework of a two-stage optimization problem, where the participation of EHs in the hourly DA (5 min RT) energy market is considered in the first (second) stage. The EH in this scheme is a bilateral coordinator framework for RESs, CHP, stationary (battery) and mobile (EV) electrical energy storages, thermal energy storage (TES), boiler, and responsive loads in the DRP. Besides, to compensate for the first research gap, the first stage employs a two-layer energy management system (EMS) for the operation of the mentioned networks in the presence of EHs so that the bilateral coordination of sources and active loads with the EHO is applied in the first layer of the EMS, and the bilateral coordination between these operators and ENOs is considered in the second layer of the EMS. Accordingly, the first step will be based on a bilevel optimization. In the upper-level problem, the participation of EHs in the hourly DA energy market is formulated in accordance with maximizing their profit in this market constrained to the linearized model of operation of different sources and active loads. Furthermore, optimal operation of the mentioned energy networks in the presence of EHs is considered in the lower-level problem, where its objective function deals with minimizing the cost of energy purchase from the upstream network and its constraints include linearized optimal power flow (LOPF) equations. In the second stage, the formulation of participation of EHs in the 5 min RT energy market is presented, which considers the minimization of the unbalance between the EHs’ profit in the mentioned two markets so that the second research gap is filled. This problem is also constrained to the operation model of elements of the EH (constraints of upper-level problem) and operation formulation of energy markets (lower-level problem) at 5 min time steps, i.e., it includes a bilevel optimization. The robust optimization is employed for modeling the uncertainties of load, energy price, the active power of RESs, and consumption energy of EVs so that the third research gap is fulfilled. The proposed method relies on the bounded uncertainty-based robust optimization (BURO) technique. In the first and second stages, the problem modeling will be based on the Benders decomposition (BD) algorithm to achieve the optimal solution in the least possible time. Eventually, the contributions of the suggested scheme are summarized as follows:(i)
Formulation of participation of grid-connected EHs in DA and RT markets in the form of a two-stage optimization(ii)
Modeling the DA operation of energy networks with EHs in accordance with a two-layer EMS in the form of a bilevel optimization based on the BD algorithm(iii)
Extraction of a robust optimal solution considering uncertainties of load, energy price, RESs, and EVs for the participation of EHs in the energy marketHypothesis of the suggested scheme include the following:(i)
Simultaneous modeling of DA and RT energy markets reduces the operation unbalance of these markets(ii)
Two-layer energy management is able to extract the optimal situation of EHs and energy networks(iii)
Optimal management of EHs can improve the operating and economic situation of energy networks(iv)
Robust optimization achieves a robust solution to the suggested scheme(v)
Benders decomposition algorithm can derive an optimal solution in less computation time


## 1.5. Paper Organization

The formulation of the proposed scheme is presented in Section2. The problem solution procedure is stated in Section 3. Finally, numerical results and conclusions are given in Sections 4 and 5, respectively.


## 2. Model of the Proposed Problem

The proposed scheme is an optimization problem. The optimization formulation includes objective function [31–34] and constraints [35–37]. Also, the optimization can be implemented on the energy network, if there are smart [38–43] and telecommunication [44–50] platform in this system.This section addresses the two-stage scheduling management of EHs in electricity, gas, and heating networks based on a two-layer EMS [51] in the deterministic and robust model, as shown in Figure 1. The hourly DA operation of EHs in the mentioned networks is presented in the first stage that has a bilevel optimization, where the upper- (lower-) level problem is according to the first (second) layer of the EMS. The first layer of EMS focuses on bilateral coordination of sources, ESSs, responsive loads, and the EH operator. Thus, the upper-level problem in the first stage deals with the optimal participation of EHs in the DA energy market, aimed at maximizing the profit of EHs in the market subject to meeting the operational constraints of EHs elements. The second layer of the EMS is related to the bilateral coordination of EH operators and the operator of mentioned energy networks. Hence, the lower-level problem in the first stage puts forward the modeling of hourly operation of networks with EHs, where the minimization of cost of energy purchase by the mentioned networks from the upstream networks constrained to LOPF equations is formulated. The second stage considers the 5 min scheduling of EHs in the RT energy market, which has a model that is same with first stage problem, but the objective function of upper-level formulation in this stage minimizes the EHs’ unbalanced profit in both DA and RT markets. Moreover, robust optimization is incorporated in this section to model the uncertainties of the suggested scheme; the details of which are presented in the remaining sections.Figure 1
Two-stage scheme of EH scheduling in energy networks based on two-layer EMS.


### 2.1. DA Operation of EHs (First-Layer EMS and Upper-Level Problem Formulation)

This subsection addresses the optimal participation of EHs in the hourly DA energy market. In this problem, as given in (1a), maximization of the profit of EHs in the mentioned market is taken into account as the objective function. As per (1b), the profit of EHs is equal to the sum of their revenues in electricity, gas, and heating markets. In each of these markets, the revenue equals the energy price multiplied by power [12]. Note that if the value of power is positive, (1b) will represent the revenue of the EH in the mentioned market, but if it is negative, (1b) will denote the cost of EHs in this market. The balances between active, reactive, heating, and gas powers between the generators and consumers in the EH are modeled in constraints (1c)–(1f) [52]. Also, (1g)–(1w) represent the model of operation of sources, storages, and responsive loads in the EH. Constraints (1g)–(1j) give the formulation of CHP’s operation [52, 53], where constraints (1g) and (1h) calculate the heating and gas power of the CHP in terms of its active power. The limits on the CHP output in both electrical and thermal sections are considered in (1i) and (1j). Constraint (1k) calculates the input gas power of the boiler in terms of its output thermal power, while the limit on its output capacity is modeled (1l). The operation model of RESs is given in (1m) and (1n), where the generated active power of these sources is equal to the maximum capacity of the sources in accordance with weather conditions and this is due to low operation cost and environmental pollution as stated by (1m) [54]. The limit on simultaneous control of active and reactive power of these sources follows (1n), where it addresses the limit on apparent power of these sources [54]. The formulation of operation of EESs is expressed in (1o)–(1r) [55], which represent the charging and discharging limits, limit on the energy stored in the ESS, and limit on the charger capacity (apparent power) of this ESS, respectively. It is worth noting that the operation model of all EVs available in the parking lot is shown in constraints (1o)–(1r) [56, 57], except that parameters CRmax, DRmax, EEES0, and S¯EES change over time [58, 59]. Thus, the subscript h is used in the model of EV parking lot for these parameters. Additionally, (1q) is stated in the form of two constraints 0≤Ei,hEES0+∑τ=1hηE,ChPi,τCh−1/ηE,DisPi,τDis and ∑h∈ΞHEi,hEES0+ηE,ChPi,hCh−1/ηE,DisPi,hDis=E¯iEES in EV model, which refer to the limit on the energy stored in EVs and their energy demand in the operation horizon. For EV parking lot, parameters CRmax, DRmax, and S¯EES at hour h are found from CRmax=∑ev=1Nhcrev, DRmax=∑ev=1Nhdrev, and S¯EES=∑ev=1NhS¯evEV, where cr, dr, S¯EV, and Nh represent the charging rate, discharging rate, the capacity of EV charger, and the number of EVs connected to the EH during hour h [56]. Also, EEES0 at hour h is equal to the sum of total arrival (initial) energy of EVs that are recently connected to the parking lot in the EH. The operation model of TES is represented in (1s)–(1u) [52, 53], which has the same model as the EES does, except for thermal power variables of charging and discharging are used in these constraints and the limit on the charger capacity, (1r), cannot be used in the operation model of TES. Finally, the operation model of responsive loads in the electrical, gas, and thermal sections are given by constraints (1v) and (1w) [1]. In this type of DRP, consumers act according to (1a) corresponding to the price signal so that the profit of EHs in the energy market is maximized. The consumers participated in the DRP receive energy from the EH during low energy price, but they inject energy into the EH during high energy prices. As inexpensive (expensive) energy price is proportional to the off-peak (peak) hours, the proposed DRP can reduce consumers’ energy demand during peak intervals and supply this reduced amount of energy during off-peak hours [1]. Hence, constraint (1v) refers to the limit on the active, gas, and thermal power control of the DRP. Also, (1w) ensures that electrical (gas or thermal) energy decreased during peak hours is supplied during off-peak hours.Constraints (1i), (1n), and (1r) denote the linearized limit on the apparent power of CHP, RES, and ESS. It should be noted that these constraints are limited to a circular plane with the center on the origin and a radius of S, P2+Q2≤S2, and has a nonlinear format [60]. Yet, this limit can be approximated by a regular polygon with a model given by (1i), (1n), and (1r) [60]. Additionally, in the case the number of sides for this plane is high, the computational error will be negligible compared to that of the circular plane.
(1a)maxF1=∑h∈ΞHPRhEHssubject to(1b)PRhEHs=∑i∈ΞEHγhEPi,hEH+γhGGi,hEH+γhHHi,hEH∀h(1c)Pi,hEH=Pi,hCHP+Pi,hRES+Pi,hDR+Pi,hDis−Pi,hCh−Pi,hD∀i,h(1d)Qi,hEH=Qi,hCHP+Qi,hRES+Qi,hEES−Pi,hD∀i,h(1e)Hi,hEH=Hi,hCHP+Hi,hBO+Hi,hDR+Hi,hDis−Hi,hCh−Hi,hD∀i,h(1f)Gi,hEH=Gi,hDR−Gi,hCHP−Gi,hBO−Gi,hD∀i,h(1g)Hi,hCHP=1−ηt−ηlηhηtPi,hCHP∀i,h(1h)Gi,hCHP=1ηtPi,hCHP∀i,h(1i)Pi,hCHPcosm·Δβ+Qi,hCHPsinm·Δβ≤S¯iCHP∀i,h,m∈ΞM=1,2,⋯,nM,Δβ=2πnM(1j)0≤Hi,hCHP≤H¯iCHP∀i,h(1k)Gi,hBO=1ηbPi,hBO∀i,h(1l)0≤Hi,hBO≤H¯iBO∀i,h(1m)Pi,hRES=P¯iRES∀i,h(1n)Pi,hREScosm·Δβ+Qi,hRESsinm·Δβ≤S¯iRES∀i,h,m(1o)0≤Pi,hCh≤CRiEESxi,h∀i,h(1p)0≤Pi,hDis≤DRiEES1−xi,h∀i,h(1q)E¯iEES≤EiEES0+∑τ=1hηE,ChPi,τCh−1ηE,DisPi,τDis≤E¯iEES∀i,h(1r)Pi,hDis−Pi,hChcosm·Δβ+Qi,hEESsinm·Δβ≤S¯iEES∀i,h,m(1s)0≤Hi,hCh≤CRiTESyi,h∀i,h(1t)0≤Hi,hDis≤DRiTES1−yi,h∀i,h(1u)E¯iTES≤EiTES0+∑τ=1hηT,ChHi,τCh−1ηT,DisHi,τDis≤E¯iTES∀i,h(1v)−ξiPi,hDorHi,hDorGi,hD≤Pi,hDRorHi,hDRorGi,hDR≤ξiPi,hDorHi,hDorGi,hD∀i,h(1w)∑h∈ΞHPi,hDRorHi,hDRorGi,hDR=0∀i


### 2.2. DA Operation of Energy Networks (Second-Layer EMS and Lower-Level Formulation)

This subsection presents the hourly DA operation of electricity, gas, and heating networks in the presence of EHs. The problem is given in (2a)–(2u) and has an objective function aimed at minimizing the cost of energy purchased from the upstream electrical (first part), gas (second part), and heating (third section) networks [15]. Constraints of this problem as given by (2b)–(2u) refer to LOPF equations in the mentioned networks [12, 15]. Power flow equations are given in (2b)–(2l), where constraints (2b)–(2e) introduce power flow of the electricity network and refer to active and reactive power balance at the buses of the mentioned network and active and reactive power flowing through the line. Equations (2d) and (2e) are nonlinear in the forms of Pe,jEL=Ge,jVe2−VeVjGe,jcosφe−φj+Be,jsinφe−φj and Qe,jEL=−Be,jVe2+VeVjBe,jcosφe−φj−Ge,jsinφe−φj [60–64]. However, since the voltage phase angle between any two buses, φe−φj, in the distribution is less than 6° or 0.105 rad, the terms cosφe−φj and sinφe−φj are approximated by 1 and φe−φj, respectively. Moreover, since the magnitude of bus voltage should be between 0.9 and 1.1 p.u., it will be always close to 1 p.u. Hence, the voltage magnitude, V, is approximated as 1+ΔV, where ΔV represents the voltage deviation. According to these hypotheses and the small values of ΔV2 and ΔVφe−φj because of the small values of φe−φj and ΔV, constraints (2d) and (2e) replace the mentioned nonlinear constraints so that the model of LOPF is achieved [12]. Relations (2f)–(2j) concern power flow equations of the gas network, where the balance of gas power at different nodes is modeled in (2f). Also, (2g) calculates the gas power flow through the pipeline. In this equation, the squared power of the auxiliary variable σ, that is σ2, is in the form of signπg,h,πj,hπg,h2−πj,h2 [12, 15]. Hence, (2h)–(2j) present the linearized model of this equation based on the conventional piecewise linearization technique [60], where (2h) is the mentioned equation and (2i) and (2j) show the variables σ and π in terms of deviation variables of σ and π that are Δσ and Δπ, respectively. The formulation of power flow of the heating network is stated in (2k) and (2l), which, respectively, denote the balance of thermal power at different nodes and thermal power flowing through the pipeline.Constraints (2m)–(2u) present the operation limits of the different energy network [12, 15] and, respectively, refer to the limit on the voltage deviation of buses, limit on the capacity of the line (apparent power) and distribution substation, limit on the gas pressure at different nodes, limit on the capacity of pipeline (gas power) and gas station, limit on the temperature of nodes, and limit on the capacity (thermal power) of the pipeline and heating station. Note that, in the power flow model of the electricity network, given by (2b)–(2e), a voltage deviation variable is used to reach the LOPF; consequently, the limit on the voltage deviation given by (2m) substitutes the limit on voltage magnitude, V¯≤ΔVe≤V¯. Besides, the status of constraints (2n) and (2o) is according to the rule used in (1i).
(2a)minF2=∑h∈ΞHρhEPref,hES+ρhGGref,hES+ρhHHref,hESsubject to(2b)Pe,hES+∑i∈ΞEHIe,iE,EHPi,hEH−∑j∈ΞENIe,jELPe,j,hEL=Pe,hD∀e,h(2c)Qe,hES+∑i∈ΞEHIe,iE,EHQi,hEH−∑j∈ΞENIe,jELQe,j,hEL=Qe,hD∀e,h(2d)Pe,j,hEL=Ge,jΔVe,h−ΔVj,h−Be,jφe,h−φj,h∀e,j,h(2e)Qe,j,hEL=−Be,jΔVe,h−ΔVj,h−Ge,jφe,h−φj,h∀e,j,h(2f)Gg,hES+∑i∈ΞEHIg,iG,EHGi,hEH−∑j∈ΞGNIg,jGLGg,j,hGL=Gg,hD∀g,h(2g)Gg,j,hGL=κg,jsignπg,h,πj,hσg,j,h∀g,j,h(2h)Δσg,j,h,l=κg,jsignπg,h,πj,hSlπSlσΔπg,h,l−Δπj,h,l∀g,j,h,l(2i)σg,j,h=∑l∈ΞLΔσg,j,h,l∀g,j,h(2j)πg,h=π¯+∑l∈ΞLΔπg,h,l∀g,h(2k)Ht,hES+∑i∈ΞEHIt,iH,EHHi,hEH−∑j∈ΞHNIt,jHLHt,j,hHL=Ht,hD∀t,h(2l)Ht,j,hHL=ct,jṁt,jTt,h−Tj,h∀t,h(2m)1−V¯≤ΔVe,h≤V¯−1∀e,h(2n)Pe,j,hELcosm·Δβ+Qe,j,hELsinm·Δβ≤S¯e,jEL∀e,j,h,m(2o)Pe,hEScosm·Δβ+Qe,hESsinm·Δβ≤S¯eES∀e=ref,h,m(2p)π¯≤πg,h≤π¯∀g,h(2q)−G¯g,jGL≤Gg,j,hGL≤G¯g,jGL∀g,j,h(2r)−G¯gGS≤Gg,hGS≤G¯gGS∀g=ref,h(2s)T¯≤Tt,h≤T¯∀t,h(2t)−H¯t,jHL≤Ht,j,hHL≤H¯t,jHL∀t,j,h(2u)−H¯tHS≤Ht,hHS≤H¯tHS∀t=ref,h


### 2.3. Five-Minute RT Scheduling of EHs in Energy Networks

This subsection presents the optimal participation of EHs in the 5 min RT energy market subject to the condition that the unbalance between RT and DA operations is the minimum. It is worth noting that EHs consist of RESs. Since the active power generation by these sources depends on weather conditions and it is a challenging task with generally forecast errors, the operation results are expected to differ due to the uncertainty in the active power generation. This can lead to unbalance between supply and demand in the RT operation, and this is considered for low flexibility of the system [65, 66]. To compensate for this issue, the exploitation of flexibility sources, storage and DRP, has been recommended by many scholars [65]. Yet, this situation can be extracted to establish a proper optimization framework; so, this subsection, as given by (3a), the minimization of unbalance between the EHs’ profit in DA and RT energy markets is considered as the objective function [66]. The EHs’ profit in the RT market can be found from (3b). Also, the constraints of this problem as given in (3c) include the operation model of EHs, constraints of problem (1a)–(1w), and operation model of energy networks (problem 2). Note that the objective function in this problem is quadratic and its constraints have a linearized formulated framework. According to [52, 67], this model can be solved using linear programming (LP) solvers such as CPLEX in the least possible computation time.
(3a)minF3=∑h′∈ΞH′PRh′EHs−PRh′′EHs2subject to(3b)PRh′′EHs=∑i∈ΞEHγh′′EPi,h′′EH+γh′′GGi,h′′EH+γh′′HHi,h′′EH∀h′(3c)Constraints1b‐1w,andmodel2a‐2uconsideringh′insteadofh,andusing.′forparametersandvariables


### 2.4. Robust Optimization of the Proposed Scheme

Parameters of the energy price,γ and ρ; load, PD, QD, GD, and HD; renewable power, P¯RES; and parameters of EVs, including CREES, DRESS, EEES0, E¯EES, and S¯EES, in the problems described by (1)–(3) are uncertain. In general, most studies including [6–8] adopt stochastic programming to model uncertainties. Nonetheless, as obtaining a guaranteed solution requires the analysis of numerous scenarios [68–72], the computational time will be long. As a result, this subsection uses robust programming to model these parameters. It is worth noting that this programming finds the optimal solution for the proposed scheme in the worst-case scenario caused by the mentioned uncertainties [15, 73]. So, it needs a scenario so that the optimal operation point of EHs at different energy networks is obtained at the least possible time compared to the stochastic modeling of uncertainties [73]. As is observed in Sections 2.1 and 2.2, the time step for the proposed scheme in DA and RT operation cases is, respectively, 1 hour and 5 minutes; thus, achieving a low computation time for the suggested problem is of high importance. Another benefit of robust planning is that the status of the optimal point in the mentioned scheme is more suitable in other scenarios in comparison to that in the worst-case scenario. In other words, this type of modeling finds an optimal robust point for the proposed scheme [74, 75]. In the following, the BURO technique is incorporated to model the mentioned uncertainties [76]. This method chooses an integer value of the uncertainty parameter u between its upper and lower boundaries, i.e., between 1+χu¯ and 1−χu¯. In the worst-case scenario, the value of u based on its position in the optimization problem is equal to its upper or lower boundary limit so that the solution space of the problem becomes smaller than that of the deterministic model and a more suboptimal value of the objective function can be available in the robust model than the deterministic model [76]. For instance, as per (1a) and (1b), according to which the profit of EHs should be maximized, by choosing the lower boundary value for energy price γ in the worst-case scenario, it is expected that the profit in the robust model is less than that in the deterministic model. The inverse situation holds for energy price ρ. If this parameter is set at its upper boundary limit in the worst-case scenario, the energy cost of networks in the robust model will be higher than that in the deterministic model. As the objective function (2a) should be minimized, the terms 1−χ·γ, 1+χ·ρ, 1+χ·PD, 1+χ·QD, 1+χ·GD, 1+χ·HD, 1−χP¯RES, 1−χ·CREES, 1−χ·DRESS, 1−χ·EEES0, 1+χE¯EES, and 1−χS¯EES substitute γ, ρ, PD, QD, GD, HD, P¯RES, CREES, DRESS, EEES0, E¯EES, and S¯EES, respectively, in the BURO model. According to this situation, the solution space of the proposed problem in the worst-case scenario will be smaller than that in a scenario corresponding to the deterministic model.


## 2.1. DA Operation of EHs (First-Layer EMS and Upper-Level Problem Formulation)

This subsection addresses the optimal participation of EHs in the hourly DA energy market. In this problem, as given in (1a), maximization of the profit of EHs in the mentioned market is taken into account as the objective function. As per (1b), the profit of EHs is equal to the sum of their revenues in electricity, gas, and heating markets. In each of these markets, the revenue equals the energy price multiplied by power [12]. Note that if the value of power is positive, (1b) will represent the revenue of the EH in the mentioned market, but if it is negative, (1b) will denote the cost of EHs in this market. The balances between active, reactive, heating, and gas powers between the generators and consumers in the EH are modeled in constraints (1c)–(1f) [52]. Also, (1g)–(1w) represent the model of operation of sources, storages, and responsive loads in the EH. Constraints (1g)–(1j) give the formulation of CHP’s operation [52, 53], where constraints (1g) and (1h) calculate the heating and gas power of the CHP in terms of its active power. The limits on the CHP output in both electrical and thermal sections are considered in (1i) and (1j). Constraint (1k) calculates the input gas power of the boiler in terms of its output thermal power, while the limit on its output capacity is modeled (1l). The operation model of RESs is given in (1m) and (1n), where the generated active power of these sources is equal to the maximum capacity of the sources in accordance with weather conditions and this is due to low operation cost and environmental pollution as stated by (1m) [54]. The limit on simultaneous control of active and reactive power of these sources follows (1n), where it addresses the limit on apparent power of these sources [54]. The formulation of operation of EESs is expressed in (1o)–(1r) [55], which represent the charging and discharging limits, limit on the energy stored in the ESS, and limit on the charger capacity (apparent power) of this ESS, respectively. It is worth noting that the operation model of all EVs available in the parking lot is shown in constraints (1o)–(1r) [56, 57], except that parameters CRmax, DRmax, EEES0, and S¯EES change over time [58, 59]. Thus, the subscript h is used in the model of EV parking lot for these parameters. Additionally, (1q) is stated in the form of two constraints 0≤Ei,hEES0+∑τ=1hηE,ChPi,τCh−1/ηE,DisPi,τDis and ∑h∈ΞHEi,hEES0+ηE,ChPi,hCh−1/ηE,DisPi,hDis=E¯iEES in EV model, which refer to the limit on the energy stored in EVs and their energy demand in the operation horizon. For EV parking lot, parameters CRmax, DRmax, and S¯EES at hour h are found from CRmax=∑ev=1Nhcrev, DRmax=∑ev=1Nhdrev, and S¯EES=∑ev=1NhS¯evEV, where cr, dr, S¯EV, and Nh represent the charging rate, discharging rate, the capacity of EV charger, and the number of EVs connected to the EH during hour h [56]. Also, EEES0 at hour h is equal to the sum of total arrival (initial) energy of EVs that are recently connected to the parking lot in the EH. The operation model of TES is represented in (1s)–(1u) [52, 53], which has the same model as the EES does, except for thermal power variables of charging and discharging are used in these constraints and the limit on the charger capacity, (1r), cannot be used in the operation model of TES. Finally, the operation model of responsive loads in the electrical, gas, and thermal sections are given by constraints (1v) and (1w) [1]. In this type of DRP, consumers act according to (1a) corresponding to the price signal so that the profit of EHs in the energy market is maximized. The consumers participated in the DRP receive energy from the EH during low energy price, but they inject energy into the EH during high energy prices. As inexpensive (expensive) energy price is proportional to the off-peak (peak) hours, the proposed DRP can reduce consumers’ energy demand during peak intervals and supply this reduced amount of energy during off-peak hours [1]. Hence, constraint (1v) refers to the limit on the active, gas, and thermal power control of the DRP. Also, (1w) ensures that electrical (gas or thermal) energy decreased during peak hours is supplied during off-peak hours.Constraints (1i), (1n), and (1r) denote the linearized limit on the apparent power of CHP, RES, and ESS. It should be noted that these constraints are limited to a circular plane with the center on the origin and a radius of S, P2+Q2≤S2, and has a nonlinear format [60]. Yet, this limit can be approximated by a regular polygon with a model given by (1i), (1n), and (1r) [60]. Additionally, in the case the number of sides for this plane is high, the computational error will be negligible compared to that of the circular plane.
(1a)maxF1=∑h∈ΞHPRhEHssubject to(1b)PRhEHs=∑i∈ΞEHγhEPi,hEH+γhGGi,hEH+γhHHi,hEH∀h(1c)Pi,hEH=Pi,hCHP+Pi,hRES+Pi,hDR+Pi,hDis−Pi,hCh−Pi,hD∀i,h(1d)Qi,hEH=Qi,hCHP+Qi,hRES+Qi,hEES−Pi,hD∀i,h(1e)Hi,hEH=Hi,hCHP+Hi,hBO+Hi,hDR+Hi,hDis−Hi,hCh−Hi,hD∀i,h(1f)Gi,hEH=Gi,hDR−Gi,hCHP−Gi,hBO−Gi,hD∀i,h(1g)Hi,hCHP=1−ηt−ηlηhηtPi,hCHP∀i,h(1h)Gi,hCHP=1ηtPi,hCHP∀i,h(1i)Pi,hCHPcosm·Δβ+Qi,hCHPsinm·Δβ≤S¯iCHP∀i,h,m∈ΞM=1,2,⋯,nM,Δβ=2πnM(1j)0≤Hi,hCHP≤H¯iCHP∀i,h(1k)Gi,hBO=1ηbPi,hBO∀i,h(1l)0≤Hi,hBO≤H¯iBO∀i,h(1m)Pi,hRES=P¯iRES∀i,h(1n)Pi,hREScosm·Δβ+Qi,hRESsinm·Δβ≤S¯iRES∀i,h,m(1o)0≤Pi,hCh≤CRiEESxi,h∀i,h(1p)0≤Pi,hDis≤DRiEES1−xi,h∀i,h(1q)E¯iEES≤EiEES0+∑τ=1hηE,ChPi,τCh−1ηE,DisPi,τDis≤E¯iEES∀i,h(1r)Pi,hDis−Pi,hChcosm·Δβ+Qi,hEESsinm·Δβ≤S¯iEES∀i,h,m(1s)0≤Hi,hCh≤CRiTESyi,h∀i,h(1t)0≤Hi,hDis≤DRiTES1−yi,h∀i,h(1u)E¯iTES≤EiTES0+∑τ=1hηT,ChHi,τCh−1ηT,DisHi,τDis≤E¯iTES∀i,h(1v)−ξiPi,hDorHi,hDorGi,hD≤Pi,hDRorHi,hDRorGi,hDR≤ξiPi,hDorHi,hDorGi,hD∀i,h(1w)∑h∈ΞHPi,hDRorHi,hDRorGi,hDR=0∀i


## 2.2. DA Operation of Energy Networks (Second-Layer EMS and Lower-Level Formulation)

This subsection presents the hourly DA operation of electricity, gas, and heating networks in the presence of EHs. The problem is given in (2a)–(2u) and has an objective function aimed at minimizing the cost of energy purchased from the upstream electrical (first part), gas (second part), and heating (third section) networks [15]. Constraints of this problem as given by (2b)–(2u) refer to LOPF equations in the mentioned networks [12, 15]. Power flow equations are given in (2b)–(2l), where constraints (2b)–(2e) introduce power flow of the electricity network and refer to active and reactive power balance at the buses of the mentioned network and active and reactive power flowing through the line. Equations (2d) and (2e) are nonlinear in the forms of Pe,jEL=Ge,jVe2−VeVjGe,jcosφe−φj+Be,jsinφe−φj and Qe,jEL=−Be,jVe2+VeVjBe,jcosφe−φj−Ge,jsinφe−φj [60–64]. However, since the voltage phase angle between any two buses, φe−φj, in the distribution is less than 6° or 0.105 rad, the terms cosφe−φj and sinφe−φj are approximated by 1 and φe−φj, respectively. Moreover, since the magnitude of bus voltage should be between 0.9 and 1.1 p.u., it will be always close to 1 p.u. Hence, the voltage magnitude, V, is approximated as 1+ΔV, where ΔV represents the voltage deviation. According to these hypotheses and the small values of ΔV2 and ΔVφe−φj because of the small values of φe−φj and ΔV, constraints (2d) and (2e) replace the mentioned nonlinear constraints so that the model of LOPF is achieved [12]. Relations (2f)–(2j) concern power flow equations of the gas network, where the balance of gas power at different nodes is modeled in (2f). Also, (2g) calculates the gas power flow through the pipeline. In this equation, the squared power of the auxiliary variable σ, that is σ2, is in the form of signπg,h,πj,hπg,h2−πj,h2 [12, 15]. Hence, (2h)–(2j) present the linearized model of this equation based on the conventional piecewise linearization technique [60], where (2h) is the mentioned equation and (2i) and (2j) show the variables σ and π in terms of deviation variables of σ and π that are Δσ and Δπ, respectively. The formulation of power flow of the heating network is stated in (2k) and (2l), which, respectively, denote the balance of thermal power at different nodes and thermal power flowing through the pipeline.Constraints (2m)–(2u) present the operation limits of the different energy network [12, 15] and, respectively, refer to the limit on the voltage deviation of buses, limit on the capacity of the line (apparent power) and distribution substation, limit on the gas pressure at different nodes, limit on the capacity of pipeline (gas power) and gas station, limit on the temperature of nodes, and limit on the capacity (thermal power) of the pipeline and heating station. Note that, in the power flow model of the electricity network, given by (2b)–(2e), a voltage deviation variable is used to reach the LOPF; consequently, the limit on the voltage deviation given by (2m) substitutes the limit on voltage magnitude, V¯≤ΔVe≤V¯. Besides, the status of constraints (2n) and (2o) is according to the rule used in (1i).
(2a)minF2=∑h∈ΞHρhEPref,hES+ρhGGref,hES+ρhHHref,hESsubject to(2b)Pe,hES+∑i∈ΞEHIe,iE,EHPi,hEH−∑j∈ΞENIe,jELPe,j,hEL=Pe,hD∀e,h(2c)Qe,hES+∑i∈ΞEHIe,iE,EHQi,hEH−∑j∈ΞENIe,jELQe,j,hEL=Qe,hD∀e,h(2d)Pe,j,hEL=Ge,jΔVe,h−ΔVj,h−Be,jφe,h−φj,h∀e,j,h(2e)Qe,j,hEL=−Be,jΔVe,h−ΔVj,h−Ge,jφe,h−φj,h∀e,j,h(2f)Gg,hES+∑i∈ΞEHIg,iG,EHGi,hEH−∑j∈ΞGNIg,jGLGg,j,hGL=Gg,hD∀g,h(2g)Gg,j,hGL=κg,jsignπg,h,πj,hσg,j,h∀g,j,h(2h)Δσg,j,h,l=κg,jsignπg,h,πj,hSlπSlσΔπg,h,l−Δπj,h,l∀g,j,h,l(2i)σg,j,h=∑l∈ΞLΔσg,j,h,l∀g,j,h(2j)πg,h=π¯+∑l∈ΞLΔπg,h,l∀g,h(2k)Ht,hES+∑i∈ΞEHIt,iH,EHHi,hEH−∑j∈ΞHNIt,jHLHt,j,hHL=Ht,hD∀t,h(2l)Ht,j,hHL=ct,jṁt,jTt,h−Tj,h∀t,h(2m)1−V¯≤ΔVe,h≤V¯−1∀e,h(2n)Pe,j,hELcosm·Δβ+Qe,j,hELsinm·Δβ≤S¯e,jEL∀e,j,h,m(2o)Pe,hEScosm·Δβ+Qe,hESsinm·Δβ≤S¯eES∀e=ref,h,m(2p)π¯≤πg,h≤π¯∀g,h(2q)−G¯g,jGL≤Gg,j,hGL≤G¯g,jGL∀g,j,h(2r)−G¯gGS≤Gg,hGS≤G¯gGS∀g=ref,h(2s)T¯≤Tt,h≤T¯∀t,h(2t)−H¯t,jHL≤Ht,j,hHL≤H¯t,jHL∀t,j,h(2u)−H¯tHS≤Ht,hHS≤H¯tHS∀t=ref,h


## 2.3. Five-Minute RT Scheduling of EHs in Energy Networks

This subsection presents the optimal participation of EHs in the 5 min RT energy market subject to the condition that the unbalance between RT and DA operations is the minimum. It is worth noting that EHs consist of RESs. Since the active power generation by these sources depends on weather conditions and it is a challenging task with generally forecast errors, the operation results are expected to differ due to the uncertainty in the active power generation. This can lead to unbalance between supply and demand in the RT operation, and this is considered for low flexibility of the system [65, 66]. To compensate for this issue, the exploitation of flexibility sources, storage and DRP, has been recommended by many scholars [65]. Yet, this situation can be extracted to establish a proper optimization framework; so, this subsection, as given by (3a), the minimization of unbalance between the EHs’ profit in DA and RT energy markets is considered as the objective function [66]. The EHs’ profit in the RT market can be found from (3b). Also, the constraints of this problem as given in (3c) include the operation model of EHs, constraints of problem (1a)–(1w), and operation model of energy networks (problem 2). Note that the objective function in this problem is quadratic and its constraints have a linearized formulated framework. According to [52, 67], this model can be solved using linear programming (LP) solvers such as CPLEX in the least possible computation time.
(3a)minF3=∑h′∈ΞH′PRh′EHs−PRh′′EHs2subject to(3b)PRh′′EHs=∑i∈ΞEHγh′′EPi,h′′EH+γh′′GGi,h′′EH+γh′′HHi,h′′EH∀h′(3c)Constraints1b‐1w,andmodel2a‐2uconsideringh′insteadofh,andusing.′forparametersandvariables


## 2.4. Robust Optimization of the Proposed Scheme

Parameters of the energy price,γ and ρ; load, PD, QD, GD, and HD; renewable power, P¯RES; and parameters of EVs, including CREES, DRESS, EEES0, E¯EES, and S¯EES, in the problems described by (1)–(3) are uncertain. In general, most studies including [6–8] adopt stochastic programming to model uncertainties. Nonetheless, as obtaining a guaranteed solution requires the analysis of numerous scenarios [68–72], the computational time will be long. As a result, this subsection uses robust programming to model these parameters. It is worth noting that this programming finds the optimal solution for the proposed scheme in the worst-case scenario caused by the mentioned uncertainties [15, 73]. So, it needs a scenario so that the optimal operation point of EHs at different energy networks is obtained at the least possible time compared to the stochastic modeling of uncertainties [73]. As is observed in Sections 2.1 and 2.2, the time step for the proposed scheme in DA and RT operation cases is, respectively, 1 hour and 5 minutes; thus, achieving a low computation time for the suggested problem is of high importance. Another benefit of robust planning is that the status of the optimal point in the mentioned scheme is more suitable in other scenarios in comparison to that in the worst-case scenario. In other words, this type of modeling finds an optimal robust point for the proposed scheme [74, 75]. In the following, the BURO technique is incorporated to model the mentioned uncertainties [76]. This method chooses an integer value of the uncertainty parameter u between its upper and lower boundaries, i.e., between 1+χu¯ and 1−χu¯. In the worst-case scenario, the value of u based on its position in the optimization problem is equal to its upper or lower boundary limit so that the solution space of the problem becomes smaller than that of the deterministic model and a more suboptimal value of the objective function can be available in the robust model than the deterministic model [76]. For instance, as per (1a) and (1b), according to which the profit of EHs should be maximized, by choosing the lower boundary value for energy price γ in the worst-case scenario, it is expected that the profit in the robust model is less than that in the deterministic model. The inverse situation holds for energy price ρ. If this parameter is set at its upper boundary limit in the worst-case scenario, the energy cost of networks in the robust model will be higher than that in the deterministic model. As the objective function (2a) should be minimized, the terms 1−χ·γ, 1+χ·ρ, 1+χ·PD, 1+χ·QD, 1+χ·GD, 1+χ·HD, 1−χP¯RES, 1−χ·CREES, 1−χ·DRESS, 1−χ·EEES0, 1+χE¯EES, and 1−χS¯EES substitute γ, ρ, PD, QD, GD, HD, P¯RES, CREES, DRESS, EEES0, E¯EES, and S¯EES, respectively, in the BURO model. According to this situation, the solution space of the proposed problem in the worst-case scenario will be smaller than that in a scenario corresponding to the deterministic model.


## 3. Solution Procedure

The problem presented in the previous section is a two-stage optimization problem, where the hourly DA operation of EHs in the electricity, gas, and heating networks in first stage is investigated based on Figure1. The results obtained from the first stage problem, i.e., PREHs, are applied to the second stage problem. This stage finds the 5 min RT operation of EHs in the mentioned networks. Moreover, the formulation of the first and second stages is in the form of bilevel optimization, and according to Figure 1, the EH power should be met based on the technical limits of the mentioned networks, given by (2m)–(2u). Hence, each of these problems (i.e., the first and second stages) can be expressed in the form of the model presented in
(4a)maxa1Tz+mina2Tws.to:b1w+b2z≤b3:μc1w≤c2:⏞LowerlevelSecondlayersubject to(4b)d1z≤d2.The first part of (4a) and constraint (4b) refer to the formulation of the upper-level problem (Section 2.1) and represent (1a) (or (3a)) and (1b)–(1w), respectively. In these equations, variable z is the vector of variables of the upper-level problem, i.e., PREHsPEHQEHGEHHEHPCHPQCHPGCHPHCHPPRESQRESPDisPChQEESPDRGDRHDRGBOHBOHDisHCh. The second part of (4a) refers to the formulation of the lower-level problem (Section 2.2), where hourly DA (RT) operation of energy networks is considered in the first (second) stage problem. In this equation, variable w denotes that vector of variables of the lower-level problem, i.e., PESQESPELQELGGSGGLHHSHHLΔVφσΔσπΔπT. Furthermore, λ and μ show the dual variables of the lower-level problem constraints, formulation (2a)–(2u). To solve the problem described by (4), first of all, the dual model of the lower-level problem, i.e., the second part of (4a), should be obtained. It should be mentioned that the value of z is fixed in this problem, and the values of dual variables in the upper-level problem are also fixed. Accordingly, the formulation of the dual of the second part in (4a) will be as follows [73, 77]:
(5)Ψz=maxΓμ,λb3−b2zTμ+c2Tλ∀Γμ,λ=μ,λb1μ+c1λ≤=≥a2,μ≥0,λ≥0,where Ψz and Γμ,λ represent the optimal solution and solution space for the dual formulation of the lower-level problem, respectively. Also, choosing the terms ≤/=/≥ in the set Γμ,λ depends on the type of variable w; if it is a positive/free/negative variable, the signs ≥/=/≤ will be used, respectively [77]. Considering (5) in the model described by (4), the formulation (4a)–(4b) can be stated as follows [72]:
(6)maxΓza1Tz+Ψz∀Γz=zd1z≤d2.Equation (6) can calculate variable z, while the function can be extracted from (5). Hence, the BD method, which has generally the least possible computation time [78], is employed to solve this problem. The steps for solving the problem described by (5) and (6) using the BD method are summarized as follows [73, 77, 78].Regarding the solution process adopted for solving the first stage problem, the initial values of variablesPEH, QEH, GEH, and HEH are set zero. The optimal values of dual variables of constraints (2b)–(2u) are extracted by solving the dual formulation of the lower-level problem (Section 2.2), at the first iteration, n=1. In the next step, the Benders cut, i.e., ω≤b3−b2zμkT+c2λkT, is added to the DA operation model of EHs, problem (1a)–(1w). Note that Benders cut for the proposed scheme is as given in (7). In this case, the new optimal values of PEH, QEH, GEH, and HEH are calculated. Problem (5) is then solved based on these new values. This process continues until the convergence points are found for the BD algorithm.
(7)∑e,hdualvaribleofEq.2b×Pe,hD−∑i∈ΞEHIe,iE,EHPi,hEH+∑e,hdualvaribleofEq.2c×Qe,hD−∑i∈ΞEHIe,iE,EHQi,hEH+∑g,hdualvaribleofEq.2f×Gg,hD−∑i∈ΞEHIg,iG,EHGi,hEH+∑t,hdualvaribleofEq.2k×Ht,hD−∑i∈ΞEHIt,iH,EHHi,hEH+∑e,g,thdualvarible×constantparameterforEqs.2j,2m–2u.The results of the first stage are applied to the formulation of the second stage. Then, the model of the second stage obtains the optimal solution using a solution procedure similar to that used for the problem of the first stage. Finally, the BD process for the proposed scheme is summarized in Algorithm1.Algorithm 1: BD process for the proposed scheme.
Initialization: choose a value (z0) for z and solve Ψz0 using Equation (5) to find the values of μ1 and λ1. Set the lower/upper-level boundary of BD algorithm, zlower/zupper, to −∞/+∞. The number of iterations is set n=1.Forn≥1Step 1: solve the master problem of the BD approach modeled as mina1Tz+ω and s.tod1z≤d2 and ω≤b3−b2zμkT+c2λkT∀k≤n. The problem described finds the optimal values of zn and ωn. Set zlower=c1Tzn+ωn.Step 2: solve the subproblem of the BD approach to find the optimal situation of Ψzn, μn+1, and λn+1 modeled by (5). Then, set zupper=c1Tzn+Ψzn.Step 3: check the convergence condition of the BD technique. Problem (4a)–(4b) includes the optimal solution if zupper–zlower≤ε, where ε>0 denotes the BD tolerance. Otherwise, set n=n+1, and go to Step 1.


## 4. Numerical Results and Discussion

### 4.1. Case Study

The problem described in previous sections is applied to a sample system that consists of a 9-bus electrical network, a 4-node gas network, and a 7-node heating network, as shown in Figure2 [12]. The base power is 1 MVA, 1 MW, and 1 MW for these networks, respectively. Base voltage, pressure, and temperature are 1 kV, 10 bar, and 100°C. The upper and lower limits of these three parameters are 1.1 p.u. and 0.9 p.u. Characteristics of distribution lines and pipers together with electricity, gas, and heating network are reported in [12]. The data of peak electrical and thermal loads are extracted from [12], and the gas network excludes passive loads. The amount of load during other hours is equal to the multiplication of load factor and peak load. Figure 3(a) illustrates the daily load factor curve for electricity and gas networks [12, 15].Figure 2
Proposed test system in the presence of EH (hub) [12].Figure 3
Daily curve of (a) rate of load power [12], RES active generation power [77–81], and EV penetration [12] and (b) price of energy market [12].
(a)(b)The system in Figure2 includes 7 EHs. EHs 1–3 and 5 possess only electrical elements such as RES, battery, EVs, and electrical responsive loads. EH 4 has a CHP, a boiler, a TES, and responsive thermal loads. EHs 6 and 7 include all the elements mentioned for problem (1a)–(1w). The places of these EHs in the mentioned energy networks are specified in Figure 2. The data of sources, batteries, and EVs are reported in [12]. The data of TES is similar to that of the battery except its charging and discharging efficiencies are 80%. Also, the participation factor of EHs’ consumers in the DRP is 30%. The hourly active power generation by RESs is equal to the multiplication of power generation rate and size of RESs. The size of RES and number of EVs are given in [12]. Figure 3(a) depicts the curve of daily power generation rate for RESs of PV and WT types and daily EV penetration rate curve. Energy price curves in electrical, gas, and thermal markets are shown in Figure 3(b) [77].


### 4.2. Results

Simulation of the proposed problem is conducted in GAMS optimization software [66] in accordance with the data provided in previous sections. CPLEX solver [67] is used to find the optimal solution. Following [12], the number of linear pieces in the conventional piecewise linear method is set 5 in this paper, and a regular 45-gon is used instead of a circular plane. According to the results in [12, 15], the calculation error for electricity power, voltage, gas network variables, and heating system variables is 2%, 0.5%, 0.9%, and zero, respectively, in the linearized model of network operation in the presence of EHs with respect to its original nonlinear formulation (as nonlinear model of problem 2).


#### 4.2.1. Evaluation of the Convergence of the BD Algorithm for the Proposed Scheme

The convergence plot of the problem of EH participation in DA and RT energy markets based on the BD algorithm is shown in Figure4 for different values of uncertainty level (UL). Referring to Figure 4(a), the difference between the upper and lower values of the BD algorithm zupper−zlower in the robust model of problems (1a)–(1w) and (2a)–(2u) for different UL values from 0 to 0.1 has a value less than the permissible tolerance (ε) of the BD at the sixth convergence iteration. The value of ε in this paper is set $0.5. Thus, this problem reaches the optimal point at the sixth iteration of the BD solver under the mentioned conditions. But, based on Figure 4(a), its computational time for different ULs differs. For UL=0, the computational time is 2.41 s, and for UL=0.1, it is 2.53 s. Based on Section 2.4, the increase in the UL means a reduction in the solution space of the problem in the worst-case scenario in comparison to a scenario corresponding to the deterministic model (with UL=0). In this case, the robust optimal solution is achieved at a longer computational time. Yet, it is worth noting that, based on Figure 4(a), increasing the UL to 0.2 leads to an increase in the convergence iteration number of the BD and computational time for the robust problems (1a)–(1w) and (2a)–(2u) compared to those of the deterministic model. These conditions are due to the reduced solution space of the mentioned problem compared to UL=0. Therefore, a small increase in the UL with respect to UL=0 only increases the computational time of the problem, while a large increase in the UL increases both computational time and the convergence iteration number of the BD. The mentioned conditions for the robust problem described in (1) and (2) are also true for achieving the optimal solution of the robust problem in which EHs participate in the RT energy market (3), as is shown in Figure 4(b). Nonetheless, this problem suffers from a high number of convergence iteration and longer computational time compared to problems (1a)–(1w) and (2a)–(2u) because its time step is 5 minutes; thus, 288 time steps are needed for a 24-hour operation horizon. The problem of EH participation in the DA energy market has a time step of 1 hour, so there are 24 time steps. Following this, the volume of problem (3a)–(3c) is greater than the formulation provided in (1) and (2) and its computational time and convergence iteration number are higher than those of formulations (1a)–(1w) and (2a)–(2u).Figure 4
BD convergence curve in the problem of (a) the first stage and (b) the second stage for different values of uncertainty level.
(a)(b)


#### 4.2.2. Evaluation of the Capability of the Two-Layer EMS Scheme for Energy Networks with EHs

The capability of three different models of EMS on electricity, gas, and heating networks with EHs is evaluated in terms of the deterministic model (UL=0). The details are provided below.
(i)
Uncoordinated model (UM) of sources and active loads in the two-layer EMS framework: it is assumed in this section that there is no coordinated form like an EH to integrate sources and active loads (ESSs and responsive loads) and they are managed separately by their operator. Therefore, the proposed scheme (1)–(3) is solved for individual elements of the EH separately. In other words, to examine the capability of the CHP, only the model of this element is used in the formulation (1a)–(1w) and models of other elements are not employed [12].(ii)
Coordinated model (CM) of sources and active loads in a two-layer EMS framework: in this design, various sources and active loads are coordinated in the form of an EH and are managed by the EH operator. Also, operators of EHs are mutually coordinated with those of energy networks to meet the second layer of the EMS. The formulation is similar to problems (1)–(3).(iii)
Coordinated model (CM) of sources and active loads in a single-layer EMS framework: it is assumed in this design that elements of EHs are mutually coordinated with operators of energy networks. The modeling of this design in the DA operation includes the objective function (1) and constraints of the problem described by (1) and (2). It means that the objective function (2a) is removed from problems (1a)–(1w) and (2a)–(2u) and the remaining relations form a single-layer problem. Such conditions hold for problem (3a)–(3c). Since the proposed problem is single-layer and linear, CPLEX is used to solve it and the BD algorithm cannot be executed on that.The results obtained in these parts of the study are tabulated in Tables2 and 3. Based on Table 2, the total profit of sources and active loads in the UM method for DA and RT energy markets is $406.4 and $409.3. But these values for the CM method in the mentioned markets are $456.12 and $457.94, meaning that the CM design for sources and active loads in the form of an EH succeeds to enhance the profit of mentioned elements in DA and RT markets by up to 12.23% and 11.88% compared to a case they were managed separately in the UM design. The value of unbalance function (3a) between these markets is $22.9 (406.4-409.3) for the UM design, while it decreases to $21.82 in the CM design. A comparison between the CM design with single- and two-layer EMSs shows that the profit of sources and active loads in the EH form have negligible differences. The difference between the profit of EHs in the mentioned two designs for DA and RT energy markets is $0.46 (456.12-455.66) and 0.41% (457.53-457.94). Referring to Table 2, this is a small value with respect to the profit obtained in the mentioned markets. Yet, note that the computational time of the proposed problem for the single-layer EMS is roughly 73.5% (8.45−4.87/4.87) higher than that of the proposed design model based on the two-layer EMS. For the two-layer EMS case, all data are divided into two problems, such as problem (1a)–(1w) and problem (2a)–(2u); but, in the case of single-layer EMS, all data is integrated into a single problem and so the computational time will be shorter for the two-layer EMS than the single-layer model. Another advantage of the two-layer EMS compared to the single-layer EMS is that it models objectives of energy networks in the objective function of the lower-level problem, while it was removed in the single-layer EMS. Consequently, according to Table 2 and the discussions made above, the CM design for sources and active loads corresponding to the two-layer EMS is superior economically and in terms of computational time to other EMS models. Table 3 reports the financial benefits of EHs for different ULs in the case of the two-layer EMS. As per this table, the increased UL reduces the EHs’ profit in both DA and RT energy markets because, based on Section 2.4, increasing the UL will increase energy demand of passive loads and EVs as well as reduce the energy generation capacity of RESs in the worst-case scenario compared to the deterministic model (UL=0). Additionally, increasing the UL increases the value of unbalance function (3a) between the two mentioned markets. As another point, note that considering the robust model of uncertainties in the worst-case scenarios (UL=0.3) imposes 107.69% (456.12-348.43) and $107.36 (457.94-350.58) in DA and RT energy markets compared to the deterministic model. This cost is caused by considering the prediction error of uncertainty parameters.Table 2
Profit of EHs in the different energy management methods forUL=0.
Energy management methodProfit ($) of EHs in the energy marketCalculation time (s)DARTTwo-layer uncoordinated model (UM)CHP149.84149.113.41RES47.8748.553.08Boiler89.4690.233.33EVs-36.71-36.023.44EES58.6459.883.43TES51.5750.883.38DRP45.7346.673.12Sum406.40409.30—Two-layer coordinated model (CM)Total456.12457.944.87Gain compared to UM12.23%11.88%One-layer coordinated modelTotal455.66457.538.45Table 3
Unbalanced situation between DA and RT energy markets for different values of uncertainty level.
UL00.10.20.3Profit ($) in energy marketDA456.12421.76385.55348.43RT457.94423.73387.61350.58Unbalanced function, (3a), in $23.3123.8814.2444.622


#### 4.2.3. Analysis of the Operation of EHs

The daily power curve for EHs participating in the markets is shown in Figure5. As per Figure 5(a), the negative value of active power generation by EHs between 1:00 and 8:00 means that EVs, batteries, and loads are charged. In parallel to this period, the electricity price is low and thus EHs benefit higher profit. Sources, ESS, and DRPs transfer power into the electricity network during 9:00-24:00 as the value of power is positive in this time interval; hence, higher benefits are gained. The power injection by EHs into the network is high from 17:00 to 00:00 in comparison to the 9:00-16:00 interval because EVs, ESSs, and loads are discharged. Note that CHPs and RESs transfer power to the electrical network during 17:00-24:00. The reason is that the price of electricity is high [51]. At 1:00-8:00 and 17:00-22:00 intervals, EHs transfer reactive power to the electricity network based on Figure 5(b). In terms of heating power, EHs inject it into the heating network during 8:00-15:00 and 23:00-24:00 or 1:00/7:00 as well as 16:00-22:00 (Figure 5(c)). Since the price of heating power is low at these intervals, the TES is charged and together with the boiler and CHP inject the heating power to the heating network at periods 8:00-15:00 and 23:00-24:00, aimed at gaining more profit for EHs. The TES is disconnected from 16:00 to 22:00, and the heating power generation by CHP and boiler is low. Gas power of EHs is based on Figure 5(d). The profit curve of EHs for individual days can be seen in Figure 6. The amount of profit relies on the prices of active power, gas power, and heating power. During 1:00-8:00, the profit obtained by EHs in the markets has a negative value because active power, gas power, and heating power are low in this interval (Figure 5). Nevertheless, profit is positive at the rest of the hours.Figure 5
Daily curve of (a) active power, (b) reactive power, (c) heating power, and (d) gas power for all EHs in the DA and RT energy markets forUL=0.2.
(a)(b)(c)(d)Figure 6
Daily EH profit curve in DA and RT energy markets forUL=0.Besides, as the figures show, the unbalance between DA and RT operations during 12:00-15:00 is higher than other hours for EH active power and their profit. The unbalance is due to the prediction error of power generation by RESs that lack power control capability. Since, based on Figure3(a), RESs inject high amounts of power into EHs, the unbalance between these operations will be higher than that in other hours. According to Figure 5(b), there is no unbalance reactive power in both DA and RT operations, and there is the same distance between DA and RT scheduling at all operation times for EH gas and heat powers based on Figures 5(c) and 5(d).


#### 4.2.4. Analysis of Operation of Energy Networks

Table4 presents operation index values for electricity, gas, and heating networks for different values of ULs concerning the two-layer EMS for two case studies: power flow study (Case I) and the proposed scheme (Case II). The table also lists the operation cost of networks (2a), the energy loss of networks (the difference between energy generation and consumption), maximum voltage drop (MVD), maximum pressure drop (MPD), maximum temperature drop (MTD), maximum overvoltage (MOV), maximum overpressure (MOP), and maximum overtemperature (MOT) in the case of RT operation. Based on the reported results in this paper, it is observed that the proposed scheme in the worst-case scenario of uncertainties with UL=0.3 succeeds to reduce the operating cost of networks by almost 18.4% (4151.7−3387.4/4151.7) compared to that in Case I. Also, it helps to degrade energy loss, MVD, and MTD values by roughly 20.5%, 26.8%, and 18.2% in comparison to Case I. However, these conditions are realized for increasing MOV, MOT, and MPD by about 0.005, 0.007, and 0.089 p.u. These values of parameters are smaller than their allowed limit of 0.1 p.u. (1.1-1). Moreover, the status of mentioned indices for the case the UL decreases becomes more suitable than the case with UL=0.3 because, under such conditions and based on Section 2.4, the capacity of energy generation by RESs increases and energy demand of passive loads and EVs decreases.Table 4
Values of operation indices in the RT scheduling of networks for different cases and various values of UL.
CaseIIIUL=0UL=0.1UL=0.2UL=0.3Operation cost, (2a), in $4151.72861.23013.53191.23387.4Energy loss (p.u.)7.245.175.355.545.76MVD (p.u.)0.1120.0780.0790.0810.082MPD (p.u.)00.0840.0860.0870.089MTD (p.u.)0.1100.0870.0880.0890.090MOV (p.u.)00.0090.0070.0060.005MOP (p.u.)00000MOT (p.u.)00.01020.0100.0080.007


## 4.1. Case Study

The problem described in previous sections is applied to a sample system that consists of a 9-bus electrical network, a 4-node gas network, and a 7-node heating network, as shown in Figure2 [12]. The base power is 1 MVA, 1 MW, and 1 MW for these networks, respectively. Base voltage, pressure, and temperature are 1 kV, 10 bar, and 100°C. The upper and lower limits of these three parameters are 1.1 p.u. and 0.9 p.u. Characteristics of distribution lines and pipers together with electricity, gas, and heating network are reported in [12]. The data of peak electrical and thermal loads are extracted from [12], and the gas network excludes passive loads. The amount of load during other hours is equal to the multiplication of load factor and peak load. Figure 3(a) illustrates the daily load factor curve for electricity and gas networks [12, 15].Figure 2
Proposed test system in the presence of EH (hub) [12].Figure 3
Daily curve of (a) rate of load power [12], RES active generation power [77–81], and EV penetration [12] and (b) price of energy market [12].
(a)(b)The system in Figure2 includes 7 EHs. EHs 1–3 and 5 possess only electrical elements such as RES, battery, EVs, and electrical responsive loads. EH 4 has a CHP, a boiler, a TES, and responsive thermal loads. EHs 6 and 7 include all the elements mentioned for problem (1a)–(1w). The places of these EHs in the mentioned energy networks are specified in Figure 2. The data of sources, batteries, and EVs are reported in [12]. The data of TES is similar to that of the battery except its charging and discharging efficiencies are 80%. Also, the participation factor of EHs’ consumers in the DRP is 30%. The hourly active power generation by RESs is equal to the multiplication of power generation rate and size of RESs. The size of RES and number of EVs are given in [12]. Figure 3(a) depicts the curve of daily power generation rate for RESs of PV and WT types and daily EV penetration rate curve. Energy price curves in electrical, gas, and thermal markets are shown in Figure 3(b) [77].


## 4.2. Results

Simulation of the proposed problem is conducted in GAMS optimization software [66] in accordance with the data provided in previous sections. CPLEX solver [67] is used to find the optimal solution. Following [12], the number of linear pieces in the conventional piecewise linear method is set 5 in this paper, and a regular 45-gon is used instead of a circular plane. According to the results in [12, 15], the calculation error for electricity power, voltage, gas network variables, and heating system variables is 2%, 0.5%, 0.9%, and zero, respectively, in the linearized model of network operation in the presence of EHs with respect to its original nonlinear formulation (as nonlinear model of problem 2).


### 4.2.1. Evaluation of the Convergence of the BD Algorithm for the Proposed Scheme

The convergence plot of the problem of EH participation in DA and RT energy markets based on the BD algorithm is shown in Figure4 for different values of uncertainty level (UL). Referring to Figure 4(a), the difference between the upper and lower values of the BD algorithm zupper−zlower in the robust model of problems (1a)–(1w) and (2a)–(2u) for different UL values from 0 to 0.1 has a value less than the permissible tolerance (ε) of the BD at the sixth convergence iteration. The value of ε in this paper is set $0.5. Thus, this problem reaches the optimal point at the sixth iteration of the BD solver under the mentioned conditions. But, based on Figure 4(a), its computational time for different ULs differs. For UL=0, the computational time is 2.41 s, and for UL=0.1, it is 2.53 s. Based on Section 2.4, the increase in the UL means a reduction in the solution space of the problem in the worst-case scenario in comparison to a scenario corresponding to the deterministic model (with UL=0). In this case, the robust optimal solution is achieved at a longer computational time. Yet, it is worth noting that, based on Figure 4(a), increasing the UL to 0.2 leads to an increase in the convergence iteration number of the BD and computational time for the robust problems (1a)–(1w) and (2a)–(2u) compared to those of the deterministic model. These conditions are due to the reduced solution space of the mentioned problem compared to UL=0. Therefore, a small increase in the UL with respect to UL=0 only increases the computational time of the problem, while a large increase in the UL increases both computational time and the convergence iteration number of the BD. The mentioned conditions for the robust problem described in (1) and (2) are also true for achieving the optimal solution of the robust problem in which EHs participate in the RT energy market (3), as is shown in Figure 4(b). Nonetheless, this problem suffers from a high number of convergence iteration and longer computational time compared to problems (1a)–(1w) and (2a)–(2u) because its time step is 5 minutes; thus, 288 time steps are needed for a 24-hour operation horizon. The problem of EH participation in the DA energy market has a time step of 1 hour, so there are 24 time steps. Following this, the volume of problem (3a)–(3c) is greater than the formulation provided in (1) and (2) and its computational time and convergence iteration number are higher than those of formulations (1a)–(1w) and (2a)–(2u).Figure 4
BD convergence curve in the problem of (a) the first stage and (b) the second stage for different values of uncertainty level.
(a)(b)


### 4.2.2. Evaluation of the Capability of the Two-Layer EMS Scheme for Energy Networks with EHs

The capability of three different models of EMS on electricity, gas, and heating networks with EHs is evaluated in terms of the deterministic model (UL=0). The details are provided below.
(i)
Uncoordinated model (UM) of sources and active loads in the two-layer EMS framework: it is assumed in this section that there is no coordinated form like an EH to integrate sources and active loads (ESSs and responsive loads) and they are managed separately by their operator. Therefore, the proposed scheme (1)–(3) is solved for individual elements of the EH separately. In other words, to examine the capability of the CHP, only the model of this element is used in the formulation (1a)–(1w) and models of other elements are not employed [12].(ii)
Coordinated model (CM) of sources and active loads in a two-layer EMS framework: in this design, various sources and active loads are coordinated in the form of an EH and are managed by the EH operator. Also, operators of EHs are mutually coordinated with those of energy networks to meet the second layer of the EMS. The formulation is similar to problems (1)–(3).(iii)
Coordinated model (CM) of sources and active loads in a single-layer EMS framework: it is assumed in this design that elements of EHs are mutually coordinated with operators of energy networks. The modeling of this design in the DA operation includes the objective function (1) and constraints of the problem described by (1) and (2). It means that the objective function (2a) is removed from problems (1a)–(1w) and (2a)–(2u) and the remaining relations form a single-layer problem. Such conditions hold for problem (3a)–(3c). Since the proposed problem is single-layer and linear, CPLEX is used to solve it and the BD algorithm cannot be executed on that.The results obtained in these parts of the study are tabulated in Tables2 and 3. Based on Table 2, the total profit of sources and active loads in the UM method for DA and RT energy markets is $406.4 and $409.3. But these values for the CM method in the mentioned markets are $456.12 and $457.94, meaning that the CM design for sources and active loads in the form of an EH succeeds to enhance the profit of mentioned elements in DA and RT markets by up to 12.23% and 11.88% compared to a case they were managed separately in the UM design. The value of unbalance function (3a) between these markets is $22.9 (406.4-409.3) for the UM design, while it decreases to $21.82 in the CM design. A comparison between the CM design with single- and two-layer EMSs shows that the profit of sources and active loads in the EH form have negligible differences. The difference between the profit of EHs in the mentioned two designs for DA and RT energy markets is $0.46 (456.12-455.66) and 0.41% (457.53-457.94). Referring to Table 2, this is a small value with respect to the profit obtained in the mentioned markets. Yet, note that the computational time of the proposed problem for the single-layer EMS is roughly 73.5% (8.45−4.87/4.87) higher than that of the proposed design model based on the two-layer EMS. For the two-layer EMS case, all data are divided into two problems, such as problem (1a)–(1w) and problem (2a)–(2u); but, in the case of single-layer EMS, all data is integrated into a single problem and so the computational time will be shorter for the two-layer EMS than the single-layer model. Another advantage of the two-layer EMS compared to the single-layer EMS is that it models objectives of energy networks in the objective function of the lower-level problem, while it was removed in the single-layer EMS. Consequently, according to Table 2 and the discussions made above, the CM design for sources and active loads corresponding to the two-layer EMS is superior economically and in terms of computational time to other EMS models. Table 3 reports the financial benefits of EHs for different ULs in the case of the two-layer EMS. As per this table, the increased UL reduces the EHs’ profit in both DA and RT energy markets because, based on Section 2.4, increasing the UL will increase energy demand of passive loads and EVs as well as reduce the energy generation capacity of RESs in the worst-case scenario compared to the deterministic model (UL=0). Additionally, increasing the UL increases the value of unbalance function (3a) between the two mentioned markets. As another point, note that considering the robust model of uncertainties in the worst-case scenarios (UL=0.3) imposes 107.69% (456.12-348.43) and $107.36 (457.94-350.58) in DA and RT energy markets compared to the deterministic model. This cost is caused by considering the prediction error of uncertainty parameters.Table 2
Profit of EHs in the different energy management methods forUL=0.
Energy management methodProfit ($) of EHs in the energy marketCalculation time (s)DARTTwo-layer uncoordinated model (UM)CHP149.84149.113.41RES47.8748.553.08Boiler89.4690.233.33EVs-36.71-36.023.44EES58.6459.883.43TES51.5750.883.38DRP45.7346.673.12Sum406.40409.30—Two-layer coordinated model (CM)Total456.12457.944.87Gain compared to UM12.23%11.88%One-layer coordinated modelTotal455.66457.538.45Table 3
Unbalanced situation between DA and RT energy markets for different values of uncertainty level.
UL00.10.20.3Profit ($) in energy marketDA456.12421.76385.55348.43RT457.94423.73387.61350.58Unbalanced function, (3a), in $23.3123.8814.2444.622


### 4.2.3. Analysis of the Operation of EHs

The daily power curve for EHs participating in the markets is shown in Figure5. As per Figure 5(a), the negative value of active power generation by EHs between 1:00 and 8:00 means that EVs, batteries, and loads are charged. In parallel to this period, the electricity price is low and thus EHs benefit higher profit. Sources, ESS, and DRPs transfer power into the electricity network during 9:00-24:00 as the value of power is positive in this time interval; hence, higher benefits are gained. The power injection by EHs into the network is high from 17:00 to 00:00 in comparison to the 9:00-16:00 interval because EVs, ESSs, and loads are discharged. Note that CHPs and RESs transfer power to the electrical network during 17:00-24:00. The reason is that the price of electricity is high [51]. At 1:00-8:00 and 17:00-22:00 intervals, EHs transfer reactive power to the electricity network based on Figure 5(b). In terms of heating power, EHs inject it into the heating network during 8:00-15:00 and 23:00-24:00 or 1:00/7:00 as well as 16:00-22:00 (Figure 5(c)). Since the price of heating power is low at these intervals, the TES is charged and together with the boiler and CHP inject the heating power to the heating network at periods 8:00-15:00 and 23:00-24:00, aimed at gaining more profit for EHs. The TES is disconnected from 16:00 to 22:00, and the heating power generation by CHP and boiler is low. Gas power of EHs is based on Figure 5(d). The profit curve of EHs for individual days can be seen in Figure 6. The amount of profit relies on the prices of active power, gas power, and heating power. During 1:00-8:00, the profit obtained by EHs in the markets has a negative value because active power, gas power, and heating power are low in this interval (Figure 5). Nevertheless, profit is positive at the rest of the hours.Figure 5
Daily curve of (a) active power, (b) reactive power, (c) heating power, and (d) gas power for all EHs in the DA and RT energy markets forUL=0.2.
(a)(b)(c)(d)Figure 6
Daily EH profit curve in DA and RT energy markets forUL=0.Besides, as the figures show, the unbalance between DA and RT operations during 12:00-15:00 is higher than other hours for EH active power and their profit. The unbalance is due to the prediction error of power generation by RESs that lack power control capability. Since, based on Figure3(a), RESs inject high amounts of power into EHs, the unbalance between these operations will be higher than that in other hours. According to Figure 5(b), there is no unbalance reactive power in both DA and RT operations, and there is the same distance between DA and RT scheduling at all operation times for EH gas and heat powers based on Figures 5(c) and 5(d).


### 4.2.4. Analysis of Operation of Energy Networks

Table4 presents operation index values for electricity, gas, and heating networks for different values of ULs concerning the two-layer EMS for two case studies: power flow study (Case I) and the proposed scheme (Case II). The table also lists the operation cost of networks (2a), the energy loss of networks (the difference between energy generation and consumption), maximum voltage drop (MVD), maximum pressure drop (MPD), maximum temperature drop (MTD), maximum overvoltage (MOV), maximum overpressure (MOP), and maximum overtemperature (MOT) in the case of RT operation. Based on the reported results in this paper, it is observed that the proposed scheme in the worst-case scenario of uncertainties with UL=0.3 succeeds to reduce the operating cost of networks by almost 18.4% (4151.7−3387.4/4151.7) compared to that in Case I. Also, it helps to degrade energy loss, MVD, and MTD values by roughly 20.5%, 26.8%, and 18.2% in comparison to Case I. However, these conditions are realized for increasing MOV, MOT, and MPD by about 0.005, 0.007, and 0.089 p.u. These values of parameters are smaller than their allowed limit of 0.1 p.u. (1.1-1). Moreover, the status of mentioned indices for the case the UL decreases becomes more suitable than the case with UL=0.3 because, under such conditions and based on Section 2.4, the capacity of energy generation by RESs increases and energy demand of passive loads and EVs decreases.Table 4
Values of operation indices in the RT scheduling of networks for different cases and various values of UL.
CaseIIIUL=0UL=0.1UL=0.2UL=0.3Operation cost, (2a), in $4151.72861.23013.53191.23387.4Energy loss (p.u.)7.245.175.355.545.76MVD (p.u.)0.1120.0780.0790.0810.082MPD (p.u.)00.0840.0860.0870.089MTD (p.u.)0.1100.0870.0880.0890.090MOV (p.u.)00.0090.0070.0060.005MOP (p.u.)00000MOT (p.u.)00.01020.0100.0080.007


## 4.2.1. Evaluation of the Convergence of the BD Algorithm for the Proposed Scheme

The convergence plot of the problem of EH participation in DA and RT energy markets based on the BD algorithm is shown in Figure4 for different values of uncertainty level (UL). Referring to Figure 4(a), the difference between the upper and lower values of the BD algorithm zupper−zlower in the robust model of problems (1a)–(1w) and (2a)–(2u) for different UL values from 0 to 0.1 has a value less than the permissible tolerance (ε) of the BD at the sixth convergence iteration. The value of ε in this paper is set $0.5. Thus, this problem reaches the optimal point at the sixth iteration of the BD solver under the mentioned conditions. But, based on Figure 4(a), its computational time for different ULs differs. For UL=0, the computational time is 2.41 s, and for UL=0.1, it is 2.53 s. Based on Section 2.4, the increase in the UL means a reduction in the solution space of the problem in the worst-case scenario in comparison to a scenario corresponding to the deterministic model (with UL=0). In this case, the robust optimal solution is achieved at a longer computational time. Yet, it is worth noting that, based on Figure 4(a), increasing the UL to 0.2 leads to an increase in the convergence iteration number of the BD and computational time for the robust problems (1a)–(1w) and (2a)–(2u) compared to those of the deterministic model. These conditions are due to the reduced solution space of the mentioned problem compared to UL=0. Therefore, a small increase in the UL with respect to UL=0 only increases the computational time of the problem, while a large increase in the UL increases both computational time and the convergence iteration number of the BD. The mentioned conditions for the robust problem described in (1) and (2) are also true for achieving the optimal solution of the robust problem in which EHs participate in the RT energy market (3), as is shown in Figure 4(b). Nonetheless, this problem suffers from a high number of convergence iteration and longer computational time compared to problems (1a)–(1w) and (2a)–(2u) because its time step is 5 minutes; thus, 288 time steps are needed for a 24-hour operation horizon. The problem of EH participation in the DA energy market has a time step of 1 hour, so there are 24 time steps. Following this, the volume of problem (3a)–(3c) is greater than the formulation provided in (1) and (2) and its computational time and convergence iteration number are higher than those of formulations (1a)–(1w) and (2a)–(2u).Figure 4
BD convergence curve in the problem of (a) the first stage and (b) the second stage for different values of uncertainty level.
(a)(b)


## 4.2.2. Evaluation of the Capability of the Two-Layer EMS Scheme for Energy Networks with EHs

The capability of three different models of EMS on electricity, gas, and heating networks with EHs is evaluated in terms of the deterministic model (UL=0). The details are provided below.
(i)
Uncoordinated model (UM) of sources and active loads in the two-layer EMS framework: it is assumed in this section that there is no coordinated form like an EH to integrate sources and active loads (ESSs and responsive loads) and they are managed separately by their operator. Therefore, the proposed scheme (1)–(3) is solved for individual elements of the EH separately. In other words, to examine the capability of the CHP, only the model of this element is used in the formulation (1a)–(1w) and models of other elements are not employed [12].(ii)
Coordinated model (CM) of sources and active loads in a two-layer EMS framework: in this design, various sources and active loads are coordinated in the form of an EH and are managed by the EH operator. Also, operators of EHs are mutually coordinated with those of energy networks to meet the second layer of the EMS. The formulation is similar to problems (1)–(3).(iii)
Coordinated model (CM) of sources and active loads in a single-layer EMS framework: it is assumed in this design that elements of EHs are mutually coordinated with operators of energy networks. The modeling of this design in the DA operation includes the objective function (1) and constraints of the problem described by (1) and (2). It means that the objective function (2a) is removed from problems (1a)–(1w) and (2a)–(2u) and the remaining relations form a single-layer problem. Such conditions hold for problem (3a)–(3c). Since the proposed problem is single-layer and linear, CPLEX is used to solve it and the BD algorithm cannot be executed on that.The results obtained in these parts of the study are tabulated in Tables2 and 3. Based on Table 2, the total profit of sources and active loads in the UM method for DA and RT energy markets is $406.4 and $409.3. But these values for the CM method in the mentioned markets are $456.12 and $457.94, meaning that the CM design for sources and active loads in the form of an EH succeeds to enhance the profit of mentioned elements in DA and RT markets by up to 12.23% and 11.88% compared to a case they were managed separately in the UM design. The value of unbalance function (3a) between these markets is $22.9 (406.4-409.3) for the UM design, while it decreases to $21.82 in the CM design. A comparison between the CM design with single- and two-layer EMSs shows that the profit of sources and active loads in the EH form have negligible differences. The difference between the profit of EHs in the mentioned two designs for DA and RT energy markets is $0.46 (456.12-455.66) and 0.41% (457.53-457.94). Referring to Table 2, this is a small value with respect to the profit obtained in the mentioned markets. Yet, note that the computational time of the proposed problem for the single-layer EMS is roughly 73.5% (8.45−4.87/4.87) higher than that of the proposed design model based on the two-layer EMS. For the two-layer EMS case, all data are divided into two problems, such as problem (1a)–(1w) and problem (2a)–(2u); but, in the case of single-layer EMS, all data is integrated into a single problem and so the computational time will be shorter for the two-layer EMS than the single-layer model. Another advantage of the two-layer EMS compared to the single-layer EMS is that it models objectives of energy networks in the objective function of the lower-level problem, while it was removed in the single-layer EMS. Consequently, according to Table 2 and the discussions made above, the CM design for sources and active loads corresponding to the two-layer EMS is superior economically and in terms of computational time to other EMS models. Table 3 reports the financial benefits of EHs for different ULs in the case of the two-layer EMS. As per this table, the increased UL reduces the EHs’ profit in both DA and RT energy markets because, based on Section 2.4, increasing the UL will increase energy demand of passive loads and EVs as well as reduce the energy generation capacity of RESs in the worst-case scenario compared to the deterministic model (UL=0). Additionally, increasing the UL increases the value of unbalance function (3a) between the two mentioned markets. As another point, note that considering the robust model of uncertainties in the worst-case scenarios (UL=0.3) imposes 107.69% (456.12-348.43) and $107.36 (457.94-350.58) in DA and RT energy markets compared to the deterministic model. This cost is caused by considering the prediction error of uncertainty parameters.Table 2
Profit of EHs in the different energy management methods forUL=0.
Energy management methodProfit ($) of EHs in the energy marketCalculation time (s)DARTTwo-layer uncoordinated model (UM)CHP149.84149.113.41RES47.8748.553.08Boiler89.4690.233.33EVs-36.71-36.023.44EES58.6459.883.43TES51.5750.883.38DRP45.7346.673.12Sum406.40409.30—Two-layer coordinated model (CM)Total456.12457.944.87Gain compared to UM12.23%11.88%One-layer coordinated modelTotal455.66457.538.45Table 3
Unbalanced situation between DA and RT energy markets for different values of uncertainty level.
UL00.10.20.3Profit ($) in energy marketDA456.12421.76385.55348.43RT457.94423.73387.61350.58Unbalanced function, (3a), in $23.3123.8814.2444.622


## 4.2.3. Analysis of the Operation of EHs

The daily power curve for EHs participating in the markets is shown in Figure5. As per Figure 5(a), the negative value of active power generation by EHs between 1:00 and 8:00 means that EVs, batteries, and loads are charged. In parallel to this period, the electricity price is low and thus EHs benefit higher profit. Sources, ESS, and DRPs transfer power into the electricity network during 9:00-24:00 as the value of power is positive in this time interval; hence, higher benefits are gained. The power injection by EHs into the network is high from 17:00 to 00:00 in comparison to the 9:00-16:00 interval because EVs, ESSs, and loads are discharged. Note that CHPs and RESs transfer power to the electrical network during 17:00-24:00. The reason is that the price of electricity is high [51]. At 1:00-8:00 and 17:00-22:00 intervals, EHs transfer reactive power to the electricity network based on Figure 5(b). In terms of heating power, EHs inject it into the heating network during 8:00-15:00 and 23:00-24:00 or 1:00/7:00 as well as 16:00-22:00 (Figure 5(c)). Since the price of heating power is low at these intervals, the TES is charged and together with the boiler and CHP inject the heating power to the heating network at periods 8:00-15:00 and 23:00-24:00, aimed at gaining more profit for EHs. The TES is disconnected from 16:00 to 22:00, and the heating power generation by CHP and boiler is low. Gas power of EHs is based on Figure 5(d). The profit curve of EHs for individual days can be seen in Figure 6. The amount of profit relies on the prices of active power, gas power, and heating power. During 1:00-8:00, the profit obtained by EHs in the markets has a negative value because active power, gas power, and heating power are low in this interval (Figure 5). Nevertheless, profit is positive at the rest of the hours.Figure 5
Daily curve of (a) active power, (b) reactive power, (c) heating power, and (d) gas power for all EHs in the DA and RT energy markets forUL=0.2.
(a)(b)(c)(d)Figure 6
Daily EH profit curve in DA and RT energy markets forUL=0.Besides, as the figures show, the unbalance between DA and RT operations during 12:00-15:00 is higher than other hours for EH active power and their profit. The unbalance is due to the prediction error of power generation by RESs that lack power control capability. Since, based on Figure3(a), RESs inject high amounts of power into EHs, the unbalance between these operations will be higher than that in other hours. According to Figure 5(b), there is no unbalance reactive power in both DA and RT operations, and there is the same distance between DA and RT scheduling at all operation times for EH gas and heat powers based on Figures 5(c) and 5(d).


## 4.2.4. Analysis of Operation of Energy Networks

Table4 presents operation index values for electricity, gas, and heating networks for different values of ULs concerning the two-layer EMS for two case studies: power flow study (Case I) and the proposed scheme (Case II). The table also lists the operation cost of networks (2a), the energy loss of networks (the difference between energy generation and consumption), maximum voltage drop (MVD), maximum pressure drop (MPD), maximum temperature drop (MTD), maximum overvoltage (MOV), maximum overpressure (MOP), and maximum overtemperature (MOT) in the case of RT operation. Based on the reported results in this paper, it is observed that the proposed scheme in the worst-case scenario of uncertainties with UL=0.3 succeeds to reduce the operating cost of networks by almost 18.4% (4151.7−3387.4/4151.7) compared to that in Case I. Also, it helps to degrade energy loss, MVD, and MTD values by roughly 20.5%, 26.8%, and 18.2% in comparison to Case I. However, these conditions are realized for increasing MOV, MOT, and MPD by about 0.005, 0.007, and 0.089 p.u. These values of parameters are smaller than their allowed limit of 0.1 p.u. (1.1-1). Moreover, the status of mentioned indices for the case the UL decreases becomes more suitable than the case with UL=0.3 because, under such conditions and based on Section 2.4, the capacity of energy generation by RESs increases and energy demand of passive loads and EVs decreases.Table 4
Values of operation indices in the RT scheduling of networks for different cases and various values of UL.
CaseIIIUL=0UL=0.1UL=0.2UL=0.3Operation cost, (2a), in $4151.72861.23013.53191.23387.4Energy loss (p.u.)7.245.175.355.545.76MVD (p.u.)0.1120.0780.0790.0810.082MPD (p.u.)00.0840.0860.0870.089MTD (p.u.)0.1100.0870.0880.0890.090MOV (p.u.)00.0090.0070.0060.005MOP (p.u.)00000MOT (p.u.)00.01020.0100.0080.007


## 5. Conclusion

This paper presented the optimal participation of grid-connected EHs in DA and RT energy markets. The problem is composed of two stages that deal with the participation of EHs in hourly DA and 5 min RT markets. Formulation of the first stage is a bilevel optimization; the upper level of which maximizes the EHs’ profit in the DA market constrained to the operation model of sources and active loads in the form of an EH. The lower level addresses the minimization of the energy cost of electricity, gas, and heating networks subject to optimized-linearized power flow equations. Formulation of the second stage is similar to that of the first stage, except for its time step which is 5 minutes and the objective function of its upper-level problem is attempting to minimize the unbalance between EHs’ profits in the markets. Then, robust programming is used to model the uncertainties of load, energy price, renewable energy, and energy demand of EVs. In the end, the BD algorithm is adopted to solve the problems of the first and second stages. Numerical results show that the BD algorithm can find the optimal robust solution at the minimum number of convergence iterations and the shortest computational time. The coordinated model of sources and active loads in the EH form obtains more financial benefit than the independent management in the mentioned markets. The proposed scheme in the form of a two-layer EMS demands less computational time than the single-layer EMS model. Moreover, by achieving an optimal DA and RT scheduling for EHs, the suggested scheme leads to a lower unbalance between DA and RT operations at different uncertainty levels. It can also enhance various operation indices of different networks even in the worst-case scenario with a high uncertainty level (0.3). Under such conditions, with the appearance of permissible values for pressure deviation, overvoltage, and overtemperature, it helps to enhance energy cost, energy loss, maximum voltage drop, and maximum temperature drop by 18.4%, 20.5%, 26.8%, and 18.2% compared to power flow studies.Therefore, the advantages of the proposed scheme compared to the previous research are as follows:(i)
This scheme used two-layer energy management so that the energy management of EHs was done in the first layer, and the energy management of networks was included in the second layer. This case made it possible to consider the goals of EHs and the network operator in the optimization problem, and it also reduced the amount of information sent to the network operator(ii)
In this scheme, modeling of both DA and RT energy markets is presented. This case makes the operation error between two markets to be measured, and it is possible to reduce the error between them by defining special penalty functions(iii)
The proposed scheme is able to improve the economic and operation situation of the network at the same time with the optimal energy management of EHs, and, in addition to this, it also extracts the optimal economic situation for EHs(iv)
BD algorithm is capable of extracting the optimal solution in low computation time. Also, robust modeling reduces the volume of the problem, even though the number of uncertainty parameters of the proposed design is high. These two cases solve the problem of low computation time. From there, the implementation step is low in the market, so low computation time is of special importance. This is observed in the proposed schemeIn the end, based on the stated advantages, the problems in the field of research have been resolved in the proposed scheme.



---

*Source: 5087872-2023-03-23.xml*