# Participation of Grid-Connected Energy Hubs and Energy Distribution Companies in the Day-Ahead Energy Wholesale and Retail Markets Constrained to Network Operation Indices

Authors: Omid Kohansal; Mahmoud Zadehbagheri; Mohammad Javad Kiani; Samad Nejatian
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

In this paper, the optimal scheduling of energy grids and networked energy hubs based on their participation in the day-ahead energy wholesale and retail markets is presented. The problem is formulated as a bilevel model. Its upper level minimizes the expected energy cost of electricity, gas, and heating grids, especially in the form of private distribution companies in the mentioned markets, in the first objective function, and it minimizes the expected energy loss of these networks in the second objective function. This problem is constrained by linearized optimal power flow equations. The lower-level formulation minimizes the expected energy cost of hubs (equal to the difference between sell and purchase of energy) as an objective function in the retail market. Constraints of this model are the operation formulation of sources and active loads and the flexibility limit of hubs. The unscented transformation approach models the uncertainties of load, renewable power, energy price, and energy demand of mobile storage. Then, the Karush–Kuhn–Tucker approach and Pareto optimization technique based onε-constraint are adopted to extract the single-level single-objective formulation. Finally, obtained results verify the capability of the present method in improving the economic status of hubs and the economic and operation situation of the mentioned networks simultaneously so that the proposed scheme by managing the power of energy hubs compared with power flow studies has been able to reduce operating costs by 8%, reduce energy losses by 10%, and improve voltage profile and temperature by 36% and 30%.


---

## Body

## 1. Introduction

### 1.1. Motivation

Nowadays, due to the advancement in power generation technologies (such as combined heat and power (CHP) systems and renewable energy sources (RES)) and energy storage systems (ESSs) (such as electric vehicles (EVs)), energy management programs such as environmentally friendly demand response programs (DRP) are adopted to save energy. Thus, these power sources and active loads (ALs) can enhance technical status of the network and reduce its cost by participating in the power system [1]. For example, RESs in the electricity grid reduce the network operating cost and energy price, which is commensurate with the increase in social welfare [2]. In addition, the increase in the number of these elements in energy networks (ENs), such as electricity, gas, and heating networks, raises the volume of data sent to the energy network operator (ENO), which complicates decision making. To address this, smart grid theory suggests that different sources and ALs be managed in the form of different aggregators such as energy hub (EH), virtual power plant (VPP), or microgrid (MG) [3]. However, since EH can manage several types of energy simultaneously, it will have higher energy efficiency than VPP and MG because VPP and MG generally produce only electrical energy [3]. In order to establish proper energy management in energy networks, sources and ALs will have bilateral coordination with the EH operator. EH operators also coordinate with ENOs, which is defined as a two-layer energy management system (EMS) [4]. Therefore, ENOs only need the information received by the EH operator, which is predicted to have a higher processing speed for ENOs than if the sources and ALs were directly related to the ENO. It is worth noting that EHs are able to inject energy into various networks within operation hours. As a result, by participating in the energy market, they can receive a good financial benefit from the market [5]. Each EN can act as a private distribution company (DisCo), purchasing energy from the upstream grid and delivering it to consumers and connected EHs. Finally, it is anticipated that ENs will also be able to enjoy adequate financial sources in this strategy. Therefore, by establishing an appropriate operating framework for ENs in the presence of EHs, steps can be taken to improve various technical and economic indices.


### 1.2. Literature Review

Various researches and works concerning energy management and optimal operation of EHs have been carried out in different ENs. In Reference [6], the stable operation of EHs containing EVs and CHP is presented. The plan considers minimizing the operating cost of EHs in electricity, gas, and heat networks as an objective function. It is also bound to the power flow equations in different networks and the operating constraints of these networks in addition to the operating model of these EHs. It obtains a linear approximation model for the problem in order to achieve the optimal solution as fast as possible. In the proposed approach, the voltage magnitude and the relationship of the gas passing through the pipeline are linearized using the linear piecewise method. Also, the capacity limit of the distribution lines, which is a circular plane, is approximated to a regular polygon. Finally, by implementing it on a test system, linear approximation model has a small computational error compared with the real problem model, as shown in the obtained results. Also, the optimal performance of CHPs and EVs in the form of EH has been able to reduce energy losses and achieve smoother voltage, pressure, and temperature profiles. In Reference [7], EH consists of RES, DRP, and ice energy storage. In this reference, the effects of the mentioned storage on the operation status of EH have been evaluated so that it has been able to peak shave in EH by storing energy during off-peak period and discharging it during peak period. In Reference [8], the stable operation of EHs connected to electrical and gas MGs is also modeled. It uses adaptive optimization to model uncertainties of power generation and loads, which ultimately results in a min-max-min problem for the design. It then uses the Benders decomposition method to achieve the optimal solution. The suggested method is divided into a master problem and a subproblem. In the master problem, the unit commitment model is solved and, in the subproblem, the operation of MGs with EHs is solved. In Reference [9], the stochastic operation of MGs with EHs including photovoltaic panel, compressed air storage, and ice is presented. According to Reference [9], the optimal operation of EHs reduces energy losses and energy costs in MGs. Reference [10] presents the optimal operation of grid-connected EHs subject to clean environmental conditions. According to Reference [10], energy hubs using RES such as solar panels and wind systems can, in addition to improving operation indices, play a very effective role in providing clean energy. In Reference [11], the capability of EHs in system reliability enhancement in an N − 1 event has been investigated. Based on the obtained results, due to its location in the consumption points, the EH is able to significantly reduce the interruption rate in the mentioned event. This is commensurate with the high ability of the EH to improve system reliability.Many studies have been presented in the field of modeling EH participation in the market. In Reference [12, 13], the participation of the electrical sector of the EH in the energy market is expressed using the bilateral contract and power pool models. In Reference [12], a deterministic model is presented for the proposed design, but in Reference [13], a robust model is incorporated for modeling uncertainties related to RES price, load, and generation capacity of RES. Based on the results obtained from Reference [12, 13], EH can manage the sources and get a good financial benefit by proper scheduling of sources and ALs even in the worst-case scenario caused by the aforementioned uncertainties. Reference [14] discusses the participation of grid-connected EHs in the electricity, gas, and heating energy market. In the proposed method, EHs generally consume gas energy due to having CHP and boiler, so they pay the gas consumption cost to the gas energy market. But they generate energy in the electricity and heating markets most of the time, so EHs benefit financially from these two markets. In Reference [12–14], energy price is considered as a parameter, but we note that energy price is different for different parts of the power system. Therefore, it is a variable that can be deduced from the market clearing price (MCP) problem. This has been implemented in [15] in the electricity and heating energy markets. In Reference [16], as in Reference [15], the optimal bidding model of EHs in the energy market is presented. Stochastic dynamic programming (SDP) is also used in Reference [17] for the optimal operation of the clearing EH.Operation of an EH that includes various types of energy resources is optimally configured in [18] using a stochastic framework. Random scenarios are produced using the historical data so that uncertain quantities like energy price in the market, speed of the wind, and irradiance are appropriately modeled. Optimal management of electrical and natural gas energy is addressed in [19] so that the cost is reduced. Power-to-gas (P2G) technology is employed in the energy hub. By doing this, natural gas is generated and sold to the gas network and, as a result, the congestion in gas transmission lines is resolved and the cost of the energy hub is significantly decreased. To make changes in the peak load periods, a demand response program is utilized in the suggested model. Furthermore, in an attempt to meet the operator’s expectations, several energy storage equipment and distributed generation are adopted. Uncertainties of wind turbine output power are modeled with the help of a scenario generation approach. Mixed-integer linear programming is also used to mathematically express the problem. Reference [20] presents a cost-based optimization model of an energy hub with storage equipment while taking into account resilience limits and some other constraints. The method divides the loads into either critical or noncritical loads. According to this method, critical loads have priority to be supplied and noncritical loads should not experience frequent interruptions. The role of storage equipment in such a model is essential. Reference [21] presents a novel arrangement for a multiresource system. In this paper, the waste heat is recovered in a regenerative gas turbine cycle to run a district heating heat exchanger, a Rankine cycle, an ejector refrigeration cycle, and a proton exchange membrane electrolyzer cycle. The authors in [22] provide a stochastic scheduling plan for multiresource hub systems so that energy hubs and distribution systems are both optimized at the same time. RES, DRPs, uncertainties, and pollution emission levels are taken into account in the given plan so that the reliable operation of the whole system is guaranteed. Moreover, probability density functions are adopted to present proper models of uncertainties including electricity, heating, and cooling demand; price of electricity; and power output of renewables. Table 1 lists some of the recent works conducted in this topic.Table 1
A review on the recent literature.
Ref.Energy market modelUncertainty modelFlexibility model[6]NoRobustNo[7]NoSBSPNo[8]NoRobustNo[9–11]NoSBSPNo[12]RetailDeterministicNo[13]RetailRobustNo[14–19]RetailSBSPNo[20]NoSBSPNo[21]NoRobustNo[22]NoSBSPNoProposed modelWholesale and retailUT methodYesSBSP : scenario-based stochastic programming; UT : unscented transformation.


### 1.3. Research Gaps

Based on the backgrounds mentioned in the previous subsection and Table1, the major research gaps for energy management of energy networks with EHs include the following:(i)
Most studies such as [6–11, 20–22] evaluated the operation of EHs and fewer have talked about their participation in the energy market. Also, some of the researchers that have considered the energy market model for this model have generally considered this model either for energy networks or have used it for EHs, as in Reference [12–19]. Note, however, that in a competitive environment that leads to financial gain for its participants and increased social welfare, it is appropriate to evaluate the simultaneous performance of DisCos and EHs in the energy market. Therefore, it is expected that all the participants of the proposed scheme, in addition to improving the operation indices of various networks, will also receive suitable financial benefits from the energy market.(ii)
The electrical part of the EH has low flexibility in the presence of RES [1]. Also, its heating part will have low flexibility in the presence of CHPs [23]. Low flexibility causes the results of real-time and day-ahead operations to be the same, which will lead to an imbalance of demand and supply for EHs in real-time operation [24]. Following this, various studies suggest incorporating flexibility sources like DRP, ESS, and nonrenewable energy sources in EH. But in less research, EH flexibility modeling has been proposed. Note, however, that in order to determine the status of an index for an element, its numerical results are required, which are available from mathematical modeling of the considered index.(iii)
In most studies such as References [7, 9–11, 14–20, 22], nonparametric methods such as scenario-based stochastic programming (SBSP) have been employed to provide modelings for uncertainties. In this scheme, access to the guaranteed optimal solution needs many scenarios, which increases the calculation period, which is of particular importance in operation problems. So, it is appropriate to incorporate parametric methods such as UT to model uncertainties. In these methods, the number of required scenarios is dependent on the number of uncertainty parameters and generally has the lowest number of scenarios. Of course, in some cases [6, 8, 13, 21], robust modeling has been used, which has only one scenario called the “worst-case scenario.” Although it has a low computational time, it is necessary to consider several different scenarios to determine the exact value of some indices such as flexibility and reliability [1]. Therefore, robust modeling of uncertainties will not be a popular technique in these conditions.


### 1.4. Contributions

Noting the research gaps early mentioned, the current study proposes wholesale and retail energy markets modeling for different energy networks such as electricity, gas, and heating networks in the presence of EHs in accordance with the two-layer EMS. In the first layer of EMS, the coordination of sources and ALs with the EH operator (EHO) is considered. In the second layer of EMS, the coordination of the EHOs with ENOs is taken into account. To cover the first research gap, in this study, these networks participate in the wholesale market in the form of a private DisCo and purchase energy from it. It also shares the energy purchased in the retail energy market environment between consumers and connected EHs. Therefore, this plan is a bilevel optimization problem considered in the upper-level problem of the operation model of energy networks. The model of operation of grid-connected EHs is formulated in the lower-level problem. In the upper-level problem, it is assumed that ENOs aim to enhance the economic situation of DisCos and operation indices of energy networks. It is therefore formulated in the form of two-objective optimization. One objective function finds the minimum expected cost of DisCos in the mentioned energy markets, and the other objective function improves operation indices. Operation of networks means minimization of expected energy losses in the mentioned networks. This problem is bounded by the linearized optimal power flow (LOPF) constraints in these networks. To provide an integrated single-objective model, theε-constraint-based Pareto optimization method is adopted. The lower-level problem considers the maximization of expected profit of EHs in the retail energy market, which is subject to the operation model of sources and ALs, and the flexibility limit of EHs to eliminate the second research gap. The Karush–Kuhn–Tucker (KKT) method is then used to develop the single-level model. In order to cover the third research gap, parametric methods such as the unscented transformation (UT) method are used to model load uncertainties, energy prices, the active power of RES, and the energy demand of EVs. The novelties of the current study include the following:(i)
Modeling the wholesale and retail energy market for different energy networks as DisCo in the presence of EHs(ii)
Modeling energy flexibility considering uncontrollable sources and flexibility sources(iii)
Gaining financial benefit for DisCos and EHs simultaneously(iv)
Simultaneous modeling of operation, economic, and flexibility indices in the two-layer energy management plan of EHs in different energy networks(v)
Stochastic modeling of uncertainties related to load, energy prices, RESs power, and EVs energy demand using a parametric method, that is, UT


### 1.5. Paper Organization

Section2 develops the mathematical model of the proposed scheme and models the uncertainties. Section 3 extracts the integrated model of the proposed design, which is extracted as a single-objective single-level model. Eventually, Section 4 provides numerical results and Section 5 presents conclusions.


## 1.1. Motivation

Nowadays, due to the advancement in power generation technologies (such as combined heat and power (CHP) systems and renewable energy sources (RES)) and energy storage systems (ESSs) (such as electric vehicles (EVs)), energy management programs such as environmentally friendly demand response programs (DRP) are adopted to save energy. Thus, these power sources and active loads (ALs) can enhance technical status of the network and reduce its cost by participating in the power system [1]. For example, RESs in the electricity grid reduce the network operating cost and energy price, which is commensurate with the increase in social welfare [2]. In addition, the increase in the number of these elements in energy networks (ENs), such as electricity, gas, and heating networks, raises the volume of data sent to the energy network operator (ENO), which complicates decision making. To address this, smart grid theory suggests that different sources and ALs be managed in the form of different aggregators such as energy hub (EH), virtual power plant (VPP), or microgrid (MG) [3]. However, since EH can manage several types of energy simultaneously, it will have higher energy efficiency than VPP and MG because VPP and MG generally produce only electrical energy [3]. In order to establish proper energy management in energy networks, sources and ALs will have bilateral coordination with the EH operator. EH operators also coordinate with ENOs, which is defined as a two-layer energy management system (EMS) [4]. Therefore, ENOs only need the information received by the EH operator, which is predicted to have a higher processing speed for ENOs than if the sources and ALs were directly related to the ENO. It is worth noting that EHs are able to inject energy into various networks within operation hours. As a result, by participating in the energy market, they can receive a good financial benefit from the market [5]. Each EN can act as a private distribution company (DisCo), purchasing energy from the upstream grid and delivering it to consumers and connected EHs. Finally, it is anticipated that ENs will also be able to enjoy adequate financial sources in this strategy. Therefore, by establishing an appropriate operating framework for ENs in the presence of EHs, steps can be taken to improve various technical and economic indices.


