ParticipationofEnergyStorage-BasedFlexibleHubs inDay-Ahead Reserve Regulation and Energy Markets Based on a Coordinated Energy Management Strategy

Department of Electrical Engineering, Bandar Deylam Branch, Islamic Azad University, Bandar Deylam, Iran Department of Electrical Engineering, Islamic Azad University Science and Research Branch, Tehran, Iran Department of Electrical Engineering, Islamic Azad University, Borujerd Branch, Borujerd, Iran Department of Electrical Engineering, Faculty of Electrical and Computer Engineering, Technical and Vocational University (TVU), Tehran, Iran Department of Engineering, Semirom Branch, Islamic Azad University, Semirom, Iran Electrical Engineering Department, Faculty of Engineering, Shahed University, Tehran, Iran


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][2][3][4][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]. erefore, 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]. erefore, 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.
ere is huge research regarding hub energy management. e 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][17][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]. e operation model of hubs in different energy networks has widely been discussed in the literature. e 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. e 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]. e 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. e 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. ermal 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. e 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. e 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]. e 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. e 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.
As per Table 1 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][27][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 2 International Transactions on Electrical Energy Systems on its active power, it cannot be controlled. erefore, 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. us, 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][15][16][17][18][19][20][21][22][25][26][27][28][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) e operation model of the energy hub in different energy networks such as NLP and MINLP is nonconvex. erefore, 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 Figure 1). e hub contains RES, CHP, ESS, EVs parking lot, boiler, and TSS to supply the demand in EGDH networks. e 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. e 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.
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.
Sections 2 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.

Proposed Original Model
e 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. erefore, 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.

Objective Function.
e objective function equation (1) refers to maximizing hubs' benefit in DA energy and reserve markets. e 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)

Network Model.
is 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. e 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][33][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 P S , Q S , G S , and H S are zero in the rest of the busses/nodes.
Up and down reserve conditions in the electrical/district heating network during time t 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.

Hub Formulation.
e 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]. ey 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. e 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]. e 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. e 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). e 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. is issue can be implemented in ESS and EVs parking lots using binary variables of es and ev, respectively, according to constraints (40), (41) BC v , respectively, where SOC and BC are related to the state of charge (SOC) of the EV's battery and its capacity and N A and N D are, respectively, the numbers of EVs at arrival and departure times [31]. e active/heating power of the hub in up and down reserve systems is calculated by (51)-(52)/(53)-(54), respectively. ese 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).

MILP Model
e 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]. e 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.

Absolute Term.
In (1), the term |Q| can be written as Q p -Q n , where Q p and Q n ≥ 0 are positive and negative components of Q. erefore, Q � Q p + Q n . (4) and (5) are nonconvex and nonlinear [34]. e 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 1 p.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 ΔV 2 , (θ e,t −θ j,t ) 2 , and (θ e,t −θ j,t ) ΔV because they have small values.

Electrical Power Flow Equations. Constraints
Q L e,j,t � −BL e,j ΔV e,t − ΔV j,t − GL e,j θ e,t − θ j,t ∀e, j, t.

Circular Plane Constraints.
Capacity limits of the distribution line and station, CHP, ESS charger, and EVs chargers in the parking lot are in the form of a circular plane ������ � P 2 + Q 2 ≤ 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 × cos(k × Δα) + Q× sin(k × Δα) � S in the PQ plane and the radius of S [31]. Hence, the linear form of the plane ������ � P 2 + Q 2 ≤ S can be written asP × cos(k × Δα) + Q × sin(k × Δα) ≤ 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. erefore, the linear format of constraints (12), (13), (30), (44), and (50) can be written as (62)-(66), respectively. P L e,j,t cos(k × Δα) + Q L e,j,t sin(k × Δα)≤ S L e,j ∀e, j, t,k, Equations of (2) (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). Figure 2 is used to test the suggested approach [20]. e system includes a 9-bus electrical network, a 4-node natural gas system, and a 7node district heating grid. e base power is 1 MW, and the base values of voltage, pressures, and temperature are 1 kV, 10 Bar, and 100 C, respectively. e allowed range of [0.9, 1.1] p.u. is considered for voltage [37][38][39][40][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. e daily load profile is given by multiplying the peak load and load factor daily curve, as shown in Figure 3(a) [23]. e 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 K Q 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. e 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]. e 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., International Transactions on Electrical Energy Systems 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. e capacity of both of these systems at each hub is 0.25 and 0.2 p.u. eir 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].

Case Study. A case study system as shown in
It should be noted that the formulation presented in Sections 2 and 3 for the proposed scheme does not have any restrictions on different related data. erefore, it can be implemented on the real data of networks, resources, and storage devices. e 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.

Results.
e 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]. Table 2, 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. erefore, 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.

Capabilities of the Proposed MILP Model. Findings are tabulated in
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 Table 2. 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  Figure 2: Illustration of the case study system [20]. 8 International Transactions on Electrical Energy Systems 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. erefore, 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.

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 Table 3. 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 P C , Q C , −G C , and H C replace P, Q, G, and H in the first term of (1) and constraints (2) networks including only CHP. erefore, 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. is 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. 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. erefore, their profit in this market is negative, which means cost. erefore, 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.

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   International Transactions on Electrical Energy Systems 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]. us, 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. Daily up and down reserve power curves of hubs in electrical and heating networks are plotted in Figures 4(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). erefore, the hourly changing form of hubs profit in mentioned markets is based on Figure 4(d). e 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 Figure 5. 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. erefore, 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.

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 Table 5. 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. e 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.
Finally, the gas energy loss (GEL) and the maximum pressure deviation (MPD) are minimized in Cases I-IV;

Constants
ESS efficiency in charge mode η EV,DC : EVs efficiency in discharge mode η EV,C : EVs efficiency in charge mode η T : Turbine efficiency in the CHP η L : Loss efficiency in the CHP η H : ermal efficiency in the CHP η TS,D : TSS efficiency in discharge mode η TS,C : TSS efficiency in charge mode λ E : Energy price in the electricity market ($/MWh) λ G : Energy price in the gas market ($/MWh) λ H : Energy price in the heat market ($/MWh) ρ UE : Reserve price in the up mode for electricity market ($/MWh) ρ DE : Reserve price in the down mode for electricity market ($/MWh) ρ UH : Reserve price in the up mode for heat market ($/MWh) ρ DH : Reserve price in the down mode for heat market ($/MWh) κ: Gas pipeline constant (p.u.) π, π: Upper and lower value of pressure (p.u.).

Data Availability
No new data were created or analyzed in this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.