# Participation of Energy Storage-Based Flexible Hubs in Day-Ahead Reserve Regulation and Energy Markets Based on a Coordinated Energy Management Strategy

Authors: Mohsen Kazemi; Saeed Yeganeh Salehpour; Farid Shahbaazy; Saeed Behzadpoor; Sasan Pirouzi; Saeid Jafarpour
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

In order to increase energy efficiency, the energy hub is considered as a form of aggregator and coordinator of various resources and storage. With the optimal performance of resources and storage generators based on a proper energy management system, it is expected that hubs can gain financial benefits from energy markets and ancillary services. So, the paper presents the participation of networked energy hubs in day-ahead (DA) reserve regulation and energy markets, where the hub operator incorporates a coordinated energy management (CEM) strategy to manage power sources and energy storage devices within the hub. Hence, this problem maximizes the total profit of hubs in the DA energy and up and down reserve markets. Also, the problem is constrained by optimal power flow (OPF) constraints in gas, electricity, and thermal networks, reserve limits, and hub constraints, including the model of the combined heat and power (CHP), renewable energy source (RES), electrical/thermal storage, parking lots of electric vehicles (EVs), and boiler. Following that, a linear format is obtained for the nonlinear equation using traditional linearization methods so that an optimal solution is found in less time considering less computational error. Eventually, a standard case system is used to test the strategy, and thus, the capabilities of the approach are investigated. The obtained findings validate the potential of the proposed design in enhancing the economic situation of power sources and storage in hub form, which can enhance operation indices by optimal management of the hub so that the energy management of resources and storage in the form of a hub based on CEM compared to their independent management plan has been able to increase the profit of these elements in energy and up and down reserve markets by about 17%, 28%, and 15%, respectively. Regarding technical indices of energy networks, the proposed scheme by creating low energy losses in the gas network and providing pressure drop, overvoltage, and overtemperature within their permissible limits succeeded in reducing the energy losses in electricity and heat networks by about 83% and 38%, respectively, compared to power flow studies. Also, in these conditions, it has reduced the maximum voltage and temperature drop by 45% and 39%, respectively.


---

## Body

## 1. Introduction

In the present era, environmentally friendly energy sources such as wind, solar, microturbine, and CHP systems are utilized at the consumption points to lessen the emissions caused by fossil fuel consumption [1–5]. Moreover, it is anticipated that ESSs such as batteries and EVs along with DRPs with energy management can be effective in reducing pollution [6]. Therefore, it is expected that the number of sources and ESSs in energy networks will increase [7, 8]. Consequently, the network management by its operator may be challenged due to the increase in data volume. To address this, the smart grid proposes the use of an aggregator and coordinator framework for sources and ESSs so that these elements are managed by the aggregator and that these operators are in bidirectional coordination with the network operator [9]. As another point, there are different types of consuming energies, including electricity, natural gas, and district heating (EGDH), at consumption points such as residential, commercial, and industrial loads, which are generally managed and operated separately. Nonetheless, it is noted that simultaneous management of all energy sources can provide reliability, security, flexibility, and economic benefits to a system that contains various networks from an energy perspective [10]. According to these issues, the energy hub theory was introduced in [11]. It should be said that the hub is a modern framework that contains a simplified model of a system able to receive, transfer, transform, and store energy with the help of power sources and storage [12, 13]. Therefore, it is predicted that utilizing a CEM strategy for different power sources and storage in the hub can lead to high benefits for these devices in different energy markets and networks.There is huge research regarding hub energy management. The nonlinear optimal power flow models of operation of the energy hub are given in [14, 15], respectively, where evolutionary algorithms (EAs), including genetic algorithm (GA), have been used to provide solutions to the problems. Moreover, the authors in [16–18] employ EAs such as Particle Swarm Optimization (PSO) algorithm, evolutionary PSO solver, and the Fruit Fly Optimization (FFO) algorithm to solve the hub operation problem in electricity, gas, and district heating network because the operation model of hubs is generally a nonlinear programming (NLP) or mixed-integer nonlinear programming (MINLP) problem [19]. The operation model of hubs in different energy networks has widely been discussed in the literature. The authors in [20] propose a robust operation method for a grid-connected hub, including CHP and electric vehicles (EVs) parking lot, in which the hub operation cost is considered as the objective function. Also, a multistep approach is used in [21] to model the large-scale energy hub operation. The economic dispatching model of the hub considering minimization of the hub cost is presented in [22]. Besides, some research expresses a market model for hubs, such as in [23], where the day-ahead (DA) electrical energy market method is used for grid-connected hubs, including RESs, an electrical storage system (ESS), and CHP. Also, the DA energy market for electrical, gas, and heating energies is presented in [24] for networked hubs. An EH containing various renewables is discussed in [25]. The paper provides modeling of this system while modeling the operation and planning of the EH and the intermittent nature of PV and wind energy generating units. The main issue with the operation of the discussed system is how to find the optimal synergy of power sources when meeting the demand. Quantum Particle Swarm Optimization (QPSO) is adopted in the paper to solve the problem and find the minimum overall cost burden of the EH. Thermal power is utilized in [26] to produce electricity by adopting an Organic Rankine Cycle (ORC). Also, a seawater desalination (SWD) unit with Reverse Osmosis (RO) technology is used as a flexible load. The research also incorporates one parking lot and one battery swapping station (BSS) to provide a more suitable interaction among the EVs. Hydrogen is supplied using an electrolyzer. The paper uses a hydrogen tank and a fuel cell to make a balance between hydrogen generation and consumption. Uncertainties are modeled using a scenario generation approach and an information gap decision theory (IGDT) method. In another study, the operation of smart EHs is simulated using a new model based on Monte Carlo and a computational approach. By adopting such methods, [27] introduces a holistic evaluation of technical and economic feasibility research on smart EHs. Energy management of EHs connected to electricity, gas, and heating networks is discussed in [28]. The EHs play the role of a coordinator among distributed generation units and storage devices. Using a deterministic model, the overall operation cost of EHs is minimized while taking into account the equations of optimal power flow and some other constraints. Uncertainties related to load, price of energy, renewables, and consumption of storage are also considered when dealing with the problem. Moreover, the proposed scheme is structured in the form of nonconvex MINLP. Noting the nonconvex nonlinear form of the problem, an adaptive robust optimization is adopted, and uncertainties are modeled. The paper uses a hybrid of Ant Lion Optimizer and Krill Herd Optimization algorithm so that optimal solutions are found. A new optimization approach is incorporated in [29] to model the flexible-reliable operation (FRO) of energy hubs (EHs) in electricity, natural gas, and district heating networks. In an attempt to provide a flexible EH, the paper also adopts storage and incentive-based demand response program (IDRP). Table 1 reports the previous pieces of research in this realm.Table 1
A report on the background research.
RefHub sources and storageHub participation in DA market ofProblem modelRESCHPESSEVBoilerTSSEnergyReserve[14]NoYesYesNoNoYesNoNoNLP[15]YesYesYesNoNoNoNoNoNLP[16]YesYesYesYesNoYesNoNoNLP[17]YesYesNoNoNoNoNoNoNLP[18]YesYesYesNoNoNoNoNoNLP[19]YesYesNoNoNoNoNoNoMINLP[20]NoYesNoYesNoNoNoNoLP[21]NoYesYesNoNoNoNoNoNLP[22]NoYesNoNoYesNoNoNoNLP[23]YesYesYesNoNoNoYesNoLP[24]YesYesYesYesYesNoYesNoMILP[25]YesNoYesNoNoNoNoNoNLP[26]YesNoYesYesNoNoNoNoNLP[27]YesNoYesYesNoNoNoNoNLP[28]YesYesYesYesYesYesNoNoNLP[29]YesYesYesNoYesYesNoNoNLPProposed modelYesYesYesYesYesYesYesYesMILPLP: linear programming; MILP: mixed-integer LP.As per Table1 and the research background, major contributions to the field of the hub operation are summarized as follows:(i)
In a few pieces of research such as [16, 20, 24, 26–28], models of EVs in the energy hub have been considered. Note, however, that EVs have the potential to reduce pollution, and they are also able to enhance system flexibility with RES accessibility, so it is anticipated that their presence in the hub will be certain in the coming years. Additionally, CHPs with RESs reduce the flexibility of the hub heating section. As the thermal power output of CHP relies on its active power, it cannot be controlled. Therefore, it is necessary to use flexibility sources such as boiler and TSS to improve the flexibility of the hub, although this has also been considered in less research.(ii)
Hubs can improve various indices of energy networks thanks to including various sources and active loads. Thus, they can increase their financial benefit by engaging in ancillary services markets like the reserve market, which has rarely been studied. Also, in most research such as [14–22, 25–29], the operation model of hubs has been considered, while considering the model of the energy market can obtain a suitable financial benefit for it, and this is also the case that has rarely been addressed [23, 24].(iii)
The operation model of the energy hub in different energy networks such as NLP and MINLP is nonconvex. Therefore, its solvers generally obtain the optimal solution in a high computational time. Yet, in operation problems, especially real-time problems, solution time is of great importance. Also, their solvers generally cannot find the optimal single solution. To compensate for this, some research works such as [20, 23, 24] have used linear modeling for the main problem, the solvers of which can calculate the unique response in the shortest possible time.To address such challenges, the current study models the participation of the networked hubs based on the CEM strategy in DA energy and reserve markets (refer to Figure1). The hub contains RES, CHP, ESS, EVs parking lot, boiler, and TSS to supply the demand in EGDH networks. The proposed strategy maximizes the total profit of the hub in DA energy and up and down reserve markets, where the energy market includes EGDH energies, but the reserve market considers electrical and heating reserve services. Also, the problem is constrained by OPF equations in the proposed energy networks and hubs and reserve constraints. The proposed strategy contains an MINLP model that finds the local optimal solution at a high computational time because of nonconvex and nonlinear equations. Hence, the paper employs a mixed-integer linear programming (MILP) problem using conventional linearization methods to address this issue. Some innovations and shares of the present study include the following:(i)
Adopting a CEM strategy to simultaneously manage different power sources and storage in the hub framework to obtain high benefits to hubs in different networks and markets.(ii)
Modeling the participation of networked hubs in day-ahead energy and reserve markets to reach higher profit for hubs in the mentioned markets.(iii)
Presenting a linear form of the proposed nonconvex and nonlinear original model to achieve the optimal solution with a low error during a shorter time.Figure 1
The suggested framework for the networked hubs partnership in the DA markets.Moreover, the main objectives of the paper include the following:(i)
Using an appropriate aggregator and coordinator framework for resources and storage to enhance energy efficiency.(ii)
Gaining suitable financial benefits for resources and storage in energy markets and auxiliary services.(iii)
Improving the technical status of energy networks in the presence of energy hubs.(iv)
Solving the operation problem in a low computational time.Sections2 and 3 express original nonlinear and linear models of management of networked hubs in DA energy and reserve markets. Sections 4 and 5 present simulation results and underlines the benefits of the suggested strategy, respectively.


## 2. Proposed Original Model

The optimization model of the energy hubs scheduling in EGDH networks based on their participation in DA energy and reserve markets is described in this part of the paper. Therefore, it maximizes the total energy hubs profit in DA energy and up reserve and down reserve markets subjected to different energy network models and energy hub and reserve formulation.


### 2.1. Objective Function