## 1.2. Literature Review

Various researches and works concerning energy management and optimal operation of EHs have been carried out in different ENs. In Reference [6], the stable operation of EHs containing EVs and CHP is presented. The plan considers minimizing the operating cost of EHs in electricity, gas, and heat networks as an objective function. It is also bound to the power flow equations in different networks and the operating constraints of these networks in addition to the operating model of these EHs. It obtains a linear approximation model for the problem in order to achieve the optimal solution as fast as possible. In the proposed approach, the voltage magnitude and the relationship of the gas passing through the pipeline are linearized using the linear piecewise method. Also, the capacity limit of the distribution lines, which is a circular plane, is approximated to a regular polygon. Finally, by implementing it on a test system, linear approximation model has a small computational error compared with the real problem model, as shown in the obtained results. Also, the optimal performance of CHPs and EVs in the form of EH has been able to reduce energy losses and achieve smoother voltage, pressure, and temperature profiles. In Reference [7], EH consists of RES, DRP, and ice energy storage. In this reference, the effects of the mentioned storage on the operation status of EH have been evaluated so that it has been able to peak shave in EH by storing energy during off-peak period and discharging it during peak period. In Reference [8], the stable operation of EHs connected to electrical and gas MGs is also modeled. It uses adaptive optimization to model uncertainties of power generation and loads, which ultimately results in a min-max-min problem for the design. It then uses the Benders decomposition method to achieve the optimal solution. The suggested method is divided into a master problem and a subproblem. In the master problem, the unit commitment model is solved and, in the subproblem, the operation of MGs with EHs is solved. In Reference [9], the stochastic operation of MGs with EHs including photovoltaic panel, compressed air storage, and ice is presented. According to Reference [9], the optimal operation of EHs reduces energy losses and energy costs in MGs. Reference [10] presents the optimal operation of grid-connected EHs subject to clean environmental conditions. According to Reference [10], energy hubs using RES such as solar panels and wind systems can, in addition to improving operation indices, play a very effective role in providing clean energy. In Reference [11], the capability of EHs in system reliability enhancement in an N − 1 event has been investigated. Based on the obtained results, due to its location in the consumption points, the EH is able to significantly reduce the interruption rate in the mentioned event. This is commensurate with the high ability of the EH to improve system reliability.Many studies have been presented in the field of modeling EH participation in the market. In Reference [12, 13], the participation of the electrical sector of the EH in the energy market is expressed using the bilateral contract and power pool models. In Reference [12], a deterministic model is presented for the proposed design, but in Reference [13], a robust model is incorporated for modeling uncertainties related to RES price, load, and generation capacity of RES. Based on the results obtained from Reference [12, 13], EH can manage the sources and get a good financial benefit by proper scheduling of sources and ALs even in the worst-case scenario caused by the aforementioned uncertainties. Reference [14] discusses the participation of grid-connected EHs in the electricity, gas, and heating energy market. In the proposed method, EHs generally consume gas energy due to having CHP and boiler, so they pay the gas consumption cost to the gas energy market. But they generate energy in the electricity and heating markets most of the time, so EHs benefit financially from these two markets. In Reference [12–14], energy price is considered as a parameter, but we note that energy price is different for different parts of the power system. Therefore, it is a variable that can be deduced from the market clearing price (MCP) problem. This has been implemented in [15] in the electricity and heating energy markets. In Reference [16], as in Reference [15], the optimal bidding model of EHs in the energy market is presented. Stochastic dynamic programming (SDP) is also used in Reference [17] for the optimal operation of the clearing EH.Operation of an EH that includes various types of energy resources is optimally configured in [18] using a stochastic framework. Random scenarios are produced using the historical data so that uncertain quantities like energy price in the market, speed of the wind, and irradiance are appropriately modeled. Optimal management of electrical and natural gas energy is addressed in [19] so that the cost is reduced. Power-to-gas (P2G) technology is employed in the energy hub. By doing this, natural gas is generated and sold to the gas network and, as a result, the congestion in gas transmission lines is resolved and the cost of the energy hub is significantly decreased. To make changes in the peak load periods, a demand response program is utilized in the suggested model. Furthermore, in an attempt to meet the operator’s expectations, several energy storage equipment and distributed generation are adopted. Uncertainties of wind turbine output power are modeled with the help of a scenario generation approach. Mixed-integer linear programming is also used to mathematically express the problem. Reference [20] presents a cost-based optimization model of an energy hub with storage equipment while taking into account resilience limits and some other constraints. The method divides the loads into either critical or noncritical loads. According to this method, critical loads have priority to be supplied and noncritical loads should not experience frequent interruptions. The role of storage equipment in such a model is essential. Reference [21] presents a novel arrangement for a multiresource system. In this paper, the waste heat is recovered in a regenerative gas turbine cycle to run a district heating heat exchanger, a Rankine cycle, an ejector refrigeration cycle, and a proton exchange membrane electrolyzer cycle. The authors in [22] provide a stochastic scheduling plan for multiresource hub systems so that energy hubs and distribution systems are both optimized at the same time. RES, DRPs, uncertainties, and pollution emission levels are taken into account in the given plan so that the reliable operation of the whole system is guaranteed. Moreover, probability density functions are adopted to present proper models of uncertainties including electricity, heating, and cooling demand; price of electricity; and power output of renewables. Table 1 lists some of the recent works conducted in this topic.Table 1
A review on the recent literature.
Ref.Energy market modelUncertainty modelFlexibility model[6]NoRobustNo[7]NoSBSPNo[8]NoRobustNo[9–11]NoSBSPNo[12]RetailDeterministicNo[13]RetailRobustNo[14–19]RetailSBSPNo[20]NoSBSPNo[21]NoRobustNo[22]NoSBSPNoProposed modelWholesale and retailUT methodYesSBSP : scenario-based stochastic programming; UT : unscented transformation.


## 1.3. Research Gaps

Based on the backgrounds mentioned in the previous subsection and Table1, the major research gaps for energy management of energy networks with EHs include the following:(i)
Most studies such as [6–11, 20–22] evaluated the operation of EHs and fewer have talked about their participation in the energy market. Also, some of the researchers that have considered the energy market model for this model have generally considered this model either for energy networks or have used it for EHs, as in Reference [12–19]. Note, however, that in a competitive environment that leads to financial gain for its participants and increased social welfare, it is appropriate to evaluate the simultaneous performance of DisCos and EHs in the energy market. Therefore, it is expected that all the participants of the proposed scheme, in addition to improving the operation indices of various networks, will also receive suitable financial benefits from the energy market.(ii)
The electrical part of the EH has low flexibility in the presence of RES [1]. Also, its heating part will have low flexibility in the presence of CHPs [23]. Low flexibility causes the results of real-time and day-ahead operations to be the same, which will lead to an imbalance of demand and supply for EHs in real-time operation [24]. Following this, various studies suggest incorporating flexibility sources like DRP, ESS, and nonrenewable energy sources in EH. But in less research, EH flexibility modeling has been proposed. Note, however, that in order to determine the status of an index for an element, its numerical results are required, which are available from mathematical modeling of the considered index.(iii)
In most studies such as References [7, 9–11, 14–20, 22], nonparametric methods such as scenario-based stochastic programming (SBSP) have been employed to provide modelings for uncertainties. In this scheme, access to the guaranteed optimal solution needs many scenarios, which increases the calculation period, which is of particular importance in operation problems. So, it is appropriate to incorporate parametric methods such as UT to model uncertainties. In these methods, the number of required scenarios is dependent on the number of uncertainty parameters and generally has the lowest number of scenarios. Of course, in some cases [6, 8, 13, 21], robust modeling has been used, which has only one scenario called the “worst-case scenario.” Although it has a low computational time, it is necessary to consider several different scenarios to determine the exact value of some indices such as flexibility and reliability [1]. Therefore, robust modeling of uncertainties will not be a popular technique in these conditions.


## 1.4. Contributions

Noting the research gaps early mentioned, the current study proposes wholesale and retail energy markets modeling for different energy networks such as electricity, gas, and heating networks in the presence of EHs in accordance with the two-layer EMS. In the first layer of EMS, the coordination of sources and ALs with the EH operator (EHO) is considered. In the second layer of EMS, the coordination of the EHOs with ENOs is taken into account. To cover the first research gap, in this study, these networks participate in the wholesale market in the form of a private DisCo and purchase energy from it. It also shares the energy purchased in the retail energy market environment between consumers and connected EHs. Therefore, this plan is a bilevel optimization problem considered in the upper-level problem of the operation model of energy networks. The model of operation of grid-connected EHs is formulated in the lower-level problem. In the upper-level problem, it is assumed that ENOs aim to enhance the economic situation of DisCos and operation indices of energy networks. It is therefore formulated in the form of two-objective optimization. One objective function finds the minimum expected cost of DisCos in the mentioned energy markets, and the other objective function improves operation indices. Operation of networks means minimization of expected energy losses in the mentioned networks. This problem is bounded by the linearized optimal power flow (LOPF) constraints in these networks. To provide an integrated single-objective model, theε-constraint-based Pareto optimization method is adopted. The lower-level problem considers the maximization of expected profit of EHs in the retail energy market, which is subject to the operation model of sources and ALs, and the flexibility limit of EHs to eliminate the second research gap. The Karush–Kuhn–Tucker (KKT) method is then used to develop the single-level model. In order to cover the third research gap, parametric methods such as the unscented transformation (UT) method are used to model load uncertainties, energy prices, the active power of RES, and the energy demand of EVs. The novelties of the current study include the following:(i)
Modeling the wholesale and retail energy market for different energy networks as DisCo in the presence of EHs(ii)
Modeling energy flexibility considering uncontrollable sources and flexibility sources(iii)
Gaining financial benefit for DisCos and EHs simultaneously(iv)
Simultaneous modeling of operation, economic, and flexibility indices in the two-layer energy management plan of EHs in different energy networks(v)
Stochastic modeling of uncertainties related to load, energy prices, RESs power, and EVs energy demand using a parametric method, that is, UT


## 1.5. Paper Organization

Section2 develops the mathematical model of the proposed scheme and models the uncertainties. Section 3 extracts the integrated model of the proposed design, which is extracted as a single-objective single-level model. Eventually, Section 4 provides numerical results and Section 5 presents conclusions.


## 2. Model of the Proposed Problem

### 2.1. Mathematical Formulation