The objective function equation (1) refers to maximizing hubs’ benefit in DA energy and reserve markets. The first term of the equation expresses hubs profit in the DA energy market obtained by selling active and reactive power, P and Q, gas power, G, and heating power, H [24]. Additionally, revenue or profit of hubs in up and down reserve markets arising from selling active and heating power is expressed by the second and third terms of (1) [30]. In this paper, it is assumed that hubs do not include local gas power sources and receive gas power from gas stations; hence, these systems cannot participate in reserve markets. Moreover, equation (1) refers to hub revenue/cost if active, gas, and heating powers in different proposed markets and is positive/negative. Finally, the reactive power exchanged between the hub and the electrical network is based on enhancing network indices, including voltage profile; therefore, revenue is obtained for hubs in both capacitive and inductive operations [31].(1)maxProfit=∑t∈ST∑m∈HλtEPm,t+KQQm,t+λtGGm,t+λtHHm,t+∑t∈ST∑m∈HρtUERUm,tE+ρtUHRUm,tH+∑t∈ST∑m∈HρtDERDm,tE+ρtDHRDm,tH.


### 2.2. Network Model

This section contains the power flow model, (2)–(10), constraints of technical indices, (11)–(19), in electrical, gas, and heating systems, and reserve requirement formulation, (20)–(23), in electrical and district heating grids. The power flow formulation in the electrical network is expressed in constraints (2)–(6) that are referred, respectively, to nodal active and reactive power balance, active and reactive power flowing from the distribution line, and the voltage angle in the slack bus [32–34]. Moreover, the gas/heating power flow model is presented in equations (7)–(8)/(9)–(10) so that constraint (7)/(9) explains the nodal gas/heating power balance model and equations (8)/(10) refer to the gas/heating flowing from a pipeline at hour t [16]. In addition, constraints of technical indices in different proposed energy networks are expressed in (11)–(19), which are presented, respectively, to limitations of the electrical bus voltage, the capacity of distribution line and station, gas pressure in different nodes, the capacity of pipeline and station in the gas network, heating nodes temperature, and heating pipeline and station capacity [16, 20]. It is considered, in this paper, that different power stations are connected to the slack bus of EGDH networks; hence, variables PS, QS, GS, and HS are zero in the rest of the busses/nodes.Up and down reserve conditions in the electrical/district heating network during timet are extracted using limits (20)-(21)/(22)-(23), respectively [30]. According to these constraints, reserve powers will be provided by energy hubs in this paper.(2)Pe,tS+∑m∈HAe,mEPm,t−∑j∈ENBe,jEPe,j,tL=Pe,tD∀e,t,(3)Qe,tS+∑m∈HAe,mEQm,t−∑j∈ENBe,jEQe,j,tL=Qe,tD∀e,t,(4)Pe,j,tL=GLe,jVe,t2−Ve,tVj,tGLe,jcosθe,t−θj,t+BLe,jsinθe,t−θj,t∀e,j,t,(5)Qe,j,tL=−BLe,jVe,t2+Ve,tVj,tBLe,jcosθe,t−θj,t−GLe,jsinθe,t−θj,t∀e,j,t,(6)θe,t=0∀e=Slack bus,t,(7)Gg,tS+∑m∈HAg,mGGm,t−∑j∈GNBg,jGGg,j,tL=Gg,tD∀g,t,(8)Gg,j,tL=κg,jsignπg,t,πj,tsignπg,t,πj,tπg,t2−πj,t2g,j,t,(9)Hh,tS+∑m∈HAh,mHHm,t−∑j∈HNBh,jHHh,j,tL=Hh,tD∀h,t,(10)Hh,j,tL=ch,jm˙h,jTh,t−Tj,t∀h,j,t,(11)V¯e≤Ve,t≤V¯e∀e,t,(12)Pe,j,tL2+Qe,j,tL2≤S¯e,jL∀e,j,t,(13)Pe,tS2+Qe,tS2≤S¯eS∀e=Slack bus,t,(14)π¯g≤πg,t≤π¯g∀g,t,(15)−G¯g,jL≤Gg,j,tL≤G¯g,jL∀g,j,t,(16)−G¯gS≤Gg,tS≤G¯gS∀g=Slack bus,t,(17)T¯h≤Th,t≤T¯h∀h,t,(18)−H¯h,jL≤Hh,j,tL≤H¯h,jL∀h,j,t,(19)−H¯hS≤Hh,tS≤H¯hS∀h=Slack bus,t,(20)∑m∈HRUm,tE≥DUtE∀t,(21)∑m∈HRDm,tE≥DDtE∀t,(22)∑m∈HRUm,tH≥DUtH∀t,(23)∑m∈HRDm,tH≥DDtH∀t.


### 2.3. Hub Formulation

The hub, in this paper, contains electrical and thermal sources and storage, such as RES, CHP, boiler,EVs parking lot, ESS, and TSS, where it connects to the different energy networks to inject and absorb the different powers. Energy hubs generally employ RESs to reduce pollution [23]. They also use CHPs to increase their energy efficiency [23]. Because of uncertain forecasting of active generation power of RESs, the hub flexibility decreases. As a solution, the use of flexibility sources like ESS and EV along with RES is suggested [24]. In addition, the operation of CHP in the presence of RES will reduce the flexibility of the heating section of the hub [14, 22]. Following this, the boiler and TSS are used as sources of flexibility in the thermal part of the hub [14, 22]. As another point, this paper assumes that CHP, ESS, and EVs are also capable of controlling reactive power, but RESs produce only active power according to IEEE 1547 [24]. In the following, the hub model is described by constraints (24)–(55) so that active, reactive, gas, and heating power balance in the hub are presented in (24)–(27), respectively. The CHP model is formulated in (28)–(31) that are referred, respectively, to the calculation of the CHP output heating and input gas power, (28) and (29), electrical and heating capacity limits of CHP, (30) and (31) [20]. The boiler is modeled based on (32) and (33) [24], where its input gas power is calculated using constraint (32), and its output heating power is limited as (33). Also, TSS equations are explained in (34)–(38) [35] that are related to the constraint of the thermal stored energy in TSS at hour t, charging and discharging heating power limits, the constraint of the initial thermal energy, and thermal stored energy limit in TSS, respectively. The ESS model is presented in (39)–(44) [36], where the format of equations (39)–(43) is the same as constraints (34)–(38) in the TSS model, but active power and electrical stored energy variables are used in this formulation. Also, constraint (44) introduces the charger capacity limit of the ESS. Moreover, the EVs parking lot formulation is modeled in (45)–(50) [31] so that the stored electrical energy of EVs available in the parking lot at the time t is expressed as in (45). The charging and discharging rate limits of EVs are based on (46) and (47), respectively. Moreover, the initial and final energy levels of EVs in the parking lot at arrival and departure hours are presented in (48) and (49), respectively, and the charger capacity of all EVs in the parking lot at the hour t is according to (50). It is noted that, in constraints (35) and (36), the binary variable of ts helps prevent charging-discharging operation at the same time in the TSS. This issue can be implemented in ESS and EVs parking lots using binary variables of es and ev, respectively, according to constraints (40), (41), (46), and (47). Besides, the terms PEV,C, PEV,D, and S¯EV in constraints of EVs parking lots are equal to ∑v=1NpCRv, ∑v=1NpDRv, and ∑v=1NpCCv, respectively. CR and DR refer to the charging and discharging rates of the EV battery, CC represents the EV’s charger capacity, and Np denotes the EVs number in the parking lot at hour t. Also, EA and ED in equations (48) and (49) are calculated by ∑v=1NASOCvBCv and∑v=1NDBCv, respectively, where SOC and BC are related to the state of charge (SOC) of the EV’s battery and its capacity and NA and ND are, respectively, the numbers of EVs at arrival and departure times [31].The active/heating power of the hub in up and down reserve systems is calculated by (51)-(52)/(53)-(54), respectively. These variables are always positive as per (55) [30]. Note that the maximum/minimum capacity of its local sources and storage determines the maximum/minimum size of the hub to inject/absorb power (51)–(54).(24)Pm,t=Pm,tC+Pm,tR+Pm,tES,D−Pm,tES,C+Pm,tEV,D−Pm,tEV,C−Pm,tD∀m,t,(25)Qm,t=Qm,tC+Qm,tES+Qm,tEV−Qm,tD∀m,t,(26)Gm,t=−Gm,tC−Gm,tB−Gm,tD∀m,t,(27)Hm,t=Hm,tC+Hm,tB+Hm,tTS,D−Hm,tTS,C−Hm,tD∀m,t,(28)Hm,tC=Pm,tC1−ηmT−ηmLηmT⋅ηmH∀m,t,(29)Gm,tC=Pm,tCηmT∀m,t,(30)Pm,tC2+Qm,tC2≤S¯mC∀m,t,(31)0≤Hm,tC≤H¯mC∀m,t,(32)Gm,tB=Hm,tBηmB∀m,t,(33)0≤Hm,tB≤H¯mB∀m,t,(34)Em,t+1TS=Em,tTS+ηTS,CHm,tTS,C−1ηTS,DHm,tTS,D∀m,t,(35)0≤Hm,tTS,C≤H¯mTS,Ctsm,t∀m,t,(36)0≤Hm,tTS,D≤H¯mTS,D1−tsm,t∀m,t,(37)Em,tTS=EmTS,I∀m,t=1,(38)E¯mTS≤Em,tTS≤E¯mTS∀m,t,(39)Em,t+1ES=Em,tES+ηES,CPm,tES,C−1ηES,DPm,tES,D∀m,t,(40)0≤Pm,tES,C≤P¯mES,Cesm,t∀m,t,(41)0≤Pm,tES,D≤P¯mES,D1−esm,t∀m,t,(42)Em,tES=EmES,I∀m,t=1,(43)E¯mES≤Em,tES≤E¯mES∀m,t,(44)Pm,tES,D−Pm,tES,C2+Qm,tES2≤S¯mES∀m,t,(45)Em,t+1EV=Em,tEV+ηEV,CPm,tEV,C−1ηEV,DPm,tEV,D∀m,t,(46)0≤Pm,tEV,C≤P¯m,tEV,Cevm,t∀m,t,(47)0≤Pm,tEV,D≤P¯m,tEV,D1−evm,t∀m,t,(48)Em,tEV=Em,tA∀m,t=Arrival time,(49)Em,tEV=Em,tD∀m,t=Departure time,(50)Pm,tEV,D−Pm,tEV,C2+Qm,tEV2≤S¯m,tEV∀m,t,(51)Pm,t+RUm,tE≤S¯mC+P¯mES,D+P¯m,tEV,D+Pm,tR∀m,t,(52)Pm,t−RDm,tE≥−P¯mES,C−P¯m,tEV,C∀m,t,(53)Hm,t+RUm,tH≤H¯mC+H¯mB+H¯mTS,D∀m,t,(54)Hm,t−RDm,tH≥−H¯mTS,C∀m,t,(55)RUm,tE,RDm,tE,RUm,tH,RDm,tH≥0∀m,t.


## 2.1. Objective Function

The objective function equation (1) refers to maximizing hubs’ benefit in DA energy and reserve markets. The first term of the equation expresses hubs profit in the DA energy market obtained by selling active and reactive power, P and Q, gas power, G, and heating power, H [24]. Additionally, revenue or profit of hubs in up and down reserve markets arising from selling active and heating power is expressed by the second and third terms of (1) [30]. In this paper, it is assumed that hubs do not include local gas power sources and receive gas power from gas stations; hence, these systems cannot participate in reserve markets. Moreover, equation (1) refers to hub revenue/cost if active, gas, and heating powers in different proposed markets and is positive/negative. Finally, the reactive power exchanged between the hub and the electrical network is based on enhancing network indices, including voltage profile; therefore, revenue is obtained for hubs in both capacitive and inductive operations [31].(1)maxProfit=∑t∈ST∑m∈HλtEPm,t+KQQm,t+λtGGm,t+λtHHm,t+∑t∈ST∑m∈HρtUERUm,tE+ρtUHRUm,tH+∑t∈ST∑m∈HρtDERDm,tE+ρtDHRDm,tH.


## 2.2. Network Model

This section contains the power flow model, (2)–(10), constraints of technical indices, (11)–(19), in electrical, gas, and heating systems, and reserve requirement formulation, (20)–(23), in electrical and district heating grids. The power flow formulation in the electrical network is expressed in constraints (2)–(6) that are referred, respectively, to nodal active and reactive power balance, active and reactive power flowing from the distribution line, and the voltage angle in the slack bus [32–34]. Moreover, the gas/heating power flow model is presented in equations (7)–(8)/(9)–(10) so that constraint (7)/(9) explains the nodal gas/heating power balance model and equations (8)/(10) refer to the gas/heating flowing from a pipeline at hour t [16]. In addition, constraints of technical indices in different proposed energy networks are expressed in (11)–(19), which are presented, respectively, to limitations of the electrical bus voltage, the capacity of distribution line and station, gas pressure in different nodes, the capacity of pipeline and station in the gas network, heating nodes temperature, and heating pipeline and station capacity [16, 20]. It is considered, in this paper, that different power stations are connected to the slack bus of EGDH networks; hence, variables PS, QS, GS, and HS are zero in the rest of the busses/nodes.Up and down reserve conditions in the electrical/district heating network during timet are extracted using limits (20)-(21)/(22)-(23), respectively [30]. According to these constraints, reserve powers will be provided by energy hubs in this paper.(2)Pe,tS+∑m∈HAe,mEPm,t−∑j∈ENBe,jEPe,j,tL=Pe,tD∀e,t,(3)Qe,tS+∑m∈HAe,mEQm,t−∑j∈ENBe,jEQe,j,tL=Qe,tD∀e,t,(4)Pe,j,tL=GLe,jVe,t2−Ve,tVj,tGLe,jcosθe,t−θj,t+BLe,jsinθe,t−θj,t∀e,j,t,(5)Qe,j,tL=−BLe,jVe,t2+Ve,tVj,tBLe,jcosθe,t−θj,t−GLe,jsinθe,t−θj,t∀e,j,t,(6)θe,t=0∀e=Slack bus,t,(7)Gg,tS+∑m∈HAg,mGGm,t−∑j∈GNBg,jGGg,j,tL=Gg,tD∀g,t,(8)Gg,j,tL=κg,jsignπg,t,πj,tsignπg,t,πj,tπg,t2−πj,t2g,j,t,(9)Hh,tS+∑m∈HAh,mHHm,t−∑j∈HNBh,jHHh,j,tL=Hh,tD∀h,t,(10)Hh,j,tL=ch,jm˙h,jTh,t−Tj,t∀h,j,t,(11)V¯e≤Ve,t≤V¯e∀e,t,(12)Pe,j,tL2+Qe,j,tL2≤S¯e,jL∀e,j,t,(13)Pe,tS2+Qe,tS2≤S¯eS∀e=Slack bus,t,(14)π¯g≤πg,t≤π¯g∀g,t,(15)−G¯g,jL≤Gg,j,tL≤G¯g,jL∀g,j,t,(16)−G¯gS≤Gg,tS≤G¯gS∀g=Slack bus,t,(17)T¯h≤Th,t≤T¯h∀h,t,(18)−H¯h,jL≤Hh,j,tL≤H¯h,jL∀h,j,t,(19)−H¯hS≤Hh,tS≤H¯hS∀h=Slack bus,t,(20)∑m∈HRUm,tE≥DUtE∀t,(21)∑m∈HRDm,tE≥DDtE∀t,(22)∑m∈HRUm,tH≥DUtH∀t,(23)∑m∈HRDm,tH≥DDtH∀t.


## 2.3. Hub Formulation

The hub, in this paper, contains electrical and thermal sources and storage, such as RES, CHP, boiler,EVs parking lot, ESS, and TSS, where it connects to the different energy networks to inject and absorb the different powers. Energy hubs generally employ RESs to reduce pollution [23]. They also use CHPs to increase their energy efficiency [23]. Because of uncertain forecasting of active generation power of RESs, the hub flexibility decreases. As a solution, the use of flexibility sources like ESS and EV along with RES is suggested [24]. In addition, the operation of CHP in the presence of RES will reduce the flexibility of the heating section of the hub [14, 22]. Following this, the boiler and TSS are used as sources of flexibility in the thermal part of the hub [14, 22]. As another point, this paper assumes that CHP, ESS, and EVs are also capable of controlling reactive power, but RESs produce only active power according to IEEE 1547 [24]. In the following, the hub model is described by constraints (24)–(55) so that active, reactive, gas, and heating power balance in the hub are presented in (24)–(27), respectively. The CHP model is formulated in (28)–(31) that are referred, respectively, to the calculation of the CHP output heating and input gas power, (28) and (29), electrical and heating capacity limits of CHP, (30) and (31) [20]. The boiler is modeled based on (32) and (33) [24], where its input gas power is calculated using constraint (32), and its output heating power is limited as (33). Also, TSS equations are explained in (34)–(38) [35] that are related to the constraint of the thermal stored energy in TSS at hour t, charging and discharging heating power limits, the constraint of the initial thermal energy, and thermal stored energy limit in TSS, respectively. The ESS model is presented in (39)–(44) [36], where the format of equations (39)–(43) is the same as constraints (34)–(38) in the TSS model, but active power and electrical stored energy variables are used in this formulation. Also, constraint (44) introduces the charger capacity limit of the ESS. Moreover, the EVs parking lot formulation is modeled in (45)–(50) [31] so that the stored electrical energy of EVs available in the parking lot at the time t is expressed as in (45). The charging and discharging rate limits of EVs are based on (46) and (47), respectively. Moreover, the initial and final energy levels of EVs in the parking lot at arrival and departure hours are presented in (48) and (49), respectively, and the charger capacity of all EVs in the parking lot at the hour t is according to (50). It is noted that, in constraints (35) and (36), the binary variable of ts helps prevent charging-discharging operation at the same time in the TSS. This issue can be implemented in ESS and EVs parking lots using binary variables of es and ev, respectively, according to constraints (40), (41), (46), and (47). Besides, the terms PEV,C, PEV,D, and S¯EV in constraints of EVs parking lots are equal to ∑v=1NpCRv, ∑v=1NpDRv, and ∑v=1NpCCv, respectively. CR and DR refer to the charging and discharging rates of the EV battery, CC represents the EV’s charger capacity, and Np denotes the EVs number in the parking lot at hour t. Also, EA and ED in equations (48) and (49) are calculated by ∑v=1NASOCvBCv and∑v=1NDBCv, respectively, where SOC and BC are related to the state of charge (SOC) of the EV’s battery and its capacity and NA and ND are, respectively, the numbers of EVs at arrival and departure times [31].The active/heating power of the hub in up and down reserve systems is calculated by (51)-(52)/(53)-(54), respectively. These variables are always positive as per (55) [30]. Note that the maximum/minimum capacity of its local sources and storage determines the maximum/minimum size of the hub to inject/absorb power (51)–(54).(24)Pm,t=Pm,tC+Pm,tR+Pm,tES,D−Pm,tES,C+Pm,tEV,D−Pm,tEV,C−Pm,tD∀m,t,(25)Qm,t=Qm,tC+Qm,tES+Qm,tEV−Qm,tD∀m,t,(26)Gm,t=−Gm,tC−Gm,tB−Gm,tD∀m,t,(27)Hm,t=Hm,tC+Hm,tB+Hm,tTS,D−Hm,tTS,C−Hm,tD∀m,t,(28)Hm,tC=Pm,tC1−ηmT−ηmLηmT⋅ηmH∀m,t,(29)Gm,tC=Pm,tCηmT∀m,t,(30)Pm,tC2+Qm,tC2≤S¯mC∀m,t,(31)0≤Hm,tC≤H¯mC∀m,t,(32)Gm,tB=Hm,tBηmB∀m,t,(33)0≤Hm,tB≤H¯mB∀m,t,(34)Em,t+1TS=Em,tTS+ηTS,CHm,tTS,C−1ηTS,DHm,tTS,D∀m,t,(35)0≤Hm,tTS,C≤H¯mTS,Ctsm,t∀m,t,(36)0≤Hm,tTS,D≤H¯mTS,D1−tsm,t∀m,t,(37)Em,tTS=EmTS,I∀m,t=1,(38)E¯mTS≤Em,tTS≤E¯mTS∀m,t,(39)Em,t+1ES=Em,tES+ηES,CPm,tES,C−1ηES,DPm,tES,D∀m,t,(40)0≤Pm,tES,C≤P¯mES,Cesm,t∀m,t,(41)0≤Pm,tES,D≤P¯mES,D1−esm,t∀m,t,(42)Em,tES=EmES,I∀m,t=1,(43)E¯mES≤Em,tES≤E¯mES∀m,t,(44)Pm,tES,D−Pm,tES,C2+Qm,tES2≤S¯mES∀m,t,(45)Em,t+1EV=Em,tEV+ηEV,CPm,tEV,C−1ηEV,DPm,tEV,D∀m,t,(46)0≤Pm,tEV,C≤P¯m,tEV,Cevm,t∀m,t,(47)0≤Pm,tEV,D≤P¯m,tEV,D1−evm,t∀m,t,(48)Em,tEV=Em,tA∀m,t=Arrival time,(49)Em,tEV=Em,tD∀m,t=Departure time,(50)Pm,tEV,D−Pm,tEV,C2+Qm,tEV2≤S¯m,tEV∀m,t,(51)Pm,t+RUm,tE≤S¯mC+P¯mES,D+P¯m,tEV,D+Pm,tR∀m,t,(52)Pm,t−RDm,tE≥−P¯mES,C−P¯m,tEV,C∀m,t,(53)Hm,t+RUm,tH≤H¯mC+H¯mB+H¯mTS,D∀m,t,(54)Hm,t−RDm,tH≥−H¯mTS,C∀m,t,(55)RUm,tE,RDm,tE,RUm,tH,RDm,tH≥0∀m,t.


## 3. MILP Model

The proposed original model of the hub scheduling, given in (1)–(55), includes a nonconvex MINLP formulation because of nonlinear constraints (4), (5), (8), (12), (13), (30), (44), and (50), the first part of (1), and nonconvex equations (4), (5), and (8) [20]. Hence, the present study converts this model to an MILP form using the classic linearization method to find the unique optimal solution with a low error at a low computational time [34]. The original method generally obtains local optimum due to nonconvex equations in a longer period arising from its solving approach that is based on the iteration-based numerical method like Newton–Raphson strategy [32]. Also, the solvers of the MINLP model fail to obtain a unique optimal solution. Finally, the suggested MILP model conforms to the following linearization methods.


### 3.1. Absolute Term

In (1), the term |Q| can be written as Qp–Qn, where Qp and Qn ≥ 0 are positive and negative components of Q. Therefore, Q = Qp + Qn.


### 3.2. Electrical Power Flow Equations

Constraints (4) and (5) are nonconvex and nonlinear [34]. The voltage angle difference between two busses across a distribution line is less than 6° or 0.105 rad [34]. Hence, terms cosθe,t−θj,t and sinθe,t−θj,t can be approximated to 1 and θe,t−θj,t, respectively. Moreover, the voltage magnitude varies between 0.9 and 1p.u. and is close to p.u. [34]. As a result, the voltage magnitude, V, can be expressed as 1 + ΔV, where ΔV<< 1 represents the voltage deviation [34]. Finally, the linear equations of (56) and (57) can be obtained for nonlinear constraints (4) and (5) considering the proposed assumption and removing terms ΔV2, (θe,t−θj,t)2, and (θe,t−θj,t) ΔV because they have small values.(56)Pe,j,tL=GLe,jΔVe,t−ΔVj,t−BLe,jθe,t−θj,t∀e,j,t,(57)Qe,j,tL=−BLe,jΔVe,t−ΔVj,t−GLe,jθe,t−θj,t∀e,j,t.


### 3.3. Gas Power Flow Equations

By rewriting (8) as constraint (58), where δ is an auxiliary variable, δg,j,t2=signπg,t,πj,tπg,t2−πj,t2 based on constraint (8). Note that δ and π can be expressed as equations (59) and (60), respectively, according to the conventional piecewise linearization method [34]. Parameters Δδ and Δπ are new auxiliary variables for linearizing (8), and also Δδ can be written as constraint (61) to obtain a linear format for δg,j,t2=signπg,t,πj,tπg,t2−πj,t2, where mπ and mδ denote line slopes for variables π and δ, respectively.(58)Gg,j,tL=κg,jsignπg,t,πj,tδg,j,t∀g,j,t,(59)δg,j,t=∑l∈LΔδg,j,t,l∀g,j,t,(60)πg,t=π¯g+∑l∈LΔπg,t,l∀g,t,(61)Δδg,j,t,l=signπg,πjmlπmlδΔπg,t,l−Δπj,t,l∀g,j,t,l.


### 3.4. Circular Plane Constraints

Capacity limits of the distribution line and station, CHP, ESS charger, andEVs chargers in the parking lot are in the form of a circular plane P2+Q2≤S, based on equations (12), (13), (30), (44), and (50). A polygon plane can be used to approximate the circular plane and obtain a linearization form for these constraints [31]. Each side of the polygon is a straight line, the equations of which are calculated using the tangents of the circle at different points as P×cosk×Δα+Q×sink×Δα=S in the PQ plane and the radius of S [31]. Hence, the linear form of the plane P2+Q2≤S can be written asP×cosk×Δα+Q×sink×Δα≤S, where Δα refers to the angle deviation and is equal to 2π/n, and k selects from a set of K = {1, 2,…,n}, where n is the total number of straight lines in a polygon. Therefore, the linear format of constraints (12), (13), (30), (44), and (50) can be written as (62)–(66), respectively.(62)Pe,j,tLcosk×Δα+Qe,j,tLsink×Δα≤S¯e,jL∀e,j,t,k,(63)Pe,tScosk×Δα+Qe,tSsink×Δα≤S¯eS∀e=Slack bus,t,k,(64)Pm,tCcosk×Δα+Qm,tCsink×Δα≤S¯mC∀m,t,k,(65)Pm,tES,D−Pm,tES,Ccosk×Δα+Qm,tESsink×Δα≤S¯mES∀m,t,k,(66)Pm,tEV,D−Pm,tEV,Ccosk×Δα+Qm,tEVsink×Δα≤S¯m,tEV∀m,t,k.The suggested MILP model based on the mentioned linearization methods will be(67)max Profit=∑t∈ST∑m∈HλtEPm,t+KQQm,tp−Qm,tn+λtGGm,t+λtHHm,t+∑t∈ST∑m∈HρtUERUm,tE+ρtUHRUm,tH+∑t∈ST∑m∈HρtDERDm,tE+ρtDHRDm,tH,subject to(68)Qm,t=Qm,tp+Qm,tn∀m,t,Qm,tp,Qm,tn≥0,(69)V¯e−1≤ΔVe,t≤V¯e−1∀e,t.Equations of (2), (3), (6), (7), (9), (10), (14)–(29), (31)–(43), (45)–(49), (51)–(66).The linear format of the objective function in (1) is expressed in (29), where it is followed by the first linearization method. Also, constraint (68) obtains according to this linearization approach. According to the second linearization approach, the voltage deviation is adopted and not voltage magnitude; thus, the voltage deviation limit (69) will be used from (64) in the proposed MILP model. In the end, the other required linear equations of the problem are presented in (70).


## 3.1. Absolute Term

In (1), the term |Q| can be written as Qp–Qn, where Qp and Qn ≥ 0 are positive and negative components of Q. Therefore, Q = Qp + Qn.


## 3.2. Electrical Power Flow Equations

Constraints (4) and (5) are nonconvex and nonlinear [34]. The voltage angle difference between two busses across a distribution line is less than 6° or 0.105 rad [34]. Hence, terms cosθe,t−θj,t and sinθe,t−θj,t can be approximated to 1 and θe,t−θj,t, respectively. Moreover, the voltage magnitude varies between 0.9 and 1p.u. and is close to p.u. [34]. As a result, the voltage magnitude, V, can be expressed as 1 + ΔV, where ΔV<< 1 represents the voltage deviation [34]. Finally, the linear equations of (56) and (57) can be obtained for nonlinear constraints (4) and (5) considering the proposed assumption and removing terms ΔV2, (θe,t−θj,t)2, and (θe,t−θj,t) ΔV because they have small values.(56)Pe,j,tL=GLe,jΔVe,t−ΔVj,t−BLe,jθe,t−θj,t∀e,j,t,(57)Qe,j,tL=−BLe,jΔVe,t−ΔVj,t−GLe,jθe,t−θj,t∀e,j,t.


## 3.3. Gas Power Flow Equations

By rewriting (8) as constraint (58), where δ is an auxiliary variable, δg,j,t2=signπg,t,πj,tπg,t2−πj,t2 based on constraint (8). Note that δ and π can be expressed as equations (59) and (60), respectively, according to the conventional piecewise linearization method [34]. Parameters Δδ and Δπ are new auxiliary variables for linearizing (8), and also Δδ can be written as constraint (61) to obtain a linear format for δg,j,t2=signπg,t,πj,tπg,t2−πj,t2, where mπ and mδ denote line slopes for variables π and δ, respectively.(58)Gg,j,tL=κg,jsignπg,t,πj,tδg,j,t∀g,j,t,(59)δg,j,t=∑l∈LΔδg,j,t,l∀g,j,t,(60)πg,t=π¯g+∑l∈LΔπg,t,l∀g,t,(61)Δδg,j,t,l=signπg,πjmlπmlδΔπg,t,l−Δπj,t,l∀g,j,t,l.


## 3.4. Circular Plane Constraints

Capacity limits of the distribution line and station, CHP, ESS charger, andEVs chargers in the parking lot are in the form of a circular plane P2+Q2≤S, based on equations (12), (13), (30), (44), and (50). A polygon plane can be used to approximate the circular plane and obtain a linearization form for these constraints [31]. Each side of the polygon is a straight line, the equations of which are calculated using the tangents of the circle at different points as P×cosk×Δα+Q×sink×Δα=S in the PQ plane and the radius of S [31]. Hence, the linear form of the plane P2+Q2≤S can be written asP×cosk×Δα+Q×sink×Δα≤S, where Δα refers to the angle deviation and is equal to 2π/n, and k selects from a set of K = {1, 2,…,n}, where n is the total number of straight lines in a polygon. Therefore, the linear format of constraints (12), (13), (30), (44), and (50) can be written as (62)–(66), respectively.(62)Pe,j,tLcosk×Δα+Qe,j,tLsink×Δα≤S¯e,jL∀e,j,t,k,(63)Pe,tScosk×Δα+Qe,tSsink×Δα≤S¯eS∀e=Slack bus,t,k,(64)Pm,tCcosk×Δα+Qm,tCsink×Δα≤S¯mC∀m,t,k,(65)Pm,tES,D−Pm,tES,Ccosk×Δα+Qm,tESsink×Δα≤S¯mES∀m,t,k,(66)Pm,tEV,D−Pm,tEV,Ccosk×Δα+Qm,tEVsink×Δα≤S¯m,tEV∀m,t,k.The suggested MILP model based on the mentioned linearization methods will be(67)max Profit=∑t∈ST∑m∈HλtEPm,t+KQQm,tp−Qm,tn+λtGGm,t+λtHHm,t+∑t∈ST∑m∈HρtUERUm,tE+ρtUHRUm,tH+∑t∈ST∑m∈HρtDERDm,tE+ρtDHRDm,tH,subject to(68)Qm,t=Qm,tp+Qm,tn∀m,t,Qm,tp,Qm,tn≥0,(69)V¯e−1≤ΔVe,t≤V¯e−1∀e,t.Equations of (2), (3), (6), (7), (9), (10), (14)–(29), (31)–(43), (45)–(49), (51)–(66).The linear format of the objective function in (1) is expressed in (29), where it is followed by the first linearization method. Also, constraint (68) obtains according to this linearization approach. According to the second linearization approach, the voltage deviation is adopted and not voltage magnitude; thus, the voltage deviation limit (69) will be used from (64) in the proposed MILP model. In the end, the other required linear equations of the problem are presented in (70).


## 4. Numerical Results

### 4.1. Case Study

A case study system as shown in Figure2 is used to test the suggested approach [20]. The system includes a 9-bus electrical network, a 4-node natural gas system, and a 7-node district heating grid. The base power is 1 MW, and the base values of voltage, pressures, and temperature are 1 kV, 10 Bar, and 100 C, respectively. The allowed range of [0.9, 1.1] p.u. is considered for voltage [37–41], pressures [24], and temperature [4], and capacities of distribution, gas, and heating stations are 7, 11, and 3 p.u., respectively [24]. Characteristics of distribution lines and pipelines in different networks are presented in [20]. Also, the data of peak electrical and thermal loads are extracted from [20], and the gas load is set to zero. The daily load profile is given by multiplying the peak load and load factor daily curve, as shown in Figure 3(a) [23]. The daily energy price curve for EGDH networks is depicted in Figure 3(b) [23]. Moreover, the daily curve for up and down reserve limits in DA electrical [42] and heating markets is illustrated in Figures 3(c) and 3(d), respectively. Coefficient KQ is set at 0.08, and the daily up and down reserve price curve is considered the same as energy price in different networks. In the proposed test system, there are 7 hubs with locations as shown in Figure 2. Hubs 1, 2, 3, and 5 contain electrical sources and storage, that is, RES, ESS, and the parking lot of EVs. Hub 4 uses CHP, boiler, and TSS, and all proposed sources and storage are located in hubs 6 and 7. The capacity of CHP in both electrical and thermal sections is assumed to be 1 p.u.; that is, the terms ηH, ηL, and ηT are 0.4, 0.08, and 0.4, respectively [24]. The output capacity of the boiler with an efficiency of 0.8 is 1 p.u. Concerning the ESS, the minimum and maximum storable energy are 0.5 p.u. and 2 p.u., respectively. Its initial energy, charging rate, discharging rate, and charger capacity are assumed to be 0.5, 0.3, 0.3, and 0.5 p.u., respectively. Characteristics of TSS are similar to ESS, except for the charging and discharging efficiencies that are 0.8. In this paper, two types of RES are used that are wind and PV systems. The capacity of both of these systems at each hub is 0.25 and 0.2 p.u. Their daily power generation profile is given by multiplying the capacity of RES and the daily curve of RES, as illustrated in Figure 3(e) [43]. Furthermore, hubs 1–3 and 5–7 consist of 60 EVs. Figure 3(f) demonstrates the daily curve of EVs’ penetration [31]. Characteristics of the rest of EVs, including SOC, charging and discharging rates of the EV, and battery capacity, are available in [24, 31].Figure 2
Illustration of the case study system [20].Figure 3
Daily curves of (a) load factor, (b) energy price, (c) electrical reserve requirement [42], (d) heating reserve requirements, (e) RES power rate, and (f) EVs penetration rate.
(a)(b)(c)(d)(e)(f)It should be noted that the formulation presented in Sections2 and 3 for the proposed scheme does not have any restrictions on different related data. Therefore, it can be implemented on the real data of networks, resources, and storage devices. The data of resources and storage reported in this paper are in accordance with the real data. In addition, the proposed scheme is presented with the assumption that there is an automation platform in networks.


### 4.2. Results

The GAMS software is employed to code the suggested approach; thus, BONMIN and CPLEX solvers help find solutions to the suggested MINLP and MILP models, respectively [44].


#### 4.2.1. Capabilities of the Proposed MILP Model

Findings are tabulated in Table2, where the numbers of linearization segments for the classic piecewise linearization method and circular plane in the proposed MILP method are 5 and 45, respectively. Among the BARON, BONMIN, and KNITRO solvers for the MINLP problem, BONMIN provides the optimal solution (maximum profit) in the shortest time (637.44 s) and the number of convergence iterations (741) compared to the other two solvers. Also, for the proposed problem, the mentioned algorithms have different solutions, so KNITRO cannot reach the convergence point. Therefore, they do not have a unique solution. However, in the proposed MILP model, its solvers, namely, CONOPT, CBC, and CPLEX, have one optimal point. Only each of them obtains the optimal solution in different convergence times and the number of iterations so that the CPLEX algorithm has the minimum computational time.Table 2
Comparison results between the proposed MINLP and MILP methods.
Convergence results in the different solversModelSolverConvergence iterationCalculation time (sec)Profit ($)MINLPBARON923845.171871.7BONMIN741637.441958.1KNITRO——InfeasibleMILPCONOPT2119.0822037.4CBC1467.1022037.4CPLEX975.8662037.4Results of different models at electrical peak load hour (20 : 00) using suitable solversParameterMINLPMILPCalculation error (%)SolverBONMINCPLEX—Distribution station active power (p.u.)3.4323.3542.33Distribution station reactive power (p.u.)2.332−2.3872.27Gas station power (p.u.)3.5303.5000.87Heating station power (p.u.)2.1842.1840Average voltage magnitude (p.u.)1.04941.0550.53Average voltage angle (rad)−0.24407−0.2430.44Average pressures (p.u.)0.974220.9750.08Average temperature (p.u.)0.9340.9340Computation time637.445.866—In the following, the computational error of different problem variables based on the results of BONMIN in the MINLP model and CPLEX in the MILP model is presented in Table2. As the table reports, the computational error in the MILP model compared to the MINLP model is about 2.3% ((3.432–3.354)/3.432), 0.9%, 0.5%, and 0.1%, respectively, for electrical power, gas power, voltage, and pressure. Nevertheless, the computational error is zero for variables in the district heating network because optimal power flow equations for this system are the same as those of the proposed MILP and MINLP approaches given in Sections 2 and 3. Moreover, the proposed MILP method finds the optimal solution within 5.866 s, while the original MINLP method can achieve it at 637.44 s, which is very large compared to that of the MILP model. Therefore, the presented MILP approach can find the optimal solution with a low computational error at a low calculation time in comparison with the originally proposed model. Hence, an MILP model can replace the MINLP model. Finally, the discussion provided in this section validates the third research gap and the third contribution given in Section 1.


#### 4.2.2. Capabilities of the Coordinated Energy Management

In this section, the capabilities of uncoordinated and coordinated energy management (UEM and CEM) methods are investigated, and their numerical results are presented in Table3. Note that, in the UEM, the energy management strategy is implemented on different networks, including CHP, boiler, TSS, RES, ESS, and EVs parking lots separately. As an instance, in the UEM approach, variables PC, QC, −GC, and HC replace P, Q, G, and H in the first term of (1) and constraints (2)-(3), (7), (9) and (51)–(54) if this approach is utilized to manage different networks including only CHP. Therefore, the proposed method adds only equations of CHP, that is, (28)–(31), to optimal power flow equations (2)–(23); thus, the hubs load is added to network loads, and matrix A presents the incidence matrix of busses and CHP. This rule is the same for energy management of UEM-based networks when all networks use one of the power sources or storage devices. However, all sources and storage are managed simultaneously in the hub if the CEM method, (56)–(70), is used for all networks. Finally, the profit values of hubs in DA energy and reserve markets are expressed in Table 3. Accordingly, the total profit of hub devices in the DA energy market is equal to $303.2 (121.1 + 57.4 + 36.7 + 74.9 + 62.1 − 49) based on the UEM method, while in the CEM approach, it is $353.7, which has increased by about 16.7% ((353.7 − 303.2)/303.2) compared to the UEM method. Moreover, the hub’s profit in DA up and down reserve markets based on the CEM method increased by about 27.7% and 15.1%, respectively, compared to the UEM method, and the results are listed in rows 9–11. Moreover, the income of source and ESS based on UEM and CEM strategies is higher in the reserve market than in the energy market. Because their income in the energy market is commensurate with the technical constraints of different networks, (11)–(19), their income in the reserve market is based on (2)−(23), independent of the technical constraints of the networks. Also, due to the high energy demand for EVs, they are unable to gain income in the energy market. Therefore, their profit in this market is negative, which means cost. Therefore, it is suitable to implement the CEM approach on different sources and storage based on a hub framework to achieve a higher profit in different markets for hub devices, where this section confirms the advantages of the first innovation mentioned in Section 1.Table 3
Hubs profit for DA markets when adopting various energy management methods.
Energy management methodDA marketEnergyUp reserveDown reserveUEMCHP121.1211.531.7Boiler57.4120.727.5TSS36.784.123.9RES74.9142.421.4ESS62.1289.137.8EVs parking lot-49301.745.2Sum303.21149.5187.5CEM353.71467.9215.8Gain (%)16.727.715.1


#### 4.2.3. Hub Power and Profit in Different DA Markets

Two cases are investigated in this section: the proposed strategy without (Case I) and with (Case II) considering system operation limits (11)−(19). Figure 4 presents the daily power and profit curves of all hubs in DA energy and up and down reserve markets for Case I. As Figure 4(a) shows, the negative active power value of hubs from hours 1 : 00 to 8 : 00 expresses that ESSs and EVs are charged at this period corresponding to low electrical energy price and load hours based on the data of [24] to obtain low energy cost according to the objective function (1). Also, all electrical sources and storage of hubs can supply high active power to the electrical network during 9 : 00–00 : 00 to provide a higher profit for hubs in the DA energy market based on (1) because this power is a positive value at these hours. Note that hubs inject high active power into the network at the period of 17 : 00–00 : 00 compared to the period of 9 : 00–16 : 00 because ESSs and EVs are discharged, and CHPs and RESs inject active power into the electrical network at hours 17 : 00–24 : 00 due to high electrical energy price at these hours based on the data of [24], while CHPs and RESs are used from 9 : 00 to 16 : 00. In addition, hubs inject high reactive power into the electrical system from 1 : 00 to 8 : 00 and 17 : 00 to 22 : 00 as the total number of reactive power sources such as ESSs, EVs parking lot, and CHPs during these periods is more than that at other hours. Moreover, hubs inject high/low heating power into the district heating network at periods of 8 : 00–15 : 00 and 23 : 00–00 : 00/1 : 00–7 : 00 and 16 : 00–22 : 00. Consequently, the TSS is charged from 1 : 00 to 7 : 00 corresponding to low heating energy price hours according to the data of [24]. Thus, TSS, boiler, and CHP inject high heating power into the heating system from 8 : 00 to 15 : 00 and 23 : 00 to 00 : 00 to obtain high profit for hubs in the DA energy market. Nonetheless, at hours 16 : 00–22 : 00, the TSS is switched off, and CHP and boiler generate low heating power in comparison to other hours due to low/high heating/gas energy price during this period based on the data in [24]. Finally, the daily gas power curve is fixed at roughly −4.2 p.u. since the gas demand depends on the operation of the boiler and CHP.Figure 4
Daily hubs power curve in (a) energy, (b) electrical reserve, and (c) heating reserve markets and (d) the daily hubs profit curve for Case I.
(a)(b)(c)(d)Daily up and down reserve power curves of hubs in electrical and heating networks are plotted in Figures4(b) and 4(c), respectively. Accordingly, the changing form of the hubs’ reserve power is the same as the daily reserve requirement curve, but the reserve power of hubs is higher than the networks’ reserve requirements due to obtaining the highest profit for hubs in DA up and down reserve markets. In addition, the daily hubs profit curve in DA energy and up and down reserve markets is shown in Figure 4(d). Note that the hub’s profit in DA energy and up and down reserve markets are based on the first, second, and third parts of (1), respectively, so that it depends on energy, reserve power, and energy price. Accordingly, the hubs profit in the DA energy market is negative at the period of 1 : 00–8 : 00 due to low active, gas, and heating power in this period based on Figure 4(a), while it is positive and high at other hours due to high heating and active and reactive power injection into electrical and heating networks at these hours. Up/down reserve profit of hubs depends on electrical and heating energy prices and up/down reserve power according to (1). Therefore, the hourly changing form of hubs profit in mentioned markets is based on Figure 4(d).The daily power and profit curve of all hubs in DA energy and up and down reserve markets for Case II, that is, considering system operation limits, is plotted in Figure5. By comparing the results of Cases I and II in Figures 4 and 5, it is seen that the daily power curve of hubs in different networks and markets in Case II shifted down with respect to Case I due to considering limits on network indices. For example, in Case II, hubs cannot absorb/inject active power from/into the electrical network as in Case I from 1 : 00 to 8 : 00/9 : 00 to 24 : 00 due to limits on bus voltages and distribution lines capacity. Therefore, the daily profit curve of hubs in different DA markets in Case II is shifted down as it depends on the hubs’ power and price of different markets according to (1). In addition, the economic results of hubs for Cases I and II are expressed in Table 4. Referring to this table, the hub’s profit in the DA energy market for Case II is reduced by about 25.5% ((474.7 − 353.7)/474.7) with respect to Case I. Also, this parameter in DA up and down reserve markets for Case II is less than that in Case I by about 7.3% and 8%, respectively. As a result, the total profit of hubs considering system operation limits is reduced by roughly 11.1% in comparison to Case I.Figure 5
Daily hubs power curves in (a) energy, (b) electrical reserve, and (c) heating reserve markets and (d) the daily hubs profit curve in Case II.
(a)(b)(c)(d)Table 4
Hubs profit in different DA markets for Cases I and II.
CaseIIIHubs profit in the DA energy market ($)474.7353.7Hubs profit in the DA up reserve market ($)1582.81467.9Hubs profit in the DA down reserve market ($)234.6215.8Total hubs profit2292.12037.4


#### 4.2.4. Capabilities of Hubs to Improve Network Indices

In this section, different cases are investigated to achieve the benefits of the proposed CEM method for networked hubs according to DA energy and reserve markets, where details of each case are presented in Table5. Also, the results of these cases are provided in Table 5, in which it is observed that Cases VI and XI can obtain the minimum electrical energy loss (EEL) and a lower maximum voltage deviation (MVD) than other cases. Note that all controllable electrical sources and storage are connected to different energy networks in these cases; hence, they can control their output power based on the minimum EEL and MVD. Moreover, EEL and MVD values are high in Case II due to the high injection power of RESs into the electrical network. Hence, the proposed scheme in Case XI can improve EEL and MVD by about 83.3% ((5.4102 − 2.9511)/5.4102) and 49.55% ((0.113 − 0.057)/0.113), respectively, with respect to Case II. Moreover, minimum values of heating energy loss (HEL) and maximum temperature deviation (MTD) are obtained in Cases X and XI because all controllable heating sources and storage are present in these cases, and hence they can improve indices of HEL and MTV in comparison with other Cases. The high value of these parameters is obtained in Case I (power flow analysis); thus, the proposed strategy in Case XI is able to improve HEL and MTD by about 37.7% and 38.7%, respectively, compared with Case I.Table 5
Technical indices for different cases.
CaseEEL (p.u.)GEL (p.u.)HEL (p.u.)I (power flow analysis)5.121102.1156II (power flow analysis considering RESs)5.410202.1156III (the hub with RES andEVs)4.397102.1156IV (the hub with RES and ESS)4.231202.1156V (the hub with RES and CHP)4.01261.11831.5832VI (the hub with RES, CHP, ESS, andEVs)2.95111.11831.5832VII (the hub with CHP)3.78051.18181.5312VIII (the hub with CHP and boiler)3.78051.98311.3703IX (the hub with CHP and TSS)3.78051.18181.4708X (the hub with CHP, boiler, and TSS)3.78051.98311.3176XI (proposed method)2.95111.98311.3176CaseMVD (p.u.)MPD (p.u.)MTD (p.u.)I (power flow analysis)0.11300.111II (power flow analysis with considering RESs)0.12800.111III (the hub with RES andEVs)0.07300.111IV (the hub with RES and ESS)0.06900.111V (the hub with RES and CHP)0.0640.0110.082VI (the hub with RES, CHP, ESS, andEVs)0.0570.0110.082VII (the hub with CHP)0.0610.0140.080VIII (the hub with CHP and boiler)0.0610.0380.070IX (the hub with CHP and TSS)0.0610.0140.074X (the hub with CHP, boiler, and TSS)0.0610.0380.068XI (proposed method)0.0570.0380.068Finally, the gas energy loss (GEL) and the maximum pressure deviation (MPD) are minimized in Cases I–IV; there are no gas customers such as CHP and boiler in these cases, and the gas load is zero; thus, GEL and MPD are zero. As a result, by increasing the GEL (MPD) by 1.9831 (0.038) p.u. in Case XI with respect to Case I, the proposed scheme can significantly reduce EEL, HEL, MVD, and MTD.


## 4.1. Case Study

A case study system as shown in Figure2 is used to test the suggested approach [20]. The system includes a 9-bus electrical network, a 4-node natural gas system, and a 7-node district heating grid. The base power is 1 MW, and the base values of voltage, pressures, and temperature are 1 kV, 10 Bar, and 100 C, respectively. The allowed range of [0.9, 1.1] p.u. is considered for voltage [37–41], pressures [24], and temperature [4], and capacities of distribution, gas, and heating stations are 7, 11, and 3 p.u., respectively [24]. Characteristics of distribution lines and pipelines in different networks are presented in [20]. Also, the data of peak electrical and thermal loads are extracted from [20], and the gas load is set to zero. The daily load profile is given by multiplying the peak load and load factor daily curve, as shown in Figure 3(a) [23]. The daily energy price curve for EGDH networks is depicted in Figure 3(b) [23]. Moreover, the daily curve for up and down reserve limits in DA electrical [42] and heating markets is illustrated in Figures 3(c) and 3(d), respectively. Coefficient KQ is set at 0.08, and the daily up and down reserve price curve is considered the same as energy price in different networks. In the proposed test system, there are 7 hubs with locations as shown in Figure 2. Hubs 1, 2, 3, and 5 contain electrical sources and storage, that is, RES, ESS, and the parking lot of EVs. Hub 4 uses CHP, boiler, and TSS, and all proposed sources and storage are located in hubs 6 and 7. The capacity of CHP in both electrical and thermal sections is assumed to be 1 p.u.; that is, the terms ηH, ηL, and ηT are 0.4, 0.08, and 0.4, respectively [24]. The output capacity of the boiler with an efficiency of 0.8 is 1 p.u. Concerning the ESS, the minimum and maximum storable energy are 0.5 p.u. and 2 p.u., respectively. Its initial energy, charging rate, discharging rate, and charger capacity are assumed to be 0.5, 0.3, 0.3, and 0.5 p.u., respectively. Characteristics of TSS are similar to ESS, except for the charging and discharging efficiencies that are 0.8. In this paper, two types of RES are used that are wind and PV systems. The capacity of both of these systems at each hub is 0.25 and 0.2 p.u. Their daily power generation profile is given by multiplying the capacity of RES and the daily curve of RES, as illustrated in Figure 3(e) [43]. Furthermore, hubs 1–3 and 5–7 consist of 60 EVs. Figure 3(f) demonstrates the daily curve of EVs’ penetration [31]. Characteristics of the rest of EVs, including SOC, charging and discharging rates of the EV, and battery capacity, are available in [24, 31].Figure 2
Illustration of the case study system [20].Figure 3
Daily curves of (a) load factor, (b) energy price, (c) electrical reserve requirement [42], (d) heating reserve requirements, (e) RES power rate, and (f) EVs penetration rate.
(a)(b)(c)(d)(e)(f)It should be noted that the formulation presented in Sections2 and 3 for the proposed scheme does not have any restrictions on different related data. Therefore, it can be implemented on the real data of networks, resources, and storage devices. The data of resources and storage reported in this paper are in accordance with the real data. In addition, the proposed scheme is presented with the assumption that there is an automation platform in networks.


## 4.2. Results

The GAMS software is employed to code the suggested approach; thus, BONMIN and CPLEX solvers help find solutions to the suggested MINLP and MILP models, respectively [44].


### 4.2.1. Capabilities of the Proposed MILP Model

Findings are tabulated in Table2, where the numbers of linearization segments for the classic piecewise linearization method and circular plane in the proposed MILP method are 5 and 45, respectively. Among the BARON, BONMIN, and KNITRO solvers for the MINLP problem, BONMIN provides the optimal solution (maximum profit) in the shortest time (637.44 s) and the number of convergence iterations (741) compared to the other two solvers. Also, for the proposed problem, the mentioned algorithms have different solutions, so KNITRO cannot reach the convergence point. Therefore, they do not have a unique solution. However, in the proposed MILP model, its solvers, namely, CONOPT, CBC, and CPLEX, have one optimal point. Only each of them obtains the optimal solution in different convergence times and the number of iterations so that the CPLEX algorithm has the minimum computational time.Table 2
Comparison results between the proposed MINLP and MILP methods.
Convergence results in the different solversModelSolverConvergence iterationCalculation time (sec)Profit ($)MINLPBARON923845.171871.7BONMIN741637.441958.1KNITRO——InfeasibleMILPCONOPT2119.0822037.4CBC1467.1022037.4CPLEX975.8662037.4Results of different models at electrical peak load hour (20 : 00) using suitable solversParameterMINLPMILPCalculation error (%)SolverBONMINCPLEX—Distribution station active power (p.u.)3.4323.3542.33Distribution station reactive power (p.u.)2.332−2.3872.27Gas station power (p.u.)3.5303.5000.87Heating station power (p.u.)2.1842.1840Average voltage magnitude (p.u.)1.04941.0550.53Average voltage angle (rad)−0.24407−0.2430.44Average pressures (p.u.)0.974220.9750.08Average temperature (p.u.)0.9340.9340Computation time637.445.866—In the following, the computational error of different problem variables based on the results of BONMIN in the MINLP model and CPLEX in the MILP model is presented in Table2. As the table reports, the computational error in the MILP model compared to the MINLP model is about 2.3% ((3.432–3.354)/3.432), 0.9%, 0.5%, and 0.1%, respectively, for electrical power, gas power, voltage, and pressure. Nevertheless, the computational error is zero for variables in the district heating network because optimal power flow equations for this system are the same as those of the proposed MILP and MINLP approaches given in Sections 2 and 3. Moreover, the proposed MILP method finds the optimal solution within 5.866 s, while the original MINLP method can achieve it at 637.44 s, which is very large compared to that of the MILP model. Therefore, the presented MILP approach can find the optimal solution with a low computational error at a low calculation time in comparison with the originally proposed model. Hence, an MILP model can replace the MINLP model. Finally, the discussion provided in this section validates the third research gap and the third contribution given in Section 1.


### 4.2.2. Capabilities of the Coordinated Energy Management

In this section, the capabilities of uncoordinated and coordinated energy management (UEM and CEM) methods are investigated, and their numerical results are presented in Table3. Note that, in the UEM, the energy management strategy is implemented on different networks, including CHP, boiler, TSS, RES, ESS, and EVs parking lots separately. As an instance, in the UEM approach, variables PC, QC, −GC, and HC replace P, Q, G, and H in the first term of (1) and constraints (2)-(3), (7), (9) and (51)–(54) if this approach is utilized to manage different networks including only CHP. Therefore, the proposed method adds only equations of CHP, that is, (28)–(31), to optimal power flow equations (2)–(23); thus, the hubs load is added to network loads, and matrix A presents the incidence matrix of busses and CHP. This rule is the same for energy management of UEM-based networks when all networks use one of the power sources or storage devices. However, all sources and storage are managed simultaneously in the hub if the CEM method, (56)–(70), is used for all networks. Finally, the profit values of hubs in DA energy and reserve markets are expressed in Table 3. Accordingly, the total profit of hub devices in the DA energy market is equal to $303.2 (121.1 + 57.4 + 36.7 + 74.9 + 62.1 − 49) based on the UEM method, while in the CEM approach, it is $353.7, which has increased by about 16.7% ((353.7 − 303.2)/303.2) compared to the UEM method. Moreover, the hub’s profit in DA up and down reserve markets based on the CEM method increased by about 27.7% and 15.1%, respectively, compared to the UEM method, and the results are listed in rows 9–11. Moreover, the income of source and ESS based on UEM and CEM strategies is higher in the reserve market than in the energy market. Because their income in the energy market is commensurate with the technical constraints of different networks, (11)–(19), their income in the reserve market is based on (2)−(23), independent of the technical constraints of the networks. Also, due to the high energy demand for EVs, they are unable to gain income in the energy market. Therefore, their profit in this market is negative, which means cost. Therefore, it is suitable to implement the CEM approach on different sources and storage based on a hub framework to achieve a higher profit in different markets for hub devices, where this section confirms the advantages of the first innovation mentioned in Section 1.Table 3
Hubs profit for DA markets when adopting various energy management methods.
Energy management methodDA marketEnergyUp reserveDown reserveUEMCHP121.1211.531.7Boiler57.4120.727.5TSS36.784.123.9RES74.9142.421.4ESS62.1289.137.8EVs parking lot-49301.745.2Sum303.21149.5187.5CEM353.71467.9215.8Gain (%)16.727.715.1


### 4.2.3. Hub Power and Profit in Different DA Markets

Two cases are investigated in this section: the proposed strategy without (Case I) and with (Case II) considering system operation limits (11)−(19). Figure 4 presents the daily power and profit curves of all hubs in DA energy and up and down reserve markets for Case I. As Figure 4(a) shows, the negative active power value of hubs from hours 1 : 00 to 8 : 00 expresses that ESSs and EVs are charged at this period corresponding to low electrical energy price and load hours based on the data of [24] to obtain low energy cost according to the objective function (1). Also, all electrical sources and storage of hubs can supply high active power to the electrical network during 9 : 00–00 : 00 to provide a higher profit for hubs in the DA energy market based on (1) because this power is a positive value at these hours. Note that hubs inject high active power into the network at the period of 17 : 00–00 : 00 compared to the period of 9 : 00–16 : 00 because ESSs and EVs are discharged, and CHPs and RESs inject active power into the electrical network at hours 17 : 00–24 : 00 due to high electrical energy price at these hours based on the data of [24], while CHPs and RESs are used from 9 : 00 to 16 : 00. In addition, hubs inject high reactive power into the electrical system from 1 : 00 to 8 : 00 and 17 : 00 to 22 : 00 as the total number of reactive power sources such as ESSs, EVs parking lot, and CHPs during these periods is more than that at other hours. Moreover, hubs inject high/low heating power into the district heating network at periods of 8 : 00–15 : 00 and 23 : 00–00 : 00/1 : 00–7 : 00 and 16 : 00–22 : 00. Consequently, the TSS is charged from 1 : 00 to 7 : 00 corresponding to low heating energy price hours according to the data of [24]. Thus, TSS, boiler, and CHP inject high heating power into the heating system from 8 : 00 to 15 : 00 and 23 : 00 to 00 : 00 to obtain high profit for hubs in the DA energy market. Nonetheless, at hours 16 : 00–22 : 00, the TSS is switched off, and CHP and boiler generate low heating power in comparison to other hours due to low/high heating/gas energy price during this period based on the data in [24]. Finally, the daily gas power curve is fixed at roughly −4.2 p.u. since the gas demand depends on the operation of the boiler and CHP.Figure 4
Daily hubs power curve in (a) energy, (b) electrical reserve, and (c) heating reserve markets and (d) the daily hubs profit curve for Case I.
(a)(b)(c)(d)Daily up and down reserve power curves of hubs in electrical and heating networks are plotted in Figures4(b) and 4(c), respectively. Accordingly, the changing form of the hubs’ reserve power is the same as the daily reserve requirement curve, but the reserve power of hubs is higher than the networks’ reserve requirements due to obtaining the highest profit for hubs in DA up and down reserve markets. In addition, the daily hubs profit curve in DA energy and up and down reserve markets is shown in Figure 4(d). Note that the hub’s profit in DA energy and up and down reserve markets are based on the first, second, and third parts of (1), respectively, so that it depends on energy, reserve power, and energy price. Accordingly, the hubs profit in the DA energy market is negative at the period of 1 : 00–8 : 00 due to low active, gas, and heating power in this period based on Figure 4(a), while it is positive and high at other hours due to high heating and active and reactive power injection into electrical and heating networks at these hours. Up/down reserve profit of hubs depends on electrical and heating energy prices and up/down reserve power according to (1). Therefore, the hourly changing form of hubs profit in mentioned markets is based on Figure 4(d).The daily power and profit curve of all hubs in DA energy and up and down reserve markets for Case II, that is, considering system operation limits, is plotted in Figure5. By comparing the results of Cases I and II in Figures 4 and 5, it is seen that the daily power curve of hubs in different networks and markets in Case II shifted down with respect to Case I due to considering limits on network indices. For example, in Case II, hubs cannot absorb/inject active power from/into the electrical network as in Case I from 1 : 00 to 8 : 00/9 : 00 to 24 : 00 due to limits on bus voltages and distribution lines capacity. Therefore, the daily profit curve of hubs in different DA markets in Case II is shifted down as it depends on the hubs’ power and price of different markets according to (1). In addition, the economic results of hubs for Cases I and II are expressed in Table 4. Referring to this table, the hub’s profit in the DA energy market for Case II is reduced by about 25.5% ((474.7 − 353.7)/474.7) with respect to Case I. Also, this parameter in DA up and down reserve markets for Case II is less than that in Case I by about 7.3% and 8%, respectively. As a result, the total profit of hubs considering system operation limits is reduced by roughly 11.1% in comparison to Case I.Figure 5
Daily hubs power curves in (a) energy, (b) electrical reserve, and (c) heating reserve markets and (d) the daily hubs profit curve in Case II.
(a)(b)(c)(d)Table 4
Hubs profit in different DA markets for Cases I and II.
CaseIIIHubs profit in the DA energy market ($)474.7353.7Hubs profit in the DA up reserve market ($)1582.81467.9Hubs profit in the DA down reserve market ($)234.6215.8Total hubs profit2292.12037.4


### 4.2.4. Capabilities of Hubs to Improve Network Indices

In this section, different cases are investigated to achieve the benefits of the proposed CEM method for networked hubs according to DA energy and reserve markets, where details of each case are presented in Table5. Also, the results of these cases are provided in Table 5, in which it is observed that Cases VI and XI can obtain the minimum electrical energy loss (EEL) and a lower maximum voltage deviation (MVD) than other cases. Note that all controllable electrical sources and storage are connected to different energy networks in these cases; hence, they can control their output power based on the minimum EEL and MVD. Moreover, EEL and MVD values are high in Case II due to the high injection power of RESs into the electrical network. Hence, the proposed scheme in Case XI can improve EEL and MVD by about 83.3% ((5.4102 − 2.9511)/5.4102) and 49.55% ((0.113 − 0.057)/0.113), respectively, with respect to Case II. Moreover, minimum values of heating energy loss (HEL) and maximum temperature deviation (MTD) are obtained in Cases X and XI because all controllable heating sources and storage are present in these cases, and hence they can improve indices of HEL and MTV in comparison with other Cases. The high value of these parameters is obtained in Case I (power flow analysis); thus, the proposed strategy in Case XI is able to improve HEL and MTD by about 37.7% and 38.7%, respectively, compared with Case I.Table 5
Technical indices for different cases.
CaseEEL (p.u.)GEL (p.u.)HEL (p.u.)I (power flow analysis)5.121102.1156II (power flow analysis considering RESs)5.410202.1156III (the hub with RES andEVs)4.397102.1156IV (the hub with RES and ESS)4.231202.1156V (the hub with RES and CHP)4.01261.11831.5832VI (the hub with RES, CHP, ESS, andEVs)2.95111.11831.5832VII (the hub with CHP)3.78051.18181.5312VIII (the hub with CHP and boiler)3.78051.98311.3703IX (the hub with CHP and TSS)3.78051.18181.4708X (the hub with CHP, boiler, and TSS)3.78051.98311.3176XI (proposed method)2.95111.98311.3176CaseMVD (p.u.)MPD (p.u.)MTD (p.u.)I (power flow analysis)0.11300.111II (power flow analysis with considering RESs)0.12800.111III (the hub with RES andEVs)0.07300.111IV (the hub with RES and ESS)0.06900.111V (the hub with RES and CHP)0.0640.0110.082VI (the hub with RES, CHP, ESS, andEVs)0.0570.0110.082VII (the hub with CHP)0.0610.0140.080VIII (the hub with CHP and boiler)0.0610.0380.070IX (the hub with CHP and TSS)0.0610.0140.074X (the hub with CHP, boiler, and TSS)0.0610.0380.068XI (proposed method)0.0570.0380.068Finally, the gas energy loss (GEL) and the maximum pressure deviation (MPD) are minimized in Cases I–IV; there are no gas customers such as CHP and boiler in these cases, and the gas load is zero; thus, GEL and MPD are zero. As a result, by increasing the GEL (MPD) by 1.9831 (0.038) p.u. in Case XI with respect to Case I, the proposed scheme can significantly reduce EEL, HEL, MVD, and MTD.


## 4.2.1. Capabilities of the Proposed MILP Model

Findings are tabulated in Table2, where the numbers of linearization segments for the classic piecewise linearization method and circular plane in the proposed MILP method are 5 and 45, respectively. Among the BARON, BONMIN, and KNITRO solvers for the MINLP problem, BONMIN provides the optimal solution (maximum profit) in the shortest time (637.44 s) and the number of convergence iterations (741) compared to the other two solvers. Also, for the proposed problem, the mentioned algorithms have different solutions, so KNITRO cannot reach the convergence point. Therefore, they do not have a unique solution. However, in the proposed MILP model, its solvers, namely, CONOPT, CBC, and CPLEX, have one optimal point. Only each of them obtains the optimal solution in different convergence times and the number of iterations so that the CPLEX algorithm has the minimum computational time.Table 2
Comparison results between the proposed MINLP and MILP methods.
Convergence results in the different solversModelSolverConvergence iterationCalculation time (sec)Profit ($)MINLPBARON923845.171871.7BONMIN741637.441958.1KNITRO——InfeasibleMILPCONOPT2119.0822037.4CBC1467.1022037.4CPLEX975.8662037.4Results of different models at electrical peak load hour (20 : 00) using suitable solversParameterMINLPMILPCalculation error (%)SolverBONMINCPLEX—Distribution station active power (p.u.)3.4323.3542.33Distribution station reactive power (p.u.)2.332−2.3872.27Gas station power (p.u.)3.5303.5000.87Heating station power (p.u.)2.1842.1840Average voltage magnitude (p.u.)1.04941.0550.53Average voltage angle (rad)−0.24407−0.2430.44Average pressures (p.u.)0.974220.9750.08Average temperature (p.u.)0.9340.9340Computation time637.445.866—In the following, the computational error of different problem variables based on the results of BONMIN in the MINLP model and CPLEX in the MILP model is presented in Table2. As the table reports, the computational error in the MILP model compared to the MINLP model is about 2.3% ((3.432–3.354)/3.432), 0.9%, 0.5%, and 0.1%, respectively, for electrical power, gas power, voltage, and pressure. Nevertheless, the computational error is zero for variables in the district heating network because optimal power flow equations for this system are the same as those of the proposed MILP and MINLP approaches given in Sections 2 and 3. Moreover, the proposed MILP method finds the optimal solution within 5.866 s, while the original MINLP method can achieve it at 637.44 s, which is very large compared to that of the MILP model. Therefore, the presented MILP approach can find the optimal solution with a low computational error at a low calculation time in comparison with the originally proposed model. Hence, an MILP model can replace the MINLP model. Finally, the discussion provided in this section validates the third research gap and the third contribution given in Section 1.


## 4.2.2. Capabilities of the Coordinated Energy Management

In this section, the capabilities of uncoordinated and coordinated energy management (UEM and CEM) methods are investigated, and their numerical results are presented in Table3. Note that, in the UEM, the energy management strategy is implemented on different networks, including CHP, boiler, TSS, RES, ESS, and EVs parking lots separately. As an instance, in the UEM approach, variables PC, QC, −GC, and HC replace P, Q, G, and H in the first term of (1) and constraints (2)-(3), (7), (9) and (51)–(54) if this approach is utilized to manage different networks including only CHP. Therefore, the proposed method adds only equations of CHP, that is, (28)–(31), to optimal power flow equations (2)–(23); thus, the hubs load is added to network loads, and matrix A presents the incidence matrix of busses and CHP. This rule is the same for energy management of UEM-based networks when all networks use one of the power sources or storage devices. However, all sources and storage are managed simultaneously in the hub if the CEM method, (56)–(70), is used for all networks. Finally, the profit values of hubs in DA energy and reserve markets are expressed in Table 3. Accordingly, the total profit of hub devices in the DA energy market is equal to $303.2 (121.1 + 57.4 + 36.7 + 74.9 + 62.1 − 49) based on the UEM method, while in the CEM approach, it is $353.7, which has increased by about 16.7% ((353.7 − 303.2)/303.2) compared to the UEM method. Moreover, the hub’s profit in DA up and down reserve markets based on the CEM method increased by about 27.7% and 15.1%, respectively, compared to the UEM method, and the results are listed in rows 9–11. Moreover, the income of source and ESS based on UEM and CEM strategies is higher in the reserve market than in the energy market. Because their income in the energy market is commensurate with the technical constraints of different networks, (11)–(19), their income in the reserve market is based on (2)−(23), independent of the technical constraints of the networks. Also, due to the high energy demand for EVs, they are unable to gain income in the energy market. Therefore, their profit in this market is negative, which means cost. Therefore, it is suitable to implement the CEM approach on different sources and storage based on a hub framework to achieve a higher profit in different markets for hub devices, where this section confirms the advantages of the first innovation mentioned in Section 1.Table 3
Hubs profit for DA markets when adopting various energy management methods.
Energy management methodDA marketEnergyUp reserveDown reserveUEMCHP121.1211.531.7Boiler57.4120.727.5TSS36.784.123.9RES74.9142.421.4ESS62.1289.137.8EVs parking lot-49301.745.2Sum303.21149.5187.5CEM353.71467.9215.8Gain (%)16.727.715.1


## 4.2.3. Hub Power and Profit in Different DA Markets

Two cases are investigated in this section: the proposed strategy without (Case I) and with (Case II) considering system operation limits (11)−(19). Figure 4 presents the daily power and profit curves of all hubs in DA energy and up and down reserve markets for Case I. As Figure 4(a) shows, the negative active power value of hubs from hours 1 : 00 to 8 : 00 expresses that ESSs and EVs are charged at this period corresponding to low electrical energy price and load hours based on the data of [24] to obtain low energy cost according to the objective function (1). Also, all electrical sources and storage of hubs can supply high active power to the electrical network during 9 : 00–00 : 00 to provide a higher profit for hubs in the DA energy market based on (1) because this power is a positive value at these hours. Note that hubs inject high active power into the network at the period of 17 : 00–00 : 00 compared to the period of 9 : 00–16 : 00 because ESSs and EVs are discharged, and CHPs and RESs inject active power into the electrical network at hours 17 : 00–24 : 00 due to high electrical energy price at these hours based on the data of [24], while CHPs and RESs are used from 9 : 00 to 16 : 00. In addition, hubs inject high reactive power into the electrical system from 1 : 00 to 8 : 00 and 17 : 00 to 22 : 00 as the total number of reactive power sources such as ESSs, EVs parking lot, and CHPs during these periods is more than that at other hours. Moreover, hubs inject high/low heating power into the district heating network at periods of 8 : 00–15 : 00 and 23 : 00–00 : 00/1 : 00–7 : 00 and 16 : 00–22 : 00. Consequently, the TSS is charged from 1 : 00 to 7 : 00 corresponding to low heating energy price hours according to the data of [24]. Thus, TSS, boiler, and CHP inject high heating power into the heating system from 8 : 00 to 15 : 00 and 23 : 00 to 00 : 00 to obtain high profit for hubs in the DA energy market. Nonetheless, at hours 16 : 00–22 : 00, the TSS is switched off, and CHP and boiler generate low heating power in comparison to other hours due to low/high heating/gas energy price during this period based on the data in [24]. Finally, the daily gas power curve is fixed at roughly −4.2 p.u. since the gas demand depends on the operation of the boiler and CHP.Figure 4
Daily hubs power curve in (a) energy, (b) electrical reserve, and (c) heating reserve markets and (d) the daily hubs profit curve for Case I.
(a)(b)(c)(d)Daily up and down reserve power curves of hubs in electrical and heating networks are plotted in Figures4(b) and 4(c), respectively. Accordingly, the changing form of the hubs’ reserve power is the same as the daily reserve requirement curve, but the reserve power of hubs is higher than the networks’ reserve requirements due to obtaining the highest profit for hubs in DA up and down reserve markets. In addition, the daily hubs profit curve in DA energy and up and down reserve markets is shown in Figure 4(d). Note that the hub’s profit in DA energy and up and down reserve markets are based on the first, second, and third parts of (1), respectively, so that it depends on energy, reserve power, and energy price. Accordingly, the hubs profit in the DA energy market is negative at the period of 1 : 00–8 : 00 due to low active, gas, and heating power in this period based on Figure 4(a), while it is positive and high at other hours due to high heating and active and reactive power injection into electrical and heating networks at these hours. Up/down reserve profit of hubs depends on electrical and heating energy prices and up/down reserve power according to (1). Therefore, the hourly changing form of hubs profit in mentioned markets is based on Figure 4(d).The daily power and profit curve of all hubs in DA energy and up and down reserve markets for Case II, that is, considering system operation limits, is plotted in Figure5. By comparing the results of Cases I and II in Figures 4 and 5, it is seen that the daily power curve of hubs in different networks and markets in Case II shifted down with respect to Case I due to considering limits on network indices. For example, in Case II, hubs cannot absorb/inject active power from/into the electrical network as in Case I from 1 : 00 to 8 : 00/9 : 00 to 24 : 00 due to limits on bus voltages and distribution lines capacity. Therefore, the daily profit curve of hubs in different DA markets in Case II is shifted down as it depends on the hubs’ power and price of different markets according to (1). In addition, the economic results of hubs for Cases I and II are expressed in Table 4. Referring to this table, the hub’s profit in the DA energy market for Case II is reduced by about 25.5% ((474.7 − 353.7)/474.7) with respect to Case I. Also, this parameter in DA up and down reserve markets for Case II is less than that in Case I by about 7.3% and 8%, respectively. As a result, the total profit of hubs considering system operation limits is reduced by roughly 11.1% in comparison to Case I.Figure 5
Daily hubs power curves in (a) energy, (b) electrical reserve, and (c) heating reserve markets and (d) the daily hubs profit curve in Case II.
(a)(b)(c)(d)Table 4
Hubs profit in different DA markets for Cases I and II.
CaseIIIHubs profit in the DA energy market ($)474.7353.7Hubs profit in the DA up reserve market ($)1582.81467.9Hubs profit in the DA down reserve market ($)234.6215.8Total hubs profit2292.12037.4


## 4.2.4. Capabilities of Hubs to Improve Network Indices

In this section, different cases are investigated to achieve the benefits of the proposed CEM method for networked hubs according to DA energy and reserve markets, where details of each case are presented in Table5. Also, the results of these cases are provided in Table 5, in which it is observed that Cases VI and XI can obtain the minimum electrical energy loss (EEL) and a lower maximum voltage deviation (MVD) than other cases. Note that all controllable electrical sources and storage are connected to different energy networks in these cases; hence, they can control their output power based on the minimum EEL and MVD. Moreover, EEL and MVD values are high in Case II due to the high injection power of RESs into the electrical network. Hence, the proposed scheme in Case XI can improve EEL and MVD by about 83.3% ((5.4102 − 2.9511)/5.4102) and 49.55% ((0.113 − 0.057)/0.113), respectively, with respect to Case II. Moreover, minimum values of heating energy loss (HEL) and maximum temperature deviation (MTD) are obtained in Cases X and XI because all controllable heating sources and storage are present in these cases, and hence they can improve indices of HEL and MTV in comparison with other Cases. The high value of these parameters is obtained in Case I (power flow analysis); thus, the proposed strategy in Case XI is able to improve HEL and MTD by about 37.7% and 38.7%, respectively, compared with Case I.Table 5
Technical indices for different cases.
CaseEEL (p.u.)GEL (p.u.)HEL (p.u.)I (power flow analysis)5.121102.1156II (power flow analysis considering RESs)5.410202.1156III (the hub with RES andEVs)4.397102.1156IV (the hub with RES and ESS)4.231202.1156V (the hub with RES and CHP)4.01261.11831.5832VI (the hub with RES, CHP, ESS, andEVs)2.95111.11831.5832VII (the hub with CHP)3.78051.18181.5312VIII (the hub with CHP and boiler)3.78051.98311.3703IX (the hub with CHP and TSS)3.78051.18181.4708X (the hub with CHP, boiler, and TSS)3.78051.98311.3176XI (proposed method)2.95111.98311.3176CaseMVD (p.u.)MPD (p.u.)MTD (p.u.)I (power flow analysis)0.11300.111II (power flow analysis with considering RESs)0.12800.111III (the hub with RES andEVs)0.07300.111IV (the hub with RES and ESS)0.06900.111V (the hub with RES and CHP)0.0640.0110.082VI (the hub with RES, CHP, ESS, andEVs)0.0570.0110.082VII (the hub with CHP)0.0610.0140.080VIII (the hub with CHP and boiler)0.0610.0380.070IX (the hub with CHP and TSS)0.0610.0140.074X (the hub with CHP, boiler, and TSS)0.0610.0380.068XI (proposed method)0.0570.0380.068Finally, the gas energy loss (GEL) and the maximum pressure deviation (MPD) are minimized in Cases I–IV; there are no gas customers such as CHP and boiler in these cases, and the gas load is zero; thus, GEL and MPD are zero. As a result, by increasing the GEL (MPD) by 1.9831 (0.038) p.u. in Case XI with respect to Case I, the proposed scheme can significantly reduce EEL, HEL, MVD, and MTD.


## 5. Conclusions

The paper presents the participation model of network hubs in DA energy and reserve markets, and the hub operator employs the CEM strategy to obtain optimal scheduling for hub devices. Therefore, this paper maximizes the total profit of hubs in DA energy and up and down reserve markets subjected to optimal power flow constraints of EGDH networks, reserve limits, and the hub model containing various power sources and storage devices. Following that, the proposed original model is converted to the MILP method. It is inferred from the findings that the suggested MILP approach finds the unique optimal solution with low computational error at a low computational time in comparison with the MINLP method. So, different solvers for the MILP provide an optimal unique solution to the proposed scheme, and this has not been obtained in the proposed MINLP. Additionally, the CPLEX algorithm for solving the MILP model has succeeded in extracting the optimal point within about 6 seconds; however, the best solver of the MINLP, that is, BONMIN, finds the solution within about 640 seconds. Moreover, the computational error of active and reactive power in the MILP model with respect to the original model of the problem, that is, MINLP, is about 2.5, and it is around 0.9%, 0.5%, and 0.1% for, respectively, gas power, voltage, and pressure. Considering the very low computational time of the MILP compared to that of the MINLP, this computational error can be neglected. In the following, using the CEM strategy, the hub operator can improve network indices, including energy loss and the maximum deviation of voltage, pressure, and temperature. Hence, the suggested scheme can decrease energy loss and voltage deviation in the electrical system by up to 83% and 45% by establishing low energy loss and pressure drop in the gas network caused by gas consumption of CHP and boiler. Furthermore, it reduces energy loss and temperature deviation in the heating network by 38% and 39%, respectively. Additionally, sources and ESSs in the hub framework and/or based on the CEM strategy can enhance their profit in energy and up and down reserve markets by 17%, 28%, and 15% compared to the UEM strategy, respectively. Also, the profit of hubs considering system operation limits is reduced by roughly 11% in comparison to case that does not consider system operation limits.



---

*Source: 6481531-2022-09-09.xml*