The two-layer energy management of EHs in electricity, gas, and heat networks is described here. The first layer deals with establishing coordination among EHOs, sources, and ALs. The second layer coordinates EHOs and ENOs. As a result, we have an optimization problem with two upper and lower levels. The former is related to the participation of ENs in retail and wholesale markets, while the latter manages the power supply and demand of the second layer of the energy management plan. The suggested plan is mathematically formulated:(1)minCostENs=∑s∈ΠSπs∑t∈ΠOHλEWt,sPESr,t,s+λHWt,sHHSr,t,s+λGWt,sGGSr,t,s+∑s∈ΠSπs∑t∈ΠOH∑i∈ΠEHλEWt,sPEHi,t,s++λHWt,sHEHi,t,s++λGWt,sGEHi,t,s+,(2)minEEL=∑s∈ΠSπs∑t∈ΠOHPESr,t,s+∑i∈ΠEHPEHi,t,s+−PEHi,t,s−−∑e∈ΠENPDe,t,s+∑s∈ΠSπs∑t∈ΠOHHHSr,t,s+∑i∈ΠEHHEHi,t,s+−HEHi,t,s−−∑h∈ΠHNHDh,t,s+∑s∈ΠSπs∑t∈ΠOHGGSr,t,s+∑i∈ΠEHGEHi,t,s+−GEHi,t,s−−∑g∈ΠGNGDg,t,s,which is subjected to(3)PESe,t,s−PDe,t,s+∑i∈ΠEHAEe,iPEHi,t,s+−PEHi,t,s−=∑j∈ΠENBEe,jPLe,j,t,s,∀e,t,s,(4)QESe,t,s−QDe,t,s+∑i∈ΠEHAEe,iQEHi,t,s=∑j∈ΠENBEe,jQLe,j,t,s,∀e,t,s,(5)PLe,j,t,s=GLe,jΔVe,t,s−ΔVj,t,s−BLe,jφe,t,s−φj,t,s,∀e,j,t,s,(6)QLe,j,t,s=−BLe,jΔVe,t,s−ΔVj,t,s−GLe,jφe,t,s−φj,t,s,∀e,j,t,s,(7)ΔVe,t,s,φe,t,s=0,∀e=r,t,s,(8)GGSg,t,s−GDg,t,s+∑i∈ΠEHAGg,iGEHi,t,s+−GEHi,t,s−=∑j∈ΠGNBGg,jGLg,j,t,s,∀g,t,s,(9)GLg,j,t,s=κg,jsignρg,t,s,ρj,t,sϑg,j,t,s,(10)Δϑg,j,t,s,p=signρg,t,s,ρj,t,sspρspϑΔρg,t,s,p−Δρj,t,s,p,∀g,j,t,s,l,(11)ϑg,j,t,s=∑p∈ΠPΔϑg,j,t,s,p,∀g,j,t,s(12)ρg,t,s=ρ¯g+∑p∈ΠPΔρg,t,s,p,∀g,t,s,(13)HESh,t,s−HDh,t,s+∑i∈ΠEHAHh,iHEHi,t,s+−HEHi,t,s−=∑j∈ΠENBHe,jHLh,j,t,s,∀h,t,s,(14)ρg,t,s=ρ¯g+∑p∈ΠPΔρg,t,s,p,∀g,t,s,(15)V¯e−1≤ΔVe,t,s≤V¯e−1,∀e,t,s,(16)PLe,j,t,scosk⋅Δβ+QLe,j,t,ssink⋅Δβ≤S¯Le,j∀e,j,t,s,k∈ΠK=1,2,…,nk,Δβ=360nk,(17)PESe,t,scosk⋅Δβ+QESe,t,ssink⋅Δβ≤S¯ESe,∀e=r,t,s,k,(18)ρ¯g≤ρg,t,s≤ρ¯g,∀g,t,s,(19)−G¯Lg,j≤GLg,j,t,s≤G¯Lg,j,∀g,j,t,s,(20)−G¯GSg≤GGSg,t,s≤G¯GSgm,∀g=r,t,s,(21)T¯h≤Th,t≤Th¯,∀h,t,(22)−H¯Lh,j≤HLh,j,t,s≤H¯Lh,j,∀h,j,t,s,(23)−H¯HSh≤HHSh,t,s≤H¯HSh,∀h=r,t,s(24)−H¯Lh,j≤HLh,j,t,s≤H¯Lh,j,∀h,j,t,s,(25)−H¯HSh≤HHSh,t,s≤H¯HSh,∀h=r,t,s,(26)PEH+,PEH−,QEH,HEH+,HEH−,GEH+,GEH−∈argminCostEHs=∑s∈ΠSπs∑t∈ΠOH∑i∈ΠEHλERt,sPEHi,t,s−+λHRt,sHEHi,t,s−+λGRt,sGEHi,t,s−−λEWt,sPEHi,t,s++λHWt,sHEHi,t,s++λGWt,sGEHi,t,s+,which is subjected to(27)PEHi,t,s+−PEHi,t,s−=PCi,t,s+PRi,t,s+PDRi,t,s+PDIS,t,s−PCHi,t,s−PDi,t,s,∀PEHi,t,s+,PEHi,t,s−≥0&i,t,s,(28)QEHi,t,s=QCi,t,s+QRi,t,s+QEi,t,s−QDi,t,s,∀i,t,s,(29)HEHi,t,s+−HEHi,t,s−=HCi,t,s+HBi,t,s+HDRi,t,s−HDi,t,s,∀HEHi,t,s+,HEHi,t,s−≥0&i,t,s,(30)GEHi,t,s+−GEHi,t,s−=GDRi,t,s−GCi,t,s−GBi,t,s−GDi,t,s,∀GEHi,t,s+,GEHi,t,s−≥0,&i,t,s,(31)HCi,t,s=PCi,t,s1−ηT−ηLηT⋅ηH,∀i,t,s,(32)GCi,t,s=1ηT.PCi,t,s,∀i,t,s,(33)PCi,t,scosk⋅Δβ+QCi,t,scosk⋅Δβ≤S¯Ci,∀i,t,s,k,(34)0≤HCi,t,s≤H¯Ci,∀i,t,s,(35)GBi,t,s=1ηB⋅HBi,t,s,∀i,t,s,(36)0≤HBi,t,s≤H¯Bi,∀i,t,s,(37)0≤PCHi,t,s≤CRi,∀i,t,s,(38)0≤PDISi,t,s≤DRi,∀i,t,s,(39)E¯i≤EIi+∑τ=1tηCHPCHi,τ,s−1ηDISPDISi,τ,s≤E¯i,∀i,t,s,(40)PDISi,t,s−PCHi,t,scosk⋅Δβ+QEi,t,scosk⋅Δβ≤S¯Ei,∀i,t,s,k,(41)PRi,t,scosk⋅Δβ+QRi,t,scosk⋅Δβ≤S¯Ri,∀i,t,s,k,(42)−υiPDi,t,sHDi,t,sorGDi,t,s≤PDRi,t,sHDRi,t,sorGDRi,t,s≤υiPDi,t,sHDi,t,sorGDi,t,s,∀i,t,s,(43)∑t∈ΠOHPDRi,t,sHDRi,t,sorGDRi,t,s=0,∀i,s,(44)−εF≤PEHi,t,s−PEHi,t,s=1≤εF,∀i,t,s,(45)−εF≤HEHi,t,s−HEHi,t,s=1≤εF,∀i,t,s.Upper-level problem: Eqs. (1)–(23) show the upper-level problem as a two-objective optimization. Equation (1) minimizes the expected energy cost of ENs (as a function of CostENs), where CostENs represents the sum of power purchased by ENs from upstream networks (the first term in Eq. (1)) [12] and EHs (the second term in Eq. (1)) [14]. Two energy markets, that is, wholesale and retail energy markets, are considered in the proposed plan, Figure 1. In the former, the wholesale market supplies energy to ENs. The first term in Eq. (1) models the wholesale market, where the positive value of PDS (HHS or GGS) means that energy is purchased from the market, while its negative value means that the wholesale market purchases energy from the network [14]. Additionally, EHs (passive customers) purchase or sell (purchase) energy through participating in a retail market. It is noteworthy that EHs are assumed to buy energy from ENs at the retail market price, λR, which is generally higher than the wholesale market price, λW. Nonetheless, if EHs can inject energy into ENs, ENs purchase this energy at a price of λW, which is generally less than λR. This is due to maintaining the favorable financial benefit of ENs in energy markets and the participation of EHs in energy production. This is due to the fact that a wholesaler generally buys his goods at a price lower than his sale price to a retailer. Note that based on Figure 1, two models of the wholesale and retail market are considered in this paper. It is assumed that the distribution company (energy network) purchases energy from the wholesale market. It then sells the energy it buys within itself as a retail market to its customers, that is, consumers and energy hubs. To encourage its customers in network energy management, it also allows its customers to sell energy in the retail market. Therefore, it can be stated that a distribution company is a wholesale buyer who provides its energy to customers in the form of retail sales. Hence, in the second part of Eq. (1), the price λW is used to buy energy by ENs from EHs in the retail market. Equation (2) deals with the operation of ENOs and minimizes the total expected energy loss in all energy networks (function of EEL) [23].Figure 1
Scheme of ENs and EHs participation in the day-ahead energy wholesale and retail markets.Constraints (3)–(23) related to the upper-level problem refer to LOPF equations in electricity, gas, and heat networks [12–14]. Equations (3)–(14) represent power flow equations, in which Eqs. (3)–(7) show power flow of the electricity network and refer to active and reactive power balance at the electrical network buses, (3) and (4), active and reactive power flow on the line (5)–(6), and slack bus voltage (7) [1, 2]. The active and reactive power of the line is nonlinear asPLe,j=GLe,jVe2−VeVjGLe,jcosφe−φj+BLe,jsinφe−φj and QLe,j=−BLe,jVe2+VeVjBLe,jcosφe−φj−GLe,jsinφe−φj, respectively [2]. Yet, because the difference of voltage angles between any two given buses, φe−φj, is smaller than 6° or 0.105 rad in the distribution network [2], the terms cosφe−φj and sinφe−φj are approximately 1 and φe−φj, respectively [1]. Regarding that the voltage magnitude of a bus is within [0.9 1.1] p.u., it will be always near 1 p.u [12–14]. The approximation of voltage magnitude, V, is 1 + ∆V, where ∆V is the voltage deviation [6]. Accordingly, considering the small values of ∆V2 and ΔV (φe - φj) because of the small values of φe−φj and ΔV, nonlinear constraints are substituted by constraints (5) and (6) and the LOPF model is obtained [1]. Equations (8)–(12) deal with power flow equations of the gas network [12]. The balance of gas power at different nodes is modeled in Reference (8). Also, Eq. (9) calculates the gas power passing in the pipeline, where the squared power of the auxiliary variable ϑ, that is, ϑ2, is in the form of signρg,ρjρg2−ρj2 [14]. Consequently, Eqs. (10)–(12) give the linearized model of this equation using the piecewise linearization technique [14]. Equation (10) represents the linear form of the mentioned equation and Eqs. (11) and (12) denote ϑ and ρ in terms of their deviation variables, that is, ∆ϑ and ∆ρ, respectively [13]. Equations (13)–(14) provide power flow of the heat network, sowing the balance of heat power at different nodes and heat power of the pipeline [6].Operation boundaries of energy networks are given in Reference (15)–(23) [12–14] and are related to the deviation of bus voltages, size (apparent power) of distribution lines and substations, the allowable range of gas pressure at various nodes, constraint of the size (gas power) of gas pipelines and stations, the limit on the temperature of nodes, and the size (heat power) of heat pipelines and stations. In the power flow model of the electricity network, Eqs. (3)–(7), a voltage deviation variable is adopted to perform LOPF [2]. As a result, the voltage deviation boundary shown by Reference (15) replaces the voltage magnitude boundary, V¯≤Ve≤V¯ [1]. Eqs. (16) and (17) represent the linearized boundary applied to the apparent power of electrical lines and substations [2]. These are bounded by a circular plane [1]. The center of the circle is in the origin, and its radius S is constrained by P2 + Q2 ≤ S2 and has a nonlinear format [14]. Nonetheless, this limit is approximated by a regular polygon (16) and (17) [14]. If the number of sides of the circular plane is significant, the computational error can be neglected in comparison with that of the circular plane [2].Lower-level problem: This problem is described by Eqs. (24)–(43), where PEH, QEH, HEH, and GEH are found. The lower-level problem deals with the optimal participation of EHs in the retail day-ahead energy market while considering flexibility constraints and the operation model of sources and ALs. As a result, the lower-level problem, that is, Eq. (24), attempts to minimize the expected cost of EHs in retail markets (function of CostEHs) [14]. According to Eq. (24), CostEHs represents the difference between the EHs energy cost and EHs’ revenue in electricity, gas, and heat retail markets. The revenue (energy cost) in each of these markets is equal to the energy price multiplied by power [14]. Note that if EH consumes energy in any of the electricity, gas, and heat markets, it purchases it at a price of λR. Yet, if it can inject energy, it sells it for less than λR, that is, for λW. Equations (25)–(43) give the constraints of the lower-level problem. The balance of active, reactive, heat, and gas powers between supply and demand in the EH is modeled by Eqs. (25)–(28). In these equations, the active, gas, or heat power of EH is divided into power consumption (P−, H−, or G−) and generation power (P+, H+, or G+), each of which has a positive value according to the constraints (25)–(28).Constraints (29)–(32) formulate the operation of the CHP [23]; the heat and gas power of the CHP are found in (29) and (30) in terms of CHP’s active power. Constraints (31)–(32) present the limits on the CHP output in both electricity and heat sections. Constraint (33) gives the input gas power of the boiler using its output heat power, and its output capacity is modeled by Eq. (34) [23]. The operation of ESSs is formulated by Reference (35)–(38) [24], showing the charging and discharging boundaries, the range of energy storable in ESSs, and the boundary of the capacity (apparent power) of the ESS charger, respectively. The operation model of EVs in the parking lot can be represented by Reference (35)–(38) [1, 2, 14], except for parameters CR, DR, EI, and S¯E that change over time [14]. As a result, subscript t is used in the model of EVs parking lot for these parameters. Equation (37) can be represented by two constraints 0≤EIi+∑τ=1tηCHPCHi,τ,s−1/ηDISPDISi,τ,s and EIi+∑τ=1tηCHPCHi,τ,s−1/ηDISPDISi,τ,s=E¯i in modeling the EVs, showing the range of energy storable in EVs and their energy demand in the operation horizon [1]. Parameters CR, DR, and S¯E for EVs parking lot at hour t are obtained by CR=∑ev=1Ntcrev,DR=∑ev=1Ntdrev, and S¯E=∑ev=1NtS¯evEV, wherecr, dr, S¯EV, and Nt represent the charging rate, discharging rate, the size of EV charger, and the number of EVs connected to the EH during hour t [2]. Moreover, EI at hour t equals the sum of total arrival (initial) energy of EVs connected to the parking lot in the EH recently [1]. Equation (39) describes the operation model of RESs, and the boundary of the apparent power of these sources is given [2]. The operation model of responsive loads in electricity, gas, and heat networks can be shown in Reference (40)-(41) [1]. Consumers of this type of DRP act as per Eq. (24) following the price signal; consequently, the EH profit becomes maximum in the energy market [1]. Those consumers that take part in the DRP receive energy from the EH during low energy prices while injecting energy into the EH during high energy prices. Noting that inexpensive (expensive) energy price is in proportion to off-peak (peak) hours, the suggested DRP will reduce energy demanded by consumers over peak intervals and feed this energy at off-peak hours [1]. Equation (40) denotes the limit on active, gas, and heat power control of the DRP. Moreover, Eq. (41) shows that electrical (gas or heat) energy reduced at peak hours can be feed at off-peak hours [1].Constraints (42) and (43) present EH flexibility limits in the electricity and heat sections, respectively. Note that due to the uncertainty in the power generation of RESs, the results of real-time and day-ahead operation will not be the same [18, 19]. Therefore, the balance of production and consumption in real-time operation may not be observed [24]. This is known as low system flexibility [24]. To compensate for this, flexible sources such as responsive loads, energy storage, and nonrenewable sources are used alongside RESs [23]. It is therefore expected that a highly flexible system will be achieved [23]. But estimating this requires a limitation of flexibility. The goal is to minimize the deviation of the injected power of the system to the network in real-time operations compared with day-ahead operation. Therefore, considering the low deviation between the active power of EH in scenario s compared with the scenario corresponding to the deterministic model of the predicted values of uncertainty parameters (here it is assumed that s = 1) is obtained as constraint (42). In addition, according to Eq. (29), the heat power of CHP depends on its active power. Since changes in RESs power in real-time operation compared with day-ahead operation cause changes in ACIT CHP power due to uncertainty in the power generation of RESs, it is expected that the flexibility of the EH’s heat sector be low. Boilers and heat-responsive loads have been used to improve the flexibility of this sector. To estimate the high flexibility in the heat sector of EH, constraint (43) is used in the proposed problem. In Eqs. (42) and (43), the term εF represents the tolerance of flexibility so that if a value of zero is given, 100% flexibility will be available for the EH.


### 2.2. Uncertainty Modeling

Uncertainty parameters of the problem given in Eqs. (1)–(43) include the parameters of load, PD, QD, HD, and GD; energy price, λEW, λHW, λGW, λER, λHR, and λGR; renewable power, PR; charging and discharging rates and charger capacity of EVs, CR, DR, and S¯E; initial energy of EVs, EI; and energy demanded by EVs, E¯. The UT method is employed here as a potent technique for stochastic modeling of uncertainties. The number of runs in this method is very fewer than that of Monte Carlo simulation (MCS) and other analytical approaches; hence, the calculations are straightforward and computational time is less. The method assumes no mathematical assumptions concerning modeling. Besides, the UT method shows suitable behavior for nonlinear correlated/uncorrelated transitions [25, 26] and acceptable probability distribution function (PDF) estimation, and the coding is facilitated. Parameters of the stochastic uncertainties are modeled using the UT method. Parameter n denotes the dimension of the vector of input uncertain parameters (U) and is considered to be 16. The number of scenarios is assumed 2n + 1, which is set to 33 for the problem under study.Thanks to few scenarios generated by this method, scenario reduction methods are not required for facilitating the calculations and reducing the computational time. More information on the proposed method with more focus on the mathematical formulation and implementation algorithm can be extracted from Reference [25].The uncertain nonlinear stochastic problem is expressed asy = f(x). In this formula, y ∈ Rr is the uncertain output vector that has r elements, and x ∈ Rn is the uncertain input vector with mean and covariance values of μx and σx. The symmetrical and nonsymmetrical elements of σx give the variance of uncertain variables and covariance of several uncertain parameters, respectively. Mean and covariance output variables, that is, μy and σy, are obtained using the UT method [25]. The procedure of this method is elaborated as follows:Step 1. Take 2n+1 samples from the input uncertain data:(46)x0=μx,(47)xϖ=μx+n1−W0σx,∀ϖ=1,2,...,n,(48)xϖ=μx−n1−W0σx,∀ϖ=1,2,...,n.In the above equations,W0 is the weight of μx.Step 2. Check the weighting factor of samples:(49)W0=W0,(50)Wϖ=1−W02n,∀ϖ=1,2,...,n,(51)Wϖ+n=1−W02n∀ϖ+n=n+1,n+2,...,2n,(52)∑ϖ=1nWϖ=1.Step 3. Take 2n+1 points of the nonlinear function and find output samples:(53)yϖ=fxϖ.Step 4. Check σy and μy of the output variable ν:(54)μy=∑ϖ=1nWϖνϖ,(55)σy=∑ϖ=1nWϖνϖ−μy−νϖ−μyT.Figure2 illustrates the flowchart of the UT method.Figure 2
Flowchart of UT method.


## 2.1. Mathematical Formulation

The two-layer energy management of EHs in electricity, gas, and heat networks is described here. The first layer deals with establishing coordination among EHOs, sources, and ALs. The second layer coordinates EHOs and ENOs. As a result, we have an optimization problem with two upper and lower levels. The former is related to the participation of ENs in retail and wholesale markets, while the latter manages the power supply and demand of the second layer of the energy management plan. The suggested plan is mathematically formulated:(1)minCostENs=∑s∈ΠSπs∑t∈ΠOHλEWt,sPESr,t,s+λHWt,sHHSr,t,s+λGWt,sGGSr,t,s+∑s∈ΠSπs∑t∈ΠOH∑i∈ΠEHλEWt,sPEHi,t,s++λHWt,sHEHi,t,s++λGWt,sGEHi,t,s+,(2)minEEL=∑s∈ΠSπs∑t∈ΠOHPESr,t,s+∑i∈ΠEHPEHi,t,s+−PEHi,t,s−−∑e∈ΠENPDe,t,s+∑s∈ΠSπs∑t∈ΠOHHHSr,t,s+∑i∈ΠEHHEHi,t,s+−HEHi,t,s−−∑h∈ΠHNHDh,t,s+∑s∈ΠSπs∑t∈ΠOHGGSr,t,s+∑i∈ΠEHGEHi,t,s+−GEHi,t,s−−∑g∈ΠGNGDg,t,s,which is subjected to(3)PESe,t,s−PDe,t,s+∑i∈ΠEHAEe,iPEHi,t,s+−PEHi,t,s−=∑j∈ΠENBEe,jPLe,j,t,s,∀e,t,s,(4)QESe,t,s−QDe,t,s+∑i∈ΠEHAEe,iQEHi,t,s=∑j∈ΠENBEe,jQLe,j,t,s,∀e,t,s,(5)PLe,j,t,s=GLe,jΔVe,t,s−ΔVj,t,s−BLe,jφe,t,s−φj,t,s,∀e,j,t,s,(6)QLe,j,t,s=−BLe,jΔVe,t,s−ΔVj,t,s−GLe,jφe,t,s−φj,t,s,∀e,j,t,s,(7)ΔVe,t,s,φe,t,s=0,∀e=r,t,s,(8)GGSg,t,s−GDg,t,s+∑i∈ΠEHAGg,iGEHi,t,s+−GEHi,t,s−=∑j∈ΠGNBGg,jGLg,j,t,s,∀g,t,s,(9)GLg,j,t,s=κg,jsignρg,t,s,ρj,t,sϑg,j,t,s,(10)Δϑg,j,t,s,p=signρg,t,s,ρj,t,sspρspϑΔρg,t,s,p−Δρj,t,s,p,∀g,j,t,s,l,(11)ϑg,j,t,s=∑p∈ΠPΔϑg,j,t,s,p,∀g,j,t,s(12)ρg,t,s=ρ¯g+∑p∈ΠPΔρg,t,s,p,∀g,t,s,(13)HESh,t,s−HDh,t,s+∑i∈ΠEHAHh,iHEHi,t,s+−HEHi,t,s−=∑j∈ΠENBHe,jHLh,j,t,s,∀h,t,s,(14)ρg,t,s=ρ¯g+∑p∈ΠPΔρg,t,s,p,∀g,t,s,(15)V¯e−1≤ΔVe,t,s≤V¯e−1,∀e,t,s,(16)PLe,j,t,scosk⋅Δβ+QLe,j,t,ssink⋅Δβ≤S¯Le,j∀e,j,t,s,k∈ΠK=1,2,…,nk,Δβ=360nk,(17)PESe,t,scosk⋅Δβ+QESe,t,ssink⋅Δβ≤S¯ESe,∀e=r,t,s,k,(18)ρ¯g≤ρg,t,s≤ρ¯g,∀g,t,s,(19)−G¯Lg,j≤GLg,j,t,s≤G¯Lg,j,∀g,j,t,s,(20)−G¯GSg≤GGSg,t,s≤G¯GSgm,∀g=r,t,s,(21)T¯h≤Th,t≤Th¯,∀h,t,(22)−H¯Lh,j≤HLh,j,t,s≤H¯Lh,j,∀h,j,t,s,(23)−H¯HSh≤HHSh,t,s≤H¯HSh,∀h=r,t,s(24)−H¯Lh,j≤HLh,j,t,s≤H¯Lh,j,∀h,j,t,s,(25)−H¯HSh≤HHSh,t,s≤H¯HSh,∀h=r,t,s,(26)PEH+,PEH−,QEH,HEH+,HEH−,GEH+,GEH−∈argminCostEHs=∑s∈ΠSπs∑t∈ΠOH∑i∈ΠEHλERt,sPEHi,t,s−+λHRt,sHEHi,t,s−+λGRt,sGEHi,t,s−−λEWt,sPEHi,t,s++λHWt,sHEHi,t,s++λGWt,sGEHi,t,s+,which is subjected to(27)PEHi,t,s+−PEHi,t,s−=PCi,t,s+PRi,t,s+PDRi,t,s+PDIS,t,s−PCHi,t,s−PDi,t,s,∀PEHi,t,s+,PEHi,t,s−≥0&i,t,s,(28)QEHi,t,s=QCi,t,s+QRi,t,s+QEi,t,s−QDi,t,s,∀i,t,s,(29)HEHi,t,s+−HEHi,t,s−=HCi,t,s+HBi,t,s+HDRi,t,s−HDi,t,s,∀HEHi,t,s+,HEHi,t,s−≥0&i,t,s,(30)GEHi,t,s+−GEHi,t,s−=GDRi,t,s−GCi,t,s−GBi,t,s−GDi,t,s,∀GEHi,t,s+,GEHi,t,s−≥0,&i,t,s,(31)HCi,t,s=PCi,t,s1−ηT−ηLηT⋅ηH,∀i,t,s,(32)GCi,t,s=1ηT.PCi,t,s,∀i,t,s,(33)PCi,t,scosk⋅Δβ+QCi,t,scosk⋅Δβ≤S¯Ci,∀i,t,s,k,(34)0≤HCi,t,s≤H¯Ci,∀i,t,s,(35)GBi,t,s=1ηB⋅HBi,t,s,∀i,t,s,(36)0≤HBi,t,s≤H¯Bi,∀i,t,s,(37)0≤PCHi,t,s≤CRi,∀i,t,s,(38)0≤PDISi,t,s≤DRi,∀i,t,s,(39)E¯i≤EIi+∑τ=1tηCHPCHi,τ,s−1ηDISPDISi,τ,s≤E¯i,∀i,t,s,(40)PDISi,t,s−PCHi,t,scosk⋅Δβ+QEi,t,scosk⋅Δβ≤S¯Ei,∀i,t,s,k,(41)PRi,t,scosk⋅Δβ+QRi,t,scosk⋅Δβ≤S¯Ri,∀i,t,s,k,(42)−υiPDi,t,sHDi,t,sorGDi,t,s≤PDRi,t,sHDRi,t,sorGDRi,t,s≤υiPDi,t,sHDi,t,sorGDi,t,s,∀i,t,s,(43)∑t∈ΠOHPDRi,t,sHDRi,t,sorGDRi,t,s=0,∀i,s,(44)−εF≤PEHi,t,s−PEHi,t,s=1≤εF,∀i,t,s,(45)−εF≤HEHi,t,s−HEHi,t,s=1≤εF,∀i,t,s.Upper-level problem: Eqs. (1)–(23) show the upper-level problem as a two-objective optimization. Equation (1) minimizes the expected energy cost of ENs (as a function of CostENs), where CostENs represents the sum of power purchased by ENs from upstream networks (the first term in Eq. (1)) [12] and EHs (the second term in Eq. (1)) [14]. Two energy markets, that is, wholesale and retail energy markets, are considered in the proposed plan, Figure 1. In the former, the wholesale market supplies energy to ENs. The first term in Eq. (1) models the wholesale market, where the positive value of PDS (HHS or GGS) means that energy is purchased from the market, while its negative value means that the wholesale market purchases energy from the network [14]. Additionally, EHs (passive customers) purchase or sell (purchase) energy through participating in a retail market. It is noteworthy that EHs are assumed to buy energy from ENs at the retail market price, λR, which is generally higher than the wholesale market price, λW. Nonetheless, if EHs can inject energy into ENs, ENs purchase this energy at a price of λW, which is generally less than λR. This is due to maintaining the favorable financial benefit of ENs in energy markets and the participation of EHs in energy production. This is due to the fact that a wholesaler generally buys his goods at a price lower than his sale price to a retailer. Note that based on Figure 1, two models of the wholesale and retail market are considered in this paper. It is assumed that the distribution company (energy network) purchases energy from the wholesale market. It then sells the energy it buys within itself as a retail market to its customers, that is, consumers and energy hubs. To encourage its customers in network energy management, it also allows its customers to sell energy in the retail market. Therefore, it can be stated that a distribution company is a wholesale buyer who provides its energy to customers in the form of retail sales. Hence, in the second part of Eq. (1), the price λW is used to buy energy by ENs from EHs in the retail market. Equation (2) deals with the operation of ENOs and minimizes the total expected energy loss in all energy networks (function of EEL) [23].Figure 1
Scheme of ENs and EHs participation in the day-ahead energy wholesale and retail markets.Constraints (3)–(23) related to the upper-level problem refer to LOPF equations in electricity, gas, and heat networks [12–14]. Equations (3)–(14) represent power flow equations, in which Eqs. (3)–(7) show power flow of the electricity network and refer to active and reactive power balance at the electrical network buses, (3) and (4), active and reactive power flow on the line (5)–(6), and slack bus voltage (7) [1, 2]. The active and reactive power of the line is nonlinear asPLe,j=GLe,jVe2−VeVjGLe,jcosφe−φj+BLe,jsinφe−φj and QLe,j=−BLe,jVe2+VeVjBLe,jcosφe−φj−GLe,jsinφe−φj, respectively [2]. Yet, because the difference of voltage angles between any two given buses, φe−φj, is smaller than 6° or 0.105 rad in the distribution network [2], the terms cosφe−φj and sinφe−φj are approximately 1 and φe−φj, respectively [1]. Regarding that the voltage magnitude of a bus is within [0.9 1.1] p.u., it will be always near 1 p.u [12–14]. The approximation of voltage magnitude, V, is 1 + ∆V, where ∆V is the voltage deviation [6]. Accordingly, considering the small values of ∆V2 and ΔV (φe - φj) because of the small values of φe−φj and ΔV, nonlinear constraints are substituted by constraints (5) and (6) and the LOPF model is obtained [1]. Equations (8)–(12) deal with power flow equations of the gas network [12]. The balance of gas power at different nodes is modeled in Reference (8). Also, Eq. (9) calculates the gas power passing in the pipeline, where the squared power of the auxiliary variable ϑ, that is, ϑ2, is in the form of signρg,ρjρg2−ρj2 [14]. Consequently, Eqs. (10)–(12) give the linearized model of this equation using the piecewise linearization technique [14]. Equation (10) represents the linear form of the mentioned equation and Eqs. (11) and (12) denote ϑ and ρ in terms of their deviation variables, that is, ∆ϑ and ∆ρ, respectively [13]. Equations (13)–(14) provide power flow of the heat network, sowing the balance of heat power at different nodes and heat power of the pipeline [6].Operation boundaries of energy networks are given in Reference (15)–(23) [12–14] and are related to the deviation of bus voltages, size (apparent power) of distribution lines and substations, the allowable range of gas pressure at various nodes, constraint of the size (gas power) of gas pipelines and stations, the limit on the temperature of nodes, and the size (heat power) of heat pipelines and stations. In the power flow model of the electricity network, Eqs. (3)–(7), a voltage deviation variable is adopted to perform LOPF [2]. As a result, the voltage deviation boundary shown by Reference (15) replaces the voltage magnitude boundary, V¯≤Ve≤V¯ [1]. Eqs. (16) and (17) represent the linearized boundary applied to the apparent power of electrical lines and substations [2]. These are bounded by a circular plane [1]. The center of the circle is in the origin, and its radius S is constrained by P2 + Q2 ≤ S2 and has a nonlinear format [14]. Nonetheless, this limit is approximated by a regular polygon (16) and (17) [14]. If the number of sides of the circular plane is significant, the computational error can be neglected in comparison with that of the circular plane [2].Lower-level problem: This problem is described by Eqs. (24)–(43), where PEH, QEH, HEH, and GEH are found. The lower-level problem deals with the optimal participation of EHs in the retail day-ahead energy market while considering flexibility constraints and the operation model of sources and ALs. As a result, the lower-level problem, that is, Eq. (24), attempts to minimize the expected cost of EHs in retail markets (function of CostEHs) [14]. According to Eq. (24), CostEHs represents the difference between the EHs energy cost and EHs’ revenue in electricity, gas, and heat retail markets. The revenue (energy cost) in each of these markets is equal to the energy price multiplied by power [14]. Note that if EH consumes energy in any of the electricity, gas, and heat markets, it purchases it at a price of λR. Yet, if it can inject energy, it sells it for less than λR, that is, for λW. Equations (25)–(43) give the constraints of the lower-level problem. The balance of active, reactive, heat, and gas powers between supply and demand in the EH is modeled by Eqs. (25)–(28). In these equations, the active, gas, or heat power of EH is divided into power consumption (P−, H−, or G−) and generation power (P+, H+, or G+), each of which has a positive value according to the constraints (25)–(28).Constraints (29)–(32) formulate the operation of the CHP [23]; the heat and gas power of the CHP are found in (29) and (30) in terms of CHP’s active power. Constraints (31)–(32) present the limits on the CHP output in both electricity and heat sections. Constraint (33) gives the input gas power of the boiler using its output heat power, and its output capacity is modeled by Eq. (34) [23]. The operation of ESSs is formulated by Reference (35)–(38) [24], showing the charging and discharging boundaries, the range of energy storable in ESSs, and the boundary of the capacity (apparent power) of the ESS charger, respectively. The operation model of EVs in the parking lot can be represented by Reference (35)–(38) [1, 2, 14], except for parameters CR, DR, EI, and S¯E that change over time [14]. As a result, subscript t is used in the model of EVs parking lot for these parameters. Equation (37) can be represented by two constraints 0≤EIi+∑τ=1tηCHPCHi,τ,s−1/ηDISPDISi,τ,s and EIi+∑τ=1tηCHPCHi,τ,s−1/ηDISPDISi,τ,s=E¯i in modeling the EVs, showing the range of energy storable in EVs and their energy demand in the operation horizon [1]. Parameters CR, DR, and S¯E for EVs parking lot at hour t are obtained by CR=∑ev=1Ntcrev,DR=∑ev=1Ntdrev, and S¯E=∑ev=1NtS¯evEV, wherecr, dr, S¯EV, and Nt represent the charging rate, discharging rate, the size of EV charger, and the number of EVs connected to the EH during hour t [2]. Moreover, EI at hour t equals the sum of total arrival (initial) energy of EVs connected to the parking lot in the EH recently [1]. Equation (39) describes the operation model of RESs, and the boundary of the apparent power of these sources is given [2]. The operation model of responsive loads in electricity, gas, and heat networks can be shown in Reference (40)-(41) [1]. Consumers of this type of DRP act as per Eq. (24) following the price signal; consequently, the EH profit becomes maximum in the energy market [1]. Those consumers that take part in the DRP receive energy from the EH during low energy prices while injecting energy into the EH during high energy prices. Noting that inexpensive (expensive) energy price is in proportion to off-peak (peak) hours, the suggested DRP will reduce energy demanded by consumers over peak intervals and feed this energy at off-peak hours [1]. Equation (40) denotes the limit on active, gas, and heat power control of the DRP. Moreover, Eq. (41) shows that electrical (gas or heat) energy reduced at peak hours can be feed at off-peak hours [1].Constraints (42) and (43) present EH flexibility limits in the electricity and heat sections, respectively. Note that due to the uncertainty in the power generation of RESs, the results of real-time and day-ahead operation will not be the same [18, 19]. Therefore, the balance of production and consumption in real-time operation may not be observed [24]. This is known as low system flexibility [24]. To compensate for this, flexible sources such as responsive loads, energy storage, and nonrenewable sources are used alongside RESs [23]. It is therefore expected that a highly flexible system will be achieved [23]. But estimating this requires a limitation of flexibility. The goal is to minimize the deviation of the injected power of the system to the network in real-time operations compared with day-ahead operation. Therefore, considering the low deviation between the active power of EH in scenario s compared with the scenario corresponding to the deterministic model of the predicted values of uncertainty parameters (here it is assumed that s = 1) is obtained as constraint (42). In addition, according to Eq. (29), the heat power of CHP depends on its active power. Since changes in RESs power in real-time operation compared with day-ahead operation cause changes in ACIT CHP power due to uncertainty in the power generation of RESs, it is expected that the flexibility of the EH’s heat sector be low. Boilers and heat-responsive loads have been used to improve the flexibility of this sector. To estimate the high flexibility in the heat sector of EH, constraint (43) is used in the proposed problem. In Eqs. (42) and (43), the term εF represents the tolerance of flexibility so that if a value of zero is given, 100% flexibility will be available for the EH.


## 2.2. Uncertainty Modeling

Uncertainty parameters of the problem given in Eqs. (1)–(43) include the parameters of load, PD, QD, HD, and GD; energy price, λEW, λHW, λGW, λER, λHR, and λGR; renewable power, PR; charging and discharging rates and charger capacity of EVs, CR, DR, and S¯E; initial energy of EVs, EI; and energy demanded by EVs, E¯. The UT method is employed here as a potent technique for stochastic modeling of uncertainties. The number of runs in this method is very fewer than that of Monte Carlo simulation (MCS) and other analytical approaches; hence, the calculations are straightforward and computational time is less. The method assumes no mathematical assumptions concerning modeling. Besides, the UT method shows suitable behavior for nonlinear correlated/uncorrelated transitions [25, 26] and acceptable probability distribution function (PDF) estimation, and the coding is facilitated. Parameters of the stochastic uncertainties are modeled using the UT method. Parameter n denotes the dimension of the vector of input uncertain parameters (U) and is considered to be 16. The number of scenarios is assumed 2n + 1, which is set to 33 for the problem under study.Thanks to few scenarios generated by this method, scenario reduction methods are not required for facilitating the calculations and reducing the computational time. More information on the proposed method with more focus on the mathematical formulation and implementation algorithm can be extracted from Reference [25].The uncertain nonlinear stochastic problem is expressed asy = f(x). In this formula, y ∈ Rr is the uncertain output vector that has r elements, and x ∈ Rn is the uncertain input vector with mean and covariance values of μx and σx. The symmetrical and nonsymmetrical elements of σx give the variance of uncertain variables and covariance of several uncertain parameters, respectively. Mean and covariance output variables, that is, μy and σy, are obtained using the UT method [25]. The procedure of this method is elaborated as follows:Step 1. Take 2n+1 samples from the input uncertain data:(46)x0=μx,(47)xϖ=μx+n1−W0σx,∀ϖ=1,2,...,n,(48)xϖ=μx−n1−W0σx,∀ϖ=1,2,...,n.In the above equations,W0 is the weight of μx.Step 2. Check the weighting factor of samples:(49)W0=W0,(50)Wϖ=1−W02n,∀ϖ=1,2,...,n,(51)Wϖ+n=1−W02n∀ϖ+n=n+1,n+2,...,2n,(52)∑ϖ=1nWϖ=1.Step 3. Take 2n+1 points of the nonlinear function and find output samples:(53)yϖ=fxϖ.Step 4. Check σy and μy of the output variable ν:(54)μy=∑ϖ=1nWϖνϖ,(55)σy=∑ϖ=1nWϖνϖ−μy−νϖ−μyT.Figure2 illustrates the flowchart of the UT method.Figure 2
Flowchart of UT method.


## 3. Integrated Model of the Proposed Formulation

### 3.1. Single-Objective Scheme

To solve the upper-level problem, the Pareto optimization is applied so that the set of desired solutions is found and presented to decision makers to select the final solution.ε-Constraint-based Pareto optimization is adopted to develop a model for the problem [27]. Function CostENs given by Eq. (1) is utilized as the objective function in the new problem, and function EEL (2) is constrained to εEEL:(56)minCostENs=∑s∈ΠSπw∑t∈ΠOHλEWt,sPESr,t,s+λHWt,sHHSr,t,s+λGWt,sGGSr,t,s+∑s∈ΠSπw∑t∈ΠOH∑i∈ΠEHλEWt,sPEHi,t,s++λHWt,sHEHi,t,s++λGWt,sGEHi,t,s+,which is subjected to(58)∑s∈ΠSπw∑t∈ΠOHPESr,t,s+∑i∈ΠEHPEHi,t,s+−PEHi,t,s−−∑e∈ΠENPDe,t,s+∑s∈ΠSπw∑t∈ΠOHHHSr,t,s+∑i∈ΠEHHEHi,t,s+−HEHi,t,s−−∑e∈ΠENHDe,t,s+∑s∈ΠSπw∑t∈ΠOHGGSr,t,s+∑i∈ΠEHGEHi,t,s+−GEHi,t,s−−∑e∈ΠENGDe,t,s≤εEEL.In Equations (3)–(43), the value of εEEL is in the range of minimum and maximum values of EEL, that is, EELmin and EELmax. So, εEEL is corresponding to a specific value of CostENs and EEL. The set of points achieved this way is named the Pareto front [27]. The fuzzy decision approach [27] is incorporated to find the best compromise solution. Details of the suggested method are provided here [27]:Step 1. Find the values of the linear fuzzy membership function (f^) of individual objective functions for different εEEL:(59)f^i=1fi≤fiminfi−fimaxfimin−fimax0fimin≤f≤fimaxf≥fimax,∀i=CostENs,EEL.Step 2. Find the value of ar=minf^CostENsr,f^EELr, r = 1, 2, …. This shows different steps with different values of εEEL.Step 3. Relate the best compromise solution to a point containing max(a1, a2, …).


### 3.2. Single-Level Problem

Problems (54)–(56) are structured as a bilevel formulation. Traditional methods and solvers need to extract a single-level model to reach the optimal solution [28]. The present paper adopts the KKT.The model of the problem, given in subsection2.1, has a general form (58)–(62). Equations (58) and (59) represent the upper-level problems (1)–(23), while Eqs. (60)–(62) formulate the lower-level problems (24)–(43). x (y) is the vector of variables of the upper-level (lower-level) problem. Parameters ρ and μ represent the Lagrange multipliers.(60)minF1=aTx+bTy,which is subjected to(61)c1.x+d1y≤/=/≥e1,(62)y∈argminF2=fTy,which is subjected to(63)g1y=h1:ρ,(64)g2y≤h2:μ.Constraints found by KKT of the lower-level problem are added to the upper-level problem to reach the single-objective model for the problem [28]. The Lagrange function (L) of the lower-level problem should be determined (63) to find these constraints. The Lagrange function of a problem is the sum of its objective function and penalty functions. The penalty function for constraints a ≤ b and a = b is represented by μ.max(0, a - b) and ρ.(b - a), respectively [28].(65)L=F2+ρ⋅h1−g1y+μ⋅max0,g2y−h2.Constraints obtained from KKT are related to taking the derivative of the Lagrange function with respect to its variables (y,μ, and ρ) equal to zero [28]. As a result, a single-level formulation for problem (58)–(62) will be given by Eqs. (64)–(69). In the newly formed problem, Eqs. (64) and (65), which are the same as Eqs. (58) and (59), model the upper-level problem. Constraint (66) can be found by taking the derivative of the Lagrange function with respect to the primal variable of the lower-level problem (y) equal to zero. Constraint (67) is formed from ∂L/∂ρ=0, which will be the same as constraint (61). The result of ∂L/∂μ=0 (where μ represents the Lagrange multiplier of an inequality constraint) is bounded by two conditions (68), where constraint (62) is resulted based on its first condition, and constraint μ⋅g2y−h2=0 is found based on the second condition. This constraint forms a nonlinear equation that needs to be linearized. To this end, constraints −M⋅z≤μ≤M⋅zand −M⋅1−z≤g2y−h2≤M⋅1−zare substituted for μ⋅g2y−h2=0. Parameter M denotes a large fixed number, for instance, 106, and z is a binary variable [28]. Equation (69) gives the boundary on Lagrange multipliers. Finally, Figure 3 shows the flowchart of solution process.(66)minF1=aTx+bTy,which is subjected to(67)g1y=h1:ρ,(68)g2y≤h2:μ,(69)Constraint(70)∂L∂y=0⇒g1ρ+g2μ=f,(71)∂L∂ρ=0⇒Constraint61,(72)∂L∂μ=0⇒Constraint62∀First conditionμ⋅g2y−h2=0∀Second condition,(73)ρ∈−∞,+∞,μ∈0,+∞,Figure 3
Flowchart of the solution process.


## 3.1. Single-Objective Scheme

To solve the upper-level problem, the Pareto optimization is applied so that the set of desired solutions is found and presented to decision makers to select the final solution.ε-Constraint-based Pareto optimization is adopted to develop a model for the problem [27]. Function CostENs given by Eq. (1) is utilized as the objective function in the new problem, and function EEL (2) is constrained to εEEL:(56)minCostENs=∑s∈ΠSπw∑t∈ΠOHλEWt,sPESr,t,s+λHWt,sHHSr,t,s+λGWt,sGGSr,t,s+∑s∈ΠSπw∑t∈ΠOH∑i∈ΠEHλEWt,sPEHi,t,s++λHWt,sHEHi,t,s++λGWt,sGEHi,t,s+,which is subjected to(58)∑s∈ΠSπw∑t∈ΠOHPESr,t,s+∑i∈ΠEHPEHi,t,s+−PEHi,t,s−−∑e∈ΠENPDe,t,s+∑s∈ΠSπw∑t∈ΠOHHHSr,t,s+∑i∈ΠEHHEHi,t,s+−HEHi,t,s−−∑e∈ΠENHDe,t,s+∑s∈ΠSπw∑t∈ΠOHGGSr,t,s+∑i∈ΠEHGEHi,t,s+−GEHi,t,s−−∑e∈ΠENGDe,t,s≤εEEL.In Equations (3)–(43), the value of εEEL is in the range of minimum and maximum values of EEL, that is, EELmin and EELmax. So, εEEL is corresponding to a specific value of CostENs and EEL. The set of points achieved this way is named the Pareto front [27]. The fuzzy decision approach [27] is incorporated to find the best compromise solution. Details of the suggested method are provided here [27]:Step 1. Find the values of the linear fuzzy membership function (f^) of individual objective functions for different εEEL:(59)f^i=1fi≤fiminfi−fimaxfimin−fimax0fimin≤f≤fimaxf≥fimax,∀i=CostENs,EEL.Step 2. Find the value of ar=minf^CostENsr,f^EELr, r = 1, 2, …. This shows different steps with different values of εEEL.Step 3. Relate the best compromise solution to a point containing max(a1, a2, …).


## 3.2. Single-Level Problem

Problems (54)–(56) are structured as a bilevel formulation. Traditional methods and solvers need to extract a single-level model to reach the optimal solution [28]. The present paper adopts the KKT.The model of the problem, given in subsection2.1, has a general form (58)–(62). Equations (58) and (59) represent the upper-level problems (1)–(23), while Eqs. (60)–(62) formulate the lower-level problems (24)–(43). x (y) is the vector of variables of the upper-level (lower-level) problem. Parameters ρ and μ represent the Lagrange multipliers.(60)minF1=aTx+bTy,which is subjected to(61)c1.x+d1y≤/=/≥e1,(62)y∈argminF2=fTy,which is subjected to(63)g1y=h1:ρ,(64)g2y≤h2:μ.Constraints found by KKT of the lower-level problem are added to the upper-level problem to reach the single-objective model for the problem [28]. The Lagrange function (L) of the lower-level problem should be determined (63) to find these constraints. The Lagrange function of a problem is the sum of its objective function and penalty functions. The penalty function for constraints a ≤ b and a = b is represented by μ.max(0, a - b) and ρ.(b - a), respectively [28].(65)L=F2+ρ⋅h1−g1y+μ⋅max0,g2y−h2.Constraints obtained from KKT are related to taking the derivative of the Lagrange function with respect to its variables (y,μ, and ρ) equal to zero [28]. As a result, a single-level formulation for problem (58)–(62) will be given by Eqs. (64)–(69). In the newly formed problem, Eqs. (64) and (65), which are the same as Eqs. (58) and (59), model the upper-level problem. Constraint (66) can be found by taking the derivative of the Lagrange function with respect to the primal variable of the lower-level problem (y) equal to zero. Constraint (67) is formed from ∂L/∂ρ=0, which will be the same as constraint (61). The result of ∂L/∂μ=0 (where μ represents the Lagrange multiplier of an inequality constraint) is bounded by two conditions (68), where constraint (62) is resulted based on its first condition, and constraint μ⋅g2y−h2=0 is found based on the second condition. This constraint forms a nonlinear equation that needs to be linearized. To this end, constraints −M⋅z≤μ≤M⋅zand −M⋅1−z≤g2y−h2≤M⋅1−zare substituted for μ⋅g2y−h2=0. Parameter M denotes a large fixed number, for instance, 106, and z is a binary variable [28]. Equation (69) gives the boundary on Lagrange multipliers. Finally, Figure 3 shows the flowchart of solution process.(66)minF1=aTx+bTy,which is subjected to(67)g1y=h1:ρ,(68)g2y≤h2:μ,(69)Constraint(70)∂L∂y=0⇒g1ρ+g2μ=f,(71)∂L∂ρ=0⇒Constraint61,(72)∂L∂μ=0⇒Constraint62∀First conditionμ⋅g2y−h2=0∀Second condition,(73)ρ∈−∞,+∞,μ∈0,+∞,Figure 3
Flowchart of the solution process.


## 4. Numerical Results

### 4.1. Problem Data

Figure4 demonstrates the application of the proposed scheme to a test system consisting of a 9-bus electrical network, a 4-node gas network, and a 7-node heating network [14]. The base power in the electrical network is 1 MVA, and it is 1 MW in the gas and heating networks [14]. Moreover, the base voltage, pressure, and temperature are 1 kV, 10 Bar, and 100°C, and their allowable limit is [0.9, 1.1] p.u. [13, 27, 28]. The specifications of distribution lines and gas and heat pipelines are presented in Reference [14]. Peak electrical and thermal data are reported in Reference [14]. In the gas network, CHP and boiler are the only consumers, and the passive gas load is zero. Hourly load data are found by multiplying peak load and load factor [12, 29, 30]. The expected daily load factor curve for electrical and heating networks is shown in Figure 5(a) [14]. The expected daily energy price curve in the wholesale market for various ENs is shown in Figure 5(b) [14]. ENs are assumed to increase the energy price in the retail market by 20% compared with the wholesale market to make a profit, that is, λR = 1.2×λW.Figure 4
Test system [14].Figure 5
Expected daily curves of (a) load factor [14], (b) energy price in Wholesale market [14], (c) RES power rate [31], and (d) EVs penetration rate [32].
(a)(b)(c)(d)The network includes 7 EHs whose locations in different ENs are specified in Figure4, and their load is presented in Reference [14]. EHs 1–3 and 5 have electrical sources and ALs, that is, RESs, EVs, batteries, and electrical responsive loads. The EH 4 also consists of a CHP, a boiler, and heat-responsive loads. Hubs 6 and 7 also contain all of these elements. CHP has a maximum electrical and heat output of 1 MVA and 1 MW, respectively, with ηT, ηL, and ηH being 40%, 9%, and 40%, respectively [23]. The boiler with 80% efficiency also has a maximum heat output of 0.3 MW [14]. RESs including a 0.25 MVA photovoltaic and a 0.2 MVA wind system are used in EHs [14]. The hourly power of RES is calculated by multiplying its capacity and the rate of power output [13]. The expected daily curve of power output for photovoltaics and wind systems is illustrated in Figure 5(c) [31]. The battery with a charging and discharging efficiency of 88% and capacity of 2 MWh is installed in hubs 1–3 and 5–7 [14]. The battery has a charging and discharging rate of 0.5 MW and charger capacity of 0.6 MVA, its minimum storable energy, and initial energy is 0.2 MWh [14]. In each of the EHs 1–3 and 5–7, 60 EVs can be connected, the specifications of each EV including charging/discharging rate, charger capacity, and other items are presented in Reference [1,2]. The number of EVs connected to the EH per hour is proportional to the product of the total number of EVs in the parking lot and the penetration rate of EV into the parking lot [32]. Figure 5(d) plots the expected daily curve of EVs penetration rate [32]. Hubs are assumed to take part in the DRP scheme with a penetration rate of 30% [14]. To achieve the desired flexibility in EH, a flexibility tolerance of 0.02 MW is considered.


### 4.2. Results

In this subsection, the proposed scheme applied to the data given in subsection4.1 is simulated using GAMS, and then, the CPLEX solves the problem [33]. The following are the details of the numerical results obtained from various case studies. It is noteworthy that here 5 linear pieces were used in the conventional piecewise linearization method, and the circular plane is approximated to a regular 45-gon. According to these conditions and regarding the results obtained in [14], the computational error for active, reactive, gas, and heat power in the linear approximation model used compared with the original nonlinear model is about 2%, 2%, 0.9%, and 0, respectively. Also, these values for voltage, pressure, and temperature are roughly 0.5%, 0.1%, and 0, respectively.A)Performance evaluation of EHs: The expected daily curve of active, reactive, heat, and gas power of EHs and their sources and ALs is plotted in Figure 6. Based on Figure 6(a) and according to Figure 5(c) and the data of subsection 4.1, RESs such as photovoltaics and wind systems inject active power equal to the maximum active power into EHs in all operation hours according to the weather conditions. As their operating costs are negligible (zero in this paper), EHs tend to make high use of the active power generation of RESs in order to minimize their energy costs in the retail market, as in Eq. (24). In the case of CHP, it is observed that they feed higher amounts of active power into EHs during the hours when the price of electric energy is higher than the price of gas energy according to Figure 5(b) so that at 1:00–4:00, 8:00–16:00, and 23:00-00:00, feed active power is equal to their maximum capacity, that is, 3 MW (capacity of each CHP is 1 MVA and three EH contain CHPs) into EHs. Also, according to Figure 5(a), there is peak electric load during 17:00-22:00, that is, active and reactive power consumed during these hours is high. Hence, the voltage drop may be high at this period. Therefore, in addition to injecting active power, CHPs also tend to inject reactive power during these hours, so the active power produced by CHPs is less than their maximum capacity. Also, according to Figure 5(b), at 5:00-7:00, the price of gas energy is higher than the price of electricity. Therefore, CHPs produce low active power during these hours to minimize the energy cost of EHs, but they are switched off. The reason why these sources are not turned off during the mentioned hours can be due to the following:CHPs produce thermal power at their output as well, the amount of which depends on the active power of CHPs. As shown in Figure5(b), the price of heat energy is higher than the price of gas energy at 5:00-7:00. Therefore, considering Eq. (24), CHP tends to produce heat power, which requires active power generation based on Eq. (29).CHPs are a flexible source for RESs. As shown in Figure5(c), RESs generate active power at all hours of operation. Therefore, it is expected that to improve the flexibility of EHs, the flexibility sources must control their active power at all hours of operation, so they will always be on at all hours.Figure 6
Expected daily curves of (a) active, (b) reactive, (c) heating, and (d) gas power of EHs and their sources and ALs with consideringεF = 0.02 MW.
(a)(b)(c)(d)ALs, such as batteries, EVs, and DRPs, perform charging operations when electricity price is low, that is, between 1:00-16:00 and 23:00-00:00, as shown in Figure6(a). During the peak load period when the price of electrical energy has the highest value as per Figure 5(b), ALs perform the discharging operation and inject active power into EHs. Moreover, they are switched on at all hours of operation, except for EVs, which do not inject power into EHs between 10:00 and 11:00 because, according to Figure 5(d), no EVs are connected to the EHs during these hours. This is the function model of ALs to enhance the flexibility of EHs with RESs. Finally, based on Eq. (25), the expected daily curve of the active power of EHs is like the blue curve in Figure 6(a). They appear as consumers of electrical energy in ENs because of high power consumption by ALs and lower power generation by sources during 1:00-8:00. But at other times, they are a generator of electrical energy.In the field of the expected daily curve of reactive power of EHs, it can be seen according to Figure6(b) that RESs and batteries have a time graph similar to the daily curve of reactive charge. Only at 12:00 do photovoltaics inject high active power into EHs as Figure 6(a), so their injected active power into EHs is reduced to zero. Such conditions exist for wind systems at 17:00. CHPs also inject reactive power into EHs only during peak electricity consumption interval, 17:00-22:00, so that it contributes to supplying high reactive power consumption of EHs. EVs also inject high reactive power into EHs between 1:00-7:00 and 18:00-00:00 because during these periods; according to Figure 6(a), EVs receive high active power from the electrical network, which may cause a severe voltage drop in the network. Therefore, they inject high reactive power into the electrical network during these hours until a low voltage is obtained. Also, during 18: 00-22:00, as the amount of load consumption (active and reactive) in EHs is high, EVs also produce high reactive power to reduce the reactive power requested by EHs from the upstream network at these periods. Because of this type of performance of sources and ALs, EHs can feed reactive power into the electrical system during most operation hours.The expected daily heat curve of CHPs, boilers, DRPs, and EHs can be seen in Figure6(c). As seen in this figure, the daily heat power curve of CHPs has a trend of changes similar to the time diagram of the active power of CHPs in Figure 6(a); as given in Eq. (29), thermal power is a factor of active power in CHP. Boilers also inject heat power equal to their maximum capacity, that is, 0.9 MW (capacity of each boiler is 0.3 MW and three hubs contain boilers) during all operation hours. As per Figure 5(b), the price of thermal energy is always higher than the price of gas energy. Therefore, in order to minimize Eq. (24), the boilers inject high heat power into the EHs. DRPs are also in charging mode during the hours of cheap heat energy (1:00-4:00 and 16:00-00:00) and receive heat power from the EHs. Moreover, the opposite holds during hours with expensive heat energy (5:00-15:00). Also, heat boilers and DRPs are not switched off at all times in the presence of CHPs to improve the flexibility of EHs in the heat sector because CHPs generate heat power during all operation hours. Finally, due to the type of performance of the mentioned sources and ALs, EHs can inject heat power into retail markets at all operation hours according to Figure 6(c). The following is the expected daily gas power curve of CHP, boiler, and EHs in Figure 6(d). The gas power diagram of CHPs and boilers has the same trend of changes with their thermal power in Figure 6(c), which is confirmed in Eqs. (29), (30), and (33). Additionally, according to the data in subsection 4.1, the only consumers of the gas network are CHPs and boilers, so as in Figure 6(d), the gas power of EHs is equal to the gas power consumption of these sources at all simulation hours.The expected daily curve of energy costs of EHs in retail markets (CostEHs) is plotted in Figure 7. This curve can be calculated per hour of Eq. (24), providing that the sum is not used on the t subscript. Based on this relationship, calculating CostEHs requires the retail market price, λR, and the wholesale price, λW, as presented in subsection 4.1. It also requires the power curve of the EHs, which can also be seen in Figure 6. Eventually, according to the daily power curve of EHs in Figure 6 and energy price data in different markets following subsection 4.1, the expected daily curve of CostEHs is in Figure 7. As seen in this figure, EHs appear as consumers in retail markets at 1:00-7:00 and 23:00-00:00 due to energy consumption by ALs and low energy generation by sources. Therefore, they have a positive value for CostEHs. But at other times, the total power output of sources and ALs is greater than the power consumption of EHs, so EHs appear as producers in the retail markets (albeit they are energy consumers in electricity, heat, and gas networks). As a result, CostEHs will be negative during these hours. In other words, during these hours, EHs have income in the mentioned markets. Figure 8 depicts the CostEHs curve in terms of flexibility tolerance (εF). In this equation, the term CostEHs is calculated exactly from Eq. (24) and holds for all operation hours. Figure 8 shows that decreasing εF decreases CostEHs or increases EHs’ revenue because reducing εF means reducing the importance of EH flexibility based on Eqs. (42) and (43); hence, the solution space of the proposed problem increases. In this case, CHPs can be turned off from 5:00 to 7:00, because the price of electricity is lower than the price of gas energy. Also, electrical ALs are switched off during the mid-peak electricity consumption hours, 8: 00–16:00 and 23:00-00:00, when the price of electrical energy is about 32 $ /MWh (1.2×26.4). Charging operations are generally performed during off-peak hours with the lowest electricity prices. As a result, the energy cost of sources and ALs decreases, which means that EHs increase their revenue in retail markets. However, increasing εF causes ALs and functional sources to have the same function as in Figure 6, which is commensurate with the decrease in EHs’ revenue.Figure 7
Expected daily curves of EHs cost in retail markets with consideringεF = 0.02 MW.Figure 8
Curves of EHs cost in flexibility tolerance.Table 2
EHs profit in retail markets in different cases and load levels with consideringεF = 0.02 MW.
CaseLoad level (%)Calculation time (s)90100110I47.11143.15639.26237.3II51.06346.77842.56932.1III15.02513.72912.36832.2IV48.51244.38240.12532.1V45.61441.53137.37532.0VI53.42748.93744.27432.1Note that in Figure8, low values of flexibility tolerance mean that flexibility is included in the proposal. High values also mean a decrease in the importance of flexibility, which is close to the work of references [6–22] that have not observed flexibility.Finally, Table2 reports various capabilities of the proposed design. As per Table 2, the numerical results related to following six case studies are evaluated.Case I: the proposed scheme by considering a one-layer EMS strategy. In this case study, all sources, ESSs, and responsive loads are directly under the command of the ENOs. Hence, the problem model is described by Eqs. (1)–(43) other than Eq. (24). In other words, in this case, EH objectives are not considered because ENOs cannot sense them.Case II. The proposed design assuming the presence of a two-layer EMS. This section addresses problems (1)–(43).Case III. The same Case II assuming that the purchase and sale price of energy for EHs in the retail market are the same and it is equal to λR.Case IV. The same Case II assuming that the purchase and sale price of energy for EHs in the retail market are fixed at all hours of operation and they do not change over time. In this case, the purchase price of electricity, heat, and gas are 32 $/MWh, 30 $/MWh, and 18 $/MWh, respectively, and the purchase price of energy is 20% lower than the purchase price.Case V. Case IV, assuming a 10% reduction in the purchase and sale price of energy in the retail market.Case VI. Case IV, assuming a 20% increase in the purchase and sale price of energy in the retail market.In the following, according to Table2, the proposed design for Case II has been able to obtain a higher profit than Case I. Also, the calculation time of the proposed design is less than Case I. This is the advantage of indirect management of ALs and sources by ENOs considering EHOs as intermediaries over their direct management by ENOs. In addition, considering the different purchase and sale prices of EH energy in retail markets in Case II in comparison with Case III, higher income was obtained in Case II for EHs. In Case III, since the price of sale and purchase is the same in the retail market and it generally has a higher value than the energy price in the wholesale market, ENs are reluctant to buy energy from EHs. Unless they buy energy from EHs in critical situations such as peak hours when there is a possibility of voltage drop and high heat and pressure. Hence, the income of EHs is the lowest in Case III, as in Table 2. In cases IV to VI, it is observed that for high energy prices (Case VI), it is possible to earn more income for EHs than in Case II. This, however, has a favorable social impact on EHs. Moreover, it has low social benefits for consumers because they have to buy energy at a high price. Therefore, in terms of considering the social welfare of producers and consumers at the same time, according to Table 2, energy prices should change over time. Finally, as per Table 2, the computational time in designs with two-layer EMS is less than that in those with one-layer EMS.B)Assessing the economic situation and operation of energy networks: In Table 3, the Pareto front for the suggested plan is presented. As this table shows, the minimum and maximum energy costs of ENs (CostENs) are $5995 and $9712.5, respectively. Therefore, the range of its changes is equal to $3737.9 (9712.5–5995.6). The minimum and maximum expected energy losses of the ENs (EEL) are 5.29 MWh and 8.23 MWh, respectively. Following this, the range of changes of EEL is 2.94 MWh (28.28-5.29). Hence, εEEL will be between 5.29 and 8.23, and in Table 2, it changes with steps of 0.21. According to Table 2, the direction of change of CostENs and EEL functions is not the same, so that the decrease in EEL is proportional to the increase in CostENs because to reduce CostENs, as in Eq. (1), EHs need to be energy-consuming, but to reduce EEL, local sources such as EHs need to provide a significant portion of the grid energy consumption. In other words, to reduce EEL, EHs must be energy producers. The results of fuzzy decision making are presented in columns 4–6 of Table 2 as in subsection 3.1. Finally, based on the fuzzy decision-making technique, the best compromise point between EEL and CostENs is a point at which εEEL is 5.92, which is shown in Table 2. At this point, the optimal values of EEL and CostENs functions are 5.92 MWh and $65.293, respectively. Accordingly, the EEL function is about 21.4% (5.92-5.29) ÷ (8.23-5.29) away from its minimum value. The CostENs function is also about 16% away from its minimum value. Therefore, the proposed scheme finds values for EEL and CostENs functions at the compromise point, which are close to their minimum values.Table 3
Pareto front and best compromise of the suggested method forεF = 0.02 MW.
Pareto FrontεEELCostENs ($)EEL (MWh)f^CostENsf^EELminf^CostENs,f^EEL8.235995.68.231008.025997.28.020.99960.07140.07147.816001.57.810.99840.14290.14297.606010.67.600.99600.21430.21437.396024.37.390.99230.28570.28577.186050.17.180.98530.35710.35716.976096.56.970.97290.42860.42866.766151.36.760.95810.50000.50006.556216.86.550.94050.57140.57146.346304.76.340.91680.64290.64296.136433.66.130.88220.71430.71435.926593.25.920.83920.78570.78575.717371.95.710.62970.85710.62975.508432.15.500.34450.92860.34455.299712.55.29010Best Compromise SolutionParameterUT (33 scenario)SBSP with 60 scenarioSBSP with 80 scenarioSBSP with 100 scenarioCostENs ($)6593.26632.16593.76593.3EEL (MWh)5.925.675.905.91Calculation time (s)32.141.646.151.5In Table3, the compromise point for uncertainty modeling based on UT and SBSP is reported. In SBSP, Monte Carlo simulations generate a large number of scenarios. The selection of the values of each uncertainty parameter in each scenario is based on their mean and standard deviation values. Then, the probability of each parameter is calculated. The probability of load, energy price, wind power, solar power, and EV parameters is calculated by the normal, normal, Weibull, beta, and Rayleigh probability distribution functions, respectively. Then, the probability of each scenario is equal to the product of the probability values of these parameters in this scenario. Finally, the Kantorovich method, as a scenario reduction technique, applies a small number of generated scenarios to the problem, so that their distances are close to each other. Details of this method are provided in Reference [14]. Based on the results of Table 3, it can be seen that the results of UT and SBSP are close to each other if more than 80 scenarios are considered for SBSP, while UT needs 33 scenarios. This has led to an increase in SBSP computing time compared with UT.Figure9 illustrates the expected daily curve of power through electricity, gas, and heat substations, as well as the cost of purchasing energy of all ENs, mentioned in the wholesale and retail energy markets for εF = 0.02 MW and εEEL = 5.92. According to Figure 9(a), high active power passes through the electrical substation between 1:00 and 7:00 because, according to Figure 5, most of the electrical ALs are in discharging mode and EHs are consuming electrical energy during these hours. At other hours, the EHs act similar to electricity generators as shown in Figure 6, in which case, active power of less than MW 2 passes through the substation during most operation hours. According to Figure 9(a), reactive power greater than 1.5 MVAr passes through the substation at 10:00 to 16:00, while at other times it is less than 1.5 MVAr. Because according to Figure 6(b), EHs consume reactive power during 10:00-16:00, but at other times they appear as generators of this power in the electrical network. According to [14], the peak heat load is equal to 4.3 MW. If in Figure 9, the thermal power passing through the substation in Case II is around 3.5 MW, and in most operation hours, it is less than 2.5 MW. This is because according to Figure 6(c), EHs act always as heat energy generators that inject heat into the heat network. Hence, the power demands of this network from the upstream network decrease. Finally, the expected daily gas power curve passing through the gas substation has the same trend as the daily gas power curve of EHs in Figure 6(d) because, according to subsection 4.1, the only consumers of the gas network are EHs with CHP and boiler. The expected daily energy cost curve of ENs is illustrated in Figure 9(b). This curve is calculated from Eq. (1), providing that the sum is not considered on subscript t. In the following, considering this figure, the mentioned cost is the lowest during 1:00-7:00 and 23:00-00:00, but they have high amounts in other hours. This is because the price of energy is generally lower during 1:00-7:00 and 23:00-00:00 than other hours as is seen in Figure 5(b).Figure 9
Expected daily curves of (a) station power and (b) ENs cost forεF = 0.02 MW and εEEL = 5.92.
(a)(b)Finally, in Table4, the values of economic and operation indices of ENs for flexibility tolerance of 0.02 MW and the best compromise point (εEEL = 5.92) in cases II and VII include power flow studies of ENs. Operation indices include EEL in different ENs; the maximum drop of voltage, pressure, and temperature; and maximum overvoltage, over-pressure, and over-temperature. In the context of economic indices, the phrase CostENs is also presented. According to this table, the energy cost of networks in power flow studies of ENs is roughly $7184.8, while the proposed plan for proper management of EHs has been able to reduce it to $65.293. In other words, Case II has improved the economic status of ENs by approximately 8.2% (7184.8-6593.2/7184.2)) in comparison with power flow studies. From operation angle, there is no overvoltage, overpressure, and overtemperature in Case VII for ENs. Also, since in this case gas energy consumers (EHs with CHP and boiler) are not present, EEL and pressure drop for the gas network are equal to zero. But in Case VII, the maximum voltage drop and temperature are 0.112 and 0.116 p.u., respectively, while their maximum allowable limit is 0.1 p.u. (1-0.9). The electricity and heat networks in the mentioned studies have EELs greater than 3.7 MWh and 2.8 MWh, respectively. In Case II, although the energy losses in the gas network increase to 1.47 MWh, the EEL in the electricity and heat networks decreased compared to Case VII by about 34.65% ((3.78-2.47)/3.78) and 30% ((2.83-1.98)/2.83). This has reduced the total EEL in these ENs from 6.61 MWh in Case VII to 5.92 MWh in Case II. Therefore, the suggested method can reduce the EEL by about 10.5% in comparison with power flow studies. In Case II, although the maximum overvoltage and overtemperature have increased to approximately 0.01 p.u. and the maximum pressure drop has increased to 0.039 p.u., the maximum voltage and temperature drop are 0.051 and 078, respectively. Consequently, the proposed design reduced the maximum voltage and temperature drop by about 54.5% and 32.8% compared with Case VII, respectively.Table 4
Value of economic and operation indices forεF = 0.02 MW and εEEL = 5.92 in different cases.
CaseIIVIICostENs ($)6593.27184.8EEL (MWh) in electricity network2.473.78EEL (MWh) in heating network1.982.83EEL (MWh) in gas network1.470Total EEL (MWh)5.926.61Maximum voltage drop (p.u.)0.0510.112Maximum temperature drop (p.u.)0.0780.116Maximum pressure drop (p.u.)0.0390Maximum overvoltage (p.u.)0.0110Maximum overtemperature (p.u.)0.0090Maximum overpressure (p.u.)00


## 4.1. Problem Data

Figure4 demonstrates the application of the proposed scheme to a test system consisting of a 9-bus electrical network, a 4-node gas network, and a 7-node heating network [14]. The base power in the electrical network is 1 MVA, and it is 1 MW in the gas and heating networks [14]. Moreover, the base voltage, pressure, and temperature are 1 kV, 10 Bar, and 100°C, and their allowable limit is [0.9, 1.1] p.u. [13, 27, 28]. The specifications of distribution lines and gas and heat pipelines are presented in Reference [14]. Peak electrical and thermal data are reported in Reference [14]. In the gas network, CHP and boiler are the only consumers, and the passive gas load is zero. Hourly load data are found by multiplying peak load and load factor [12, 29, 30]. The expected daily load factor curve for electrical and heating networks is shown in Figure 5(a) [14]. The expected daily energy price curve in the wholesale market for various ENs is shown in Figure 5(b) [14]. ENs are assumed to increase the energy price in the retail market by 20% compared with the wholesale market to make a profit, that is, λR = 1.2×λW.Figure 4
Test system [14].Figure 5
Expected daily curves of (a) load factor [14], (b) energy price in Wholesale market [14], (c) RES power rate [31], and (d) EVs penetration rate [32].
(a)(b)(c)(d)The network includes 7 EHs whose locations in different ENs are specified in Figure4, and their load is presented in Reference [14]. EHs 1–3 and 5 have electrical sources and ALs, that is, RESs, EVs, batteries, and electrical responsive loads. The EH 4 also consists of a CHP, a boiler, and heat-responsive loads. Hubs 6 and 7 also contain all of these elements. CHP has a maximum electrical and heat output of 1 MVA and 1 MW, respectively, with ηT, ηL, and ηH being 40%, 9%, and 40%, respectively [23]. The boiler with 80% efficiency also has a maximum heat output of 0.3 MW [14]. RESs including a 0.25 MVA photovoltaic and a 0.2 MVA wind system are used in EHs [14]. The hourly power of RES is calculated by multiplying its capacity and the rate of power output [13]. The expected daily curve of power output for photovoltaics and wind systems is illustrated in Figure 5(c) [31]. The battery with a charging and discharging efficiency of 88% and capacity of 2 MWh is installed in hubs 1–3 and 5–7 [14]. The battery has a charging and discharging rate of 0.5 MW and charger capacity of 0.6 MVA, its minimum storable energy, and initial energy is 0.2 MWh [14]. In each of the EHs 1–3 and 5–7, 60 EVs can be connected, the specifications of each EV including charging/discharging rate, charger capacity, and other items are presented in Reference [1,2]. The number of EVs connected to the EH per hour is proportional to the product of the total number of EVs in the parking lot and the penetration rate of EV into the parking lot [32]. Figure 5(d) plots the expected daily curve of EVs penetration rate [32]. Hubs are assumed to take part in the DRP scheme with a penetration rate of 30% [14]. To achieve the desired flexibility in EH, a flexibility tolerance of 0.02 MW is considered.


## 4.2. Results

In this subsection, the proposed scheme applied to the data given in subsection4.1 is simulated using GAMS, and then, the CPLEX solves the problem [33]. The following are the details of the numerical results obtained from various case studies. It is noteworthy that here 5 linear pieces were used in the conventional piecewise linearization method, and the circular plane is approximated to a regular 45-gon. According to these conditions and regarding the results obtained in [14], the computational error for active, reactive, gas, and heat power in the linear approximation model used compared with the original nonlinear model is about 2%, 2%, 0.9%, and 0, respectively. Also, these values for voltage, pressure, and temperature are roughly 0.5%, 0.1%, and 0, respectively.A)Performance evaluation of EHs: The expected daily curve of active, reactive, heat, and gas power of EHs and their sources and ALs is plotted in Figure 6. Based on Figure 6(a) and according to Figure 5(c) and the data of subsection 4.1, RESs such as photovoltaics and wind systems inject active power equal to the maximum active power into EHs in all operation hours according to the weather conditions. As their operating costs are negligible (zero in this paper), EHs tend to make high use of the active power generation of RESs in order to minimize their energy costs in the retail market, as in Eq. (24). In the case of CHP, it is observed that they feed higher amounts of active power into EHs during the hours when the price of electric energy is higher than the price of gas energy according to Figure 5(b) so that at 1:00–4:00, 8:00–16:00, and 23:00-00:00, feed active power is equal to their maximum capacity, that is, 3 MW (capacity of each CHP is 1 MVA and three EH contain CHPs) into EHs. Also, according to Figure 5(a), there is peak electric load during 17:00-22:00, that is, active and reactive power consumed during these hours is high. Hence, the voltage drop may be high at this period. Therefore, in addition to injecting active power, CHPs also tend to inject reactive power during these hours, so the active power produced by CHPs is less than their maximum capacity. Also, according to Figure 5(b), at 5:00-7:00, the price of gas energy is higher than the price of electricity. Therefore, CHPs produce low active power during these hours to minimize the energy cost of EHs, but they are switched off. The reason why these sources are not turned off during the mentioned hours can be due to the following:CHPs produce thermal power at their output as well, the amount of which depends on the active power of CHPs. As shown in Figure5(b), the price of heat energy is higher than the price of gas energy at 5:00-7:00. Therefore, considering Eq. (24), CHP tends to produce heat power, which requires active power generation based on Eq. (29).CHPs are a flexible source for RESs. As shown in Figure5(c), RESs generate active power at all hours of operation. Therefore, it is expected that to improve the flexibility of EHs, the flexibility sources must control their active power at all hours of operation, so they will always be on at all hours.Figure 6
Expected daily curves of (a) active, (b) reactive, (c) heating, and (d) gas power of EHs and their sources and ALs with consideringεF = 0.02 MW.
(a)(b)(c)(d)ALs, such as batteries, EVs, and DRPs, perform charging operations when electricity price is low, that is, between 1:00-16:00 and 23:00-00:00, as shown in Figure6(a). During the peak load period when the price of electrical energy has the highest value as per Figure 5(b), ALs perform the discharging operation and inject active power into EHs. Moreover, they are switched on at all hours of operation, except for EVs, which do not inject power into EHs between 10:00 and 11:00 because, according to Figure 5(d), no EVs are connected to the EHs during these hours. This is the function model of ALs to enhance the flexibility of EHs with RESs. Finally, based on Eq. (25), the expected daily curve of the active power of EHs is like the blue curve in Figure 6(a). They appear as consumers of electrical energy in ENs because of high power consumption by ALs and lower power generation by sources during 1:00-8:00. But at other times, they are a generator of electrical energy.In the field of the expected daily curve of reactive power of EHs, it can be seen according to Figure6(b) that RESs and batteries have a time graph similar to the daily curve of reactive charge. Only at 12:00 do photovoltaics inject high active power into EHs as Figure 6(a), so their injected active power into EHs is reduced to zero. Such conditions exist for wind systems at 17:00. CHPs also inject reactive power into EHs only during peak electricity consumption interval, 17:00-22:00, so that it contributes to supplying high reactive power consumption of EHs. EVs also inject high reactive power into EHs between 1:00-7:00 and 18:00-00:00 because during these periods; according to Figure 6(a), EVs receive high active power from the electrical network, which may cause a severe voltage drop in the network. Therefore, they inject high reactive power into the electrical network during these hours until a low voltage is obtained. Also, during 18: 00-22:00, as the amount of load consumption (active and reactive) in EHs is high, EVs also produce high reactive power to reduce the reactive power requested by EHs from the upstream network at these periods. Because of this type of performance of sources and ALs, EHs can feed reactive power into the electrical system during most operation hours.The expected daily heat curve of CHPs, boilers, DRPs, and EHs can be seen in Figure6(c). As seen in this figure, the daily heat power curve of CHPs has a trend of changes similar to the time diagram of the active power of CHPs in Figure 6(a); as given in Eq. (29), thermal power is a factor of active power in CHP. Boilers also inject heat power equal to their maximum capacity, that is, 0.9 MW (capacity of each boiler is 0.3 MW and three hubs contain boilers) during all operation hours. As per Figure 5(b), the price of thermal energy is always higher than the price of gas energy. Therefore, in order to minimize Eq. (24), the boilers inject high heat power into the EHs. DRPs are also in charging mode during the hours of cheap heat energy (1:00-4:00 and 16:00-00:00) and receive heat power from the EHs. Moreover, the opposite holds during hours with expensive heat energy (5:00-15:00). Also, heat boilers and DRPs are not switched off at all times in the presence of CHPs to improve the flexibility of EHs in the heat sector because CHPs generate heat power during all operation hours. Finally, due to the type of performance of the mentioned sources and ALs, EHs can inject heat power into retail markets at all operation hours according to Figure 6(c). The following is the expected daily gas power curve of CHP, boiler, and EHs in Figure 6(d). The gas power diagram of CHPs and boilers has the same trend of changes with their thermal power in Figure 6(c), which is confirmed in Eqs. (29), (30), and (33). Additionally, according to the data in subsection 4.1, the only consumers of the gas network are CHPs and boilers, so as in Figure 6(d), the gas power of EHs is equal to the gas power consumption of these sources at all simulation hours.The expected daily curve of energy costs of EHs in retail markets (CostEHs) is plotted in Figure 7. This curve can be calculated per hour of Eq. (24), providing that the sum is not used on the t subscript. Based on this relationship, calculating CostEHs requires the retail market price, λR, and the wholesale price, λW, as presented in subsection 4.1. It also requires the power curve of the EHs, which can also be seen in Figure 6. Eventually, according to the daily power curve of EHs in Figure 6 and energy price data in different markets following subsection 4.1, the expected daily curve of CostEHs is in Figure 7. As seen in this figure, EHs appear as consumers in retail markets at 1:00-7:00 and 23:00-00:00 due to energy consumption by ALs and low energy generation by sources. Therefore, they have a positive value for CostEHs. But at other times, the total power output of sources and ALs is greater than the power consumption of EHs, so EHs appear as producers in the retail markets (albeit they are energy consumers in electricity, heat, and gas networks). As a result, CostEHs will be negative during these hours. In other words, during these hours, EHs have income in the mentioned markets. Figure 8 depicts the CostEHs curve in terms of flexibility tolerance (εF). In this equation, the term CostEHs is calculated exactly from Eq. (24) and holds for all operation hours. Figure 8 shows that decreasing εF decreases CostEHs or increases EHs’ revenue because reducing εF means reducing the importance of EH flexibility based on Eqs. (42) and (43); hence, the solution space of the proposed problem increases. In this case, CHPs can be turned off from 5:00 to 7:00, because the price of electricity is lower than the price of gas energy. Also, electrical ALs are switched off during the mid-peak electricity consumption hours, 8: 00–16:00 and 23:00-00:00, when the price of electrical energy is about 32 $ /MWh (1.2×26.4). Charging operations are generally performed during off-peak hours with the lowest electricity prices. As a result, the energy cost of sources and ALs decreases, which means that EHs increase their revenue in retail markets. However, increasing εF causes ALs and functional sources to have the same function as in Figure 6, which is commensurate with the decrease in EHs’ revenue.Figure 7
Expected daily curves of EHs cost in retail markets with consideringεF = 0.02 MW.Figure 8
Curves of EHs cost in flexibility tolerance.Table 2
EHs profit in retail markets in different cases and load levels with consideringεF = 0.02 MW.
CaseLoad level (%)Calculation time (s)90100110I47.11143.15639.26237.3II51.06346.77842.56932.1III15.02513.72912.36832.2IV48.51244.38240.12532.1V45.61441.53137.37532.0VI53.42748.93744.27432.1Note that in Figure8, low values of flexibility tolerance mean that flexibility is included in the proposal. High values also mean a decrease in the importance of flexibility, which is close to the work of references [6–22] that have not observed flexibility.Finally, Table2 reports various capabilities of the proposed design. As per Table 2, the numerical results related to following six case studies are evaluated.Case I: the proposed scheme by considering a one-layer EMS strategy. In this case study, all sources, ESSs, and responsive loads are directly under the command of the ENOs. Hence, the problem model is described by Eqs. (1)–(43) other than Eq. (24). In other words, in this case, EH objectives are not considered because ENOs cannot sense them.Case II. The proposed design assuming the presence of a two-layer EMS. This section addresses problems (1)–(43).Case III. The same Case II assuming that the purchase and sale price of energy for EHs in the retail market are the same and it is equal to λR.Case IV. The same Case II assuming that the purchase and sale price of energy for EHs in the retail market are fixed at all hours of operation and they do not change over time. In this case, the purchase price of electricity, heat, and gas are 32 $/MWh, 30 $/MWh, and 18 $/MWh, respectively, and the purchase price of energy is 20% lower than the purchase price.Case V. Case IV, assuming a 10% reduction in the purchase and sale price of energy in the retail market.Case VI. Case IV, assuming a 20% increase in the purchase and sale price of energy in the retail market.In the following, according to Table2, the proposed design for Case II has been able to obtain a higher profit than Case I. Also, the calculation time of the proposed design is less than Case I. This is the advantage of indirect management of ALs and sources by ENOs considering EHOs as intermediaries over their direct management by ENOs. In addition, considering the different purchase and sale prices of EH energy in retail markets in Case II in comparison with Case III, higher income was obtained in Case II for EHs. In Case III, since the price of sale and purchase is the same in the retail market and it generally has a higher value than the energy price in the wholesale market, ENs are reluctant to buy energy from EHs. Unless they buy energy from EHs in critical situations such as peak hours when there is a possibility of voltage drop and high heat and pressure. Hence, the income of EHs is the lowest in Case III, as in Table 2. In cases IV to VI, it is observed that for high energy prices (Case VI), it is possible to earn more income for EHs than in Case II. This, however, has a favorable social impact on EHs. Moreover, it has low social benefits for consumers because they have to buy energy at a high price. Therefore, in terms of considering the social welfare of producers and consumers at the same time, according to Table 2, energy prices should change over time. Finally, as per Table 2, the computational time in designs with two-layer EMS is less than that in those with one-layer EMS.B)Assessing the economic situation and operation of energy networks: In Table 3, the Pareto front for the suggested plan is presented. As this table shows, the minimum and maximum energy costs of ENs (CostENs) are $5995 and $9712.5, respectively. Therefore, the range of its changes is equal to $3737.9 (9712.5–5995.6). The minimum and maximum expected energy losses of the ENs (EEL) are 5.29 MWh and 8.23 MWh, respectively. Following this, the range of changes of EEL is 2.94 MWh (28.28-5.29). Hence, εEEL will be between 5.29 and 8.23, and in Table 2, it changes with steps of 0.21. According to Table 2, the direction of change of CostENs and EEL functions is not the same, so that the decrease in EEL is proportional to the increase in CostENs because to reduce CostENs, as in Eq. (1), EHs need to be energy-consuming, but to reduce EEL, local sources such as EHs need to provide a significant portion of the grid energy consumption. In other words, to reduce EEL, EHs must be energy producers. The results of fuzzy decision making are presented in columns 4–6 of Table 2 as in subsection 3.1. Finally, based on the fuzzy decision-making technique, the best compromise point between EEL and CostENs is a point at which εEEL is 5.92, which is shown in Table 2. At this point, the optimal values of EEL and CostENs functions are 5.92 MWh and $65.293, respectively. Accordingly, the EEL function is about 21.4% (5.92-5.29) ÷ (8.23-5.29) away from its minimum value. The CostENs function is also about 16% away from its minimum value. Therefore, the proposed scheme finds values for EEL and CostENs functions at the compromise point, which are close to their minimum values.Table 3
Pareto front and best compromise of the suggested method forεF = 0.02 MW.
Pareto FrontεEELCostENs ($)EEL (MWh)f^CostENsf^EELminf^CostENs,f^EEL8.235995.68.231008.025997.28.020.99960.07140.07147.816001.57.810.99840.14290.14297.606010.67.600.99600.21430.21437.396024.37.390.99230.28570.28577.186050.17.180.98530.35710.35716.976096.56.970.97290.42860.42866.766151.36.760.95810.50000.50006.556216.86.550.94050.57140.57146.346304.76.340.91680.64290.64296.136433.66.130.88220.71430.71435.926593.25.920.83920.78570.78575.717371.95.710.62970.85710.62975.508432.15.500.34450.92860.34455.299712.55.29010Best Compromise SolutionParameterUT (33 scenario)SBSP with 60 scenarioSBSP with 80 scenarioSBSP with 100 scenarioCostENs ($)6593.26632.16593.76593.3EEL (MWh)5.925.675.905.91Calculation time (s)32.141.646.151.5In Table3, the compromise point for uncertainty modeling based on UT and SBSP is reported. In SBSP, Monte Carlo simulations generate a large number of scenarios. The selection of the values of each uncertainty parameter in each scenario is based on their mean and standard deviation values. Then, the probability of each parameter is calculated. The probability of load, energy price, wind power, solar power, and EV parameters is calculated by the normal, normal, Weibull, beta, and Rayleigh probability distribution functions, respectively. Then, the probability of each scenario is equal to the product of the probability values of these parameters in this scenario. Finally, the Kantorovich method, as a scenario reduction technique, applies a small number of generated scenarios to the problem, so that their distances are close to each other. Details of this method are provided in Reference [14]. Based on the results of Table 3, it can be seen that the results of UT and SBSP are close to each other if more than 80 scenarios are considered for SBSP, while UT needs 33 scenarios. This has led to an increase in SBSP computing time compared with UT.Figure9 illustrates the expected daily curve of power through electricity, gas, and heat substations, as well as the cost of purchasing energy of all ENs, mentioned in the wholesale and retail energy markets for εF = 0.02 MW and εEEL = 5.92. According to Figure 9(a), high active power passes through the electrical substation between 1:00 and 7:00 because, according to Figure 5, most of the electrical ALs are in discharging mode and EHs are consuming electrical energy during these hours. At other hours, the EHs act similar to electricity generators as shown in Figure 6, in which case, active power of less than MW 2 passes through the substation during most operation hours. According to Figure 9(a), reactive power greater than 1.5 MVAr passes through the substation at 10:00 to 16:00, while at other times it is less than 1.5 MVAr. Because according to Figure 6(b), EHs consume reactive power during 10:00-16:00, but at other times they appear as generators of this power in the electrical network. According to [14], the peak heat load is equal to 4.3 MW. If in Figure 9, the thermal power passing through the substation in Case II is around 3.5 MW, and in most operation hours, it is less than 2.5 MW. This is because according to Figure 6(c), EHs act always as heat energy generators that inject heat into the heat network. Hence, the power demands of this network from the upstream network decrease. Finally, the expected daily gas power curve passing through the gas substation has the same trend as the daily gas power curve of EHs in Figure 6(d) because, according to subsection 4.1, the only consumers of the gas network are EHs with CHP and boiler. The expected daily energy cost curve of ENs is illustrated in Figure 9(b). This curve is calculated from Eq. (1), providing that the sum is not considered on subscript t. In the following, considering this figure, the mentioned cost is the lowest during 1:00-7:00 and 23:00-00:00, but they have high amounts in other hours. This is because the price of energy is generally lower during 1:00-7:00 and 23:00-00:00 than other hours as is seen in Figure 5(b).Figure 9
Expected daily curves of (a) station power and (b) ENs cost forεF = 0.02 MW and εEEL = 5.92.
(a)(b)Finally, in Table4, the values of economic and operation indices of ENs for flexibility tolerance of 0.02 MW and the best compromise point (εEEL = 5.92) in cases II and VII include power flow studies of ENs. Operation indices include EEL in different ENs; the maximum drop of voltage, pressure, and temperature; and maximum overvoltage, over-pressure, and over-temperature. In the context of economic indices, the phrase CostENs is also presented. According to this table, the energy cost of networks in power flow studies of ENs is roughly $7184.8, while the proposed plan for proper management of EHs has been able to reduce it to $65.293. In other words, Case II has improved the economic status of ENs by approximately 8.2% (7184.8-6593.2/7184.2)) in comparison with power flow studies. From operation angle, there is no overvoltage, overpressure, and overtemperature in Case VII for ENs. Also, since in this case gas energy consumers (EHs with CHP and boiler) are not present, EEL and pressure drop for the gas network are equal to zero. But in Case VII, the maximum voltage drop and temperature are 0.112 and 0.116 p.u., respectively, while their maximum allowable limit is 0.1 p.u. (1-0.9). The electricity and heat networks in the mentioned studies have EELs greater than 3.7 MWh and 2.8 MWh, respectively. In Case II, although the energy losses in the gas network increase to 1.47 MWh, the EEL in the electricity and heat networks decreased compared to Case VII by about 34.65% ((3.78-2.47)/3.78) and 30% ((2.83-1.98)/2.83). This has reduced the total EEL in these ENs from 6.61 MWh in Case VII to 5.92 MWh in Case II. Therefore, the suggested method can reduce the EEL by about 10.5% in comparison with power flow studies. In Case II, although the maximum overvoltage and overtemperature have increased to approximately 0.01 p.u. and the maximum pressure drop has increased to 0.039 p.u., the maximum voltage and temperature drop are 0.051 and 078, respectively. Consequently, the proposed design reduced the maximum voltage and temperature drop by about 54.5% and 32.8% compared with Case VII, respectively.Table 4
Value of economic and operation indices forεF = 0.02 MW and εEEL = 5.92 in different cases.
CaseIIVIICostENs ($)6593.27184.8EEL (MWh) in electricity network2.473.78EEL (MWh) in heating network1.982.83EEL (MWh) in gas network1.470Total EEL (MWh)5.926.61Maximum voltage drop (p.u.)0.0510.112Maximum temperature drop (p.u.)0.0780.116Maximum pressure drop (p.u.)0.0390Maximum overvoltage (p.u.)0.0110Maximum overtemperature (p.u.)0.0090Maximum overpressure (p.u.)00


## 5. Conclusion

This paper presents the simultaneous participation of energy networks and energy hubs in day-ahead energy markets. ENs buy energy from the wholesale energy market and make it available to EHs and consumers in the retail energy markets. The proposed design was based on two-layer EMS, which, in the first layer of EMS, considered the bilateral coordination of sources and ALs with the EH operator. Also, the coordination of these operators with ENOs was considered in the second layer of EMS. Following this, the proposed plan was expressed as a bilevel optimization, the upper level of which aims to find the minimum expected energy cost of ENs in the mentioned markets and minimize the expected energy losses of ENs as a two-objective problem subject to LOPF equations. The lower-level formulation also considered minimizing the expected energy cost of EHs (equal to the difference between energy sale and purchase) in the retail energy markets, which is bound by the operation model of sources and ALs and the flexibility limitation of EHs. Then, KKT was used to extract the single-level problem, and theε-constraint-based Pareto optimization was adopted to find a single-objective problem. UT was then employed to develop modeling of the uncertainties of load, energy prices, renewable power, and energy demand of EVs. In the end, considering the obtained numerical results, the proposed design following the two-layer EMS strategy obtains the highest profit for EHs very quickly compared to the case with one-layer EMS. It has also reduced the energy cost of ENs by about 8% in comparison with power flow studies by optimal energy management of EHs. By adopting optimal scheduling for EHs, it also reduces energy loss of ENs, maximum voltage, and temperature drop by about 10%, 36%, and 30%, respectively, compared with power flow studies.∠P:A set related to the piecewise in the conventional piecewise linearization method



---

*Source: 2463003-2022-08-25.xml*