Coordinated energy management strategy in scheme of flexible grid-connected hubs participating in energy and reserve markets

Abstract

This paper uses the coordinated energy management strategy for different sources and storages in the framework of flexible grid-connected energy hubs that participate in the day-ahead (DA) energy and reserve markets. In the base problem, this method maximizes the difference between the expected revenue of hubs gained by selling energy and reserve power in the proposed markets and the expected cost of lost flexibility (COLF). Also, it is subject to linearized optimal power flow (LOPF) equations in the electricity, gas and district heating systems, as well as hub constraints including different sources, storages and reserve models. This problem contains uncertainties of load, market price, reserve requirement, renewable power and hub mobile storages parameters. Therefore, the hybrid stochastic/robust optimization (HSRO) is suitable to model these uncertain parameters and obtain robust capability for the hub to improve the system flexibility. Accordingly, the bounded uncertainty-based robust optimization (BURO) is used in this paper to model the uncertainty of hub mobile storages to achieve the hub robust potential in improving the system flexibility, and other uncertain parameters are modeled according to scenario-based stochastic programming (SBSP). Finally, the proposed strategy is implemented on a standard test system. The obtained numerical results confirm the capability of the suggested scheme in improving the economic status of sources and storages using the coordinated energy management strategy in the form of an energy hub, as well as enhancing economic conditions, operation, and flexibility of energy networks thanks to hubs for having access to optimal scheduling.

Keywords

Coordinated energy management cost of lost flexibility energy and reserve market flexible grid-connected energy hub hybrid stochastic/robust optimization

Nomenclature

Indices and sets:

h, g, e, k, l, m, t

Indices of heating, gas and electrical node, linearization parts of the circular plane, linearization parts of the conventional piecewise linearization method, hub, hour.

HN, GN, EN, K, L, H, ST

Set of heating, gas and electrical node, linearization parts of the circular plane, linearization parts of the conventional piecewise linearization method, hub, simulation time.

Variables:

E^ES, E^TS, E^EV

Stored energy in electrical storage system (ESS), thermal storage system (TSS) and electric vehicles (EVs) parking lot (p.u)

H^B, G^B

Boiler heating and gas power (p.u)

P, Q, G, H

Active, reactive, gas and heating power of hub (p.u)

P^C, Q^C, G^C, H^C

Active, reactive, gas and heating power of combined heat and power (CHP) system (p.u)

P^ES,C/P^ES,D

ESS charging/discharging power (p.u)

P^EV,C/ P^EV,D

EVs charging/discharging power in a parking lot (p.u)

P^L, Q^L, G^L, H^L

Active and reactive power of distribution line, gas and heating power of pipeline (p.u)

Profit

Hubs profit ($)

P^S, Q^S, G^S, H^S

Active, reactive, gas and heating power of station (p.u)

P^TS,C/P^TS,D

TSS charging/discharging power (p.u)

Q^ES, Q^EV

ESS reactive power, and EVs reactive power in a parking lot (p.u)

Q^p, Qⁿ

Positive and negative components of Q (p.u)

RU^E, RD^E

Up and down reserve power in the electrical network (p.u)

RU^H, RD^H

Up and down reserve power in district heating network (p.u)

st, ev, ts

ESS, TSS and EVs parking lot charging/discharging state, st, ev, ts ∈ {0, 1}

T, V, π

Temperature, voltage and pressure (p.u)

Voltage angle (rad)

Δ V, Δπ

Deviation of voltage and pressure (p.u)

δ, Δδ

Auxiliary variables (p.u)

Parameters:

A^E, B^E

Incidence matrix bus and hub, bus and distribution line.

A^G, B^G

Incidence matrix gas node and hub, gas node and pipeline.

A^H, B^H

Incidence matrix heating node and hub, heating node and heating pipeline.

BL, GL

Susceptance and conductance of a distribution line (p.u)

\dot{m}

Heating pipeline constants (p.u)

DU^E, DD^E

Up and down reserve requirements in the electrical network (p.u)

DU^H, DD^H

Up and down reserve requirements in the district heating network (p.u)

_EA,ED

Stored energy in EVs at arrival and departure hour (p.u)

{\underline{E}}^{ES}, {\bar{E}}^{ES}, E^{ES, I}

Minimum, maximum and initial stored energy of ESS (p.u)

{\underline{E}}^{TS}, {\bar{E}}^{TS}, E^{TS, I}

Minimum, maximum and initial stored energy of TSS (p.u)

{\bar{H}}^{B}

Heating capacity of the boiler (p.u)

{\bar{H}}^{TS, C} / {\bar{H}}^{TS, D}

TSS charge/discharge rate (p.u)

HV^B, HV^C

Heating capacity of natural gas in the boiler and CHP (kWh/m³)

K_Q

Rate of reactive power and electrical energy prices

P^D, Q^D, G^D, H^D

Active, reactive, gas and heating load (p.u)

{\bar{P}}^{ES, C} / {\bar{P}}^{ES, D}

ESS charge/discharge rate (p.u)

{\bar{P}}^{EV, C}, {\bar{P}}^{EV, D}

EVs charge/discharge rate in a parking lot (p.u)

P^R

Active power of renewable energy source (RES) (p.u)

{\bar{S}}^{C}, {\bar{H}}^{C}

CHP electrical and heating capacity (p.u)

{\bar{S}}^{ES}

{\bar{S}}^{EV}

ESS charger capacity, and all EVs charger capacity in a parking lot (p.u)

sign(π_g, π_j)

Gas pipeline constant, i.e. it is 1/-1 if π_g is more/less than π_j.

{\bar{S}}^{L}, {\bar{G}}^{L}, {\bar{H}}^{L}

Electrical, gas and heating line capacity (p.u)

{\bar{S}}^{S}, {\bar{G}}^{S}, {\bar{H}}^{S}

Electrical, gas and heating station capacity (p.u)

\bar{T}, \underline{T}

Maximum and minimum values of temperature (p.u)

\bar{V}, \underline{V}

Maximum and minimum values of voltage magnitude (p.u)

η ^B

Efficiency of boiler

η^ES,D/η^ES,C

The efficiency of ESS discharging/charging mode

η^EV,D/η^EV,C

The efficiency of EVs discharging/charging mode

η^T, η^L, η^H

Turbine, loss and heating efficiency of CHP

η^TS,D/η^TS,C

The efficiency of TSS discharging/charging mode

λ^E, λ^G, λ^H

DA electricity, gas and heating energy price ($/MWh)

ρ^UE, ρ^DE

DA up and down reserve market price for electrical reserve requirements ($/MWh)

ρ^UH, ρ^DH

DA up and down reserve market price for heating reserve requirements ($/MWh)

Constant of gas pipeline (p.u)

\bar{π}, \underline{π}

Maximum and minimum values of pressure (p.u)

1 Introduction

Nowadays, to reduce environmental pollution, the use of renewable and non-renewable distributed generations (DGs), and storage systems in the power system has grown significantly [1, 2]. Nevertheless, the high penetration level of sources and ESs complicates the system operator’s decision to manage the network [3]. To cope with this issue, the smart grid theory suggests using coordination schemes between sources and storages in an aggregate framework such as a microgrid, virtual power plant (VPP), and energy hub [4]. Note that the energy hub scheme includes optimal potential in the energy efficiency with respect to other outlines proposed by smart grid theory due to simultaneous management of several types of energy [5, 6]. The energy hub is defined as “a multi-carrier energy system consisting of multiple energy conversion, storage and/or network technologies, and characterized by some degree of local control” [7]. VPP and microgrid are like energy hub, except that they consider only the electrical energy equipment. Generally, microgrid is based on multi-bus network model, while VPP includes single-level network framework [4]. Therefore, it is predicted that the energy hub for coordination and energy management of different sources and storages can have suitable performance in improving the operation, flexibility, reliability and security indices of several energy networks such as electricity, gas, and district heating systems. Furthermore, it can participate in energy and auxiliary service markets to achieve optimal profit for its sources and storages.

Some various works and researches investigate the capabilities and model of the hub in energy networks and different markets. The authors of [8] incorporate several energy hubs in the distribution network, where hubs are participating in the energy market. According to the results of [8], hubs can obtain high profit for their sources and storages while reducing the distribution network operation cost. In [9], the energy hub concept is investigated in the grid-connected microgrids scheme to provide the potential capabilities of this grid in satisfying various types of energy demands. In [10], a bi-level problem is proposed for optimal interaction between energy hubs and the distribution network, in which total operation cost of the network is minimized and it limits to network’s model in the upper-level formulation, and operation cost of energy hubs connected to the distribution system is minimized in the lower level formulation. Moreover, the coupled multiple energy hubs including energy converts, renewable energy sources, and storages is modeled in [11], where it uses a peer-to-peer (P2P) transaction framework to improve economic indices of the hubs. In [12], the impacts of ice storage on energy hub operation are investigated so that it can reduce the operation cost of the hub. In addition, the optimal offering model of the price-maker energy hub is formulated in [13], where it uses the game theory to achieve the optimal bidding strategy in a competitive and equilibrium market. In [14], the authors present a cooperative trading method to a community-level energy system that consists of energy hubs and photovoltaic (PV) systems, as well as demand response (DR) program. According to [14], the cooperative trading method contains a real-time rolling horizon energy management formulation that follows the cooperative game theory considering the uncertainty of PV generation power and risk model. Also, the hub model in [14] with compressed air energy storage (CAES) is expressed in [15], where it minimizes the operating cost of different sources, storages and programs. Moreover, the market model of the energy hub in the day-ahead (DA) electricity and heat energy markets is formulated in [16] so that the results of [16] demonstrate the hub’s capability to obtain high profit for its sources and storages. Moreover, the real-time operation of hubs according to the dynamic pricing market model is presented in [17] so that hubs interaction is formulated as an accurate potential game to reduce hubs cost in the electricity and gas markets, and the customers’ satisfaction from energy consumption. In [18], the grid-connected energy hub model in the electricity, gas, and district heating systems according to their participation in the DA energy market is presented, where it employs the coordinated energy management (CEM) method to hubs sources and storages.

According to the literature review, there are three main research gaps concerning the energy hub optimization model in the energy networks and different markets as follows:

In general, the participation of hubs in the DA energy market has been modeled and investigated. However, noting that a hub includes different sources and storages, hubs are expected to enjoy different capabilities in energy networks to improve operation, flexibility, reliability, security and economic indices. Therefore, it is suitable for a hub to participate in the auxiliary service markets, such as the reserve market, to increase its profit based on its potential.

Most studies mention that an energy hub that uses the coordination between different sources and storages can improve the flexibility of the system. Yet, the flexibility model of a hub has been expressed in rare research works. However, the flexibility of the system should be measured to investigate the impacts of enhanced flexibility on other indices such as operation and economic situation in the hub or different energy networks.

It is predicted that electric vehicles (EVs) parking lot will be included in the hub in the future years; hence, impacts of management and non-management model of EVs energy on operation indices of the distribution network has been investigated in many studies. Nonetheless, EVs parking lot is expected to be a suitable flexible source for improving electrical network flexibility in the presence of renewable energy sources (RESs), a case which has rarely been examined in the literature. Furthermore, while CHPs are assumed as flexible sources in the electricity section, their thermal power depends on active power; hence, the flexibility of the thermal section of the hub may be small. To tackle this issue, thermal storage devices can be utilized as flexibility sources. This also has rarely been addressed in the literature.

To cope with the above issues, this paper models the CEM strategy for different sources and storages in the scheme of a flexible grid-connected hub that can participate in DA energy and reserve markets, as shown in Fig. 1. Note that in the proposed hub, the combined heat and power (CHP) system, electrical storage system (ESS) and the EVs parking lot are considered as flexible sources to compensate for the variation of RES generation active power. Also, boiler and thermal storage system (TSS) is used in the hub to compensate for the variations in the CHP generation heating power. Besides, the proposed problem maximizes the difference between the expected revenue of hubs in DA energy and reserve markets and the expected cost of lost flexibility (COLF). Accordingly, this study considers the COLF as a flexibility model of the hub to address the second research gap. Moreover, this problem is subject to linearized optimal power flow (LOPF) equations in the electricity, gas, and district heating systems, a hub model containing the constraints of the proposed sources and storages and formulation of the hub reserve. It should be noted that the proposed strategy includes uncertainties of the load, reserve requirement, market price, RES generation power and EVs demand. Hence, the hybrid stochastic/robust optimization (HSRO) is utilized in this paper to model these uncertainty parameters and achieve robust capability for the hub to enhance system flexibility. Based on these objectives, the bounded uncertainty-based robust optimization (BURO) models EVs uncertainty parameters because the EVs parking lot is the only flexible source subject to uncertainty. Also, the scenario-based stochastic programming (SBSP) models other uncertainty parameters according to the hybrid approach of Monte Carlo Simulation (MCS)-based scenario generation method and the synchronous backward method (SBM)-base scenario reduction technique. Finally, the main contributions of the paper are summarized as follows:

Implementing the CEM strategy in the scheme of flexible grid-connected energy hub

Modeling the flexibility index for the hub using the expected cost of lost flexibility model;

Evaluating benefits of the grid-connected hubs participating in DA energy and reserve markets; and

Achieving robust potential for a hub in the system flexibility improving.

Fig. 1

The scheme of the grid connected hubs partnership in the DA markets.

The rest of the paper is organized as follows: Section 2 expresses the linearized model to the problem of the hub participation in DA energy and reserve markets. Section 3 obtains the HSRO formulation for the proposed method. Sections 4 and 5 demonstrate numerical simulation results and main conclusions of the proposed approach, respectively.

2 Model of hub participation in the DA energy and reserve markets

In this section, the mathematical model of the grid-connected energy hubs’ participation in DA energy and reserve markets is presented. This problem maximizes the difference between the expected hubs revenue (HR) in the proposed markets and the expected COLF for hubs, and also, it is limited to electricity, gas, and district heating systems constraints and hub model. Therefore, the proposed strategy formulation is written as follows:

Objective function: The proposed problem considers the maximization of the difference between HR in energy and reserve markets and hubs COLF as an objective function that is formulated in (1). The first part refers to the expected HR, where it includes three terms. The first term obtains the HR in the DA energy market due to electrical, gas and heating power exchange between the hub and this market. The second/third term calculates the HR in the DA up/down reserve market arising from electrical and heating reserve power exchange of hub with up/down reserve market. In this part, it is assumed that the gas demand is supplied by the gas station, or local gas sources are not used in the gas network; hence, there is no reserve term for gas power in equation (1). Also, if hub power in different energy grids is positive, the hub will receive revenue from DA markets, otherwise, the hub includes the cost. Besides, the second term of equation (1) refers to the expected COLF of hubs. Note that the flexibility term is defined as “the modification of generation injection and/or consumption patterns in reaction to an external price or activation signal in order to provide a service within the electrical system” [19]. Thus, to improve the hub flexibility, the difference between hub power in different scenarios and hub power in the scenario related to the forecasted value of uncertain parameters (denoted by “1” in equation (1)) should be close to zero. Therefore, this statement can be formulated as the second part of (1), where VOLF in this equation refers to the value of lost flexibility. Also, RES/CHP eliminates the flexibility in the electrical/district heating system; hence, COLF is repeated for electrical and district heating grids. $\begin{matrix} max Profit = \\ \sum_{ω \in S} ρ_{ω} \sum_{t \in ST} \sum_{m \in H} [\begin{matrix} {\begin{matrix} λ_{t, ω}^{E} (P_{m, t, ω} + K_{Q} | Q_{m, t, ω} |) + \\ λ_{t, ω}^{G} G_{m, t, ω} + λ_{t, ω}^{H} H_{m, t, ω} \end{matrix}} + \\ {λ_{t, ω}^{UE} {RU}_{m, t, ω}^{E} + λ_{t, ω}^{UH} {RU}_{m, t, ω}^{H}} + \\ {λ_{t, ω}^{DE} {RD}_{m, t, ω}^{E} + λ_{t, ω}^{DH} {RD}_{m, t, ω}^{H}} \end{matrix}] \\ - \sum_{ω \in S} ρ_{ω} \sum_{t \in ST} \sum_{m \in H} VOLF [\begin{matrix} | P_{m, t, ω} - P_{m, t, 1} | \\ + | H_{m, t, ω} - H_{m, t, 1} | \end{matrix}] \end{matrix}$ (1)

Energy network constraints: The proposed scheme is based on the optimal power flow (OPF) model in the electricity, gas, and district heating systems. The original OPF equation in different systems is a non-convex NLP, hence, this paper uses the linearized OPF model for these networks to improve/obtain computational time/the optimal solution that is unique for all solvers. Accordingly, the linearized AC power flow equation in the electrical network is formulated in (2a)-(2e) that represent respectively the nodal active and reactive power balance, active and reactive power flowing from the distribution network, and voltage angle value on the slack bus [20 –23]. Also, the gas power flow model is written in (2f)-(2j), where equation (2f) refers to nodal gas power balance and constraint (2g) presents the amount of gas flow through a pipeline [18]. In the original gas network model, an auxiliary variable, δ, is calculated by (δ_g,j) ² = (π_g) ² - (π_j) ². The linear format of this equation is as (2 h) based on the conventional piecewise linearization method [18], where Δδ <<1 and Δπ <<1 are the deviation of δ and π with respect to 1 p.u. Hence, constraint (2i)/(2j) expresses equation δ/π according to variables Δδ/Δπ. In addition, the heating power flow equations in the district heating network are based on (2k) and (2l), which introduce the nodal heating power balance and heating power flowing from the pipeline, respectively [18].

In constraints (2 m)-(2u), operation limits of electricity, gas and district heating systems are expressed so that bus voltage deviation limit and capacity limitation of the distribution line and station are presented in (2 m)-(2o) [20], respectively. Also, gas nodes pressures, gas pipeline capacity and gas station capacity limitations are formulated in (2p)-(2r), and limits of node temperature, heating pipeline and station are modeled in (2 s)-(2u) [18]. Finally, the up and down reserve demands on the electrical and district heating network are provided by grid-connected hubs according to constraints (2v)-(2 w) and (2x)-(2y), respectively. Note that the capacity limit of the distribution line and station is based on the circular plane with a radius of S, i.e. $\sqrt{P^{2} + Q^{2}} ⩽ S$ [21]. However, a regular polygon plane that is equivalent to a circular plane is used for these limits as constraints (2n) and (2o) to present a linear formulation, where more details of this linearization method are expressed in [20]. $\begin{matrix} Constraints of Energy Networks : \\ {\begin{matrix} P_{e, t, ω}^{S} + \sum_{m \in H} A_{e, m}^{E} P_{m, t, ω} - \sum_{j \in EN} B_{e, j}^{E} P_{e, j, t, ω}^{L} = P_{e, t, ω}^{D} \\ \forall e, t, ω \end{matrix} \end{matrix}$ (2a) $\begin{matrix} Q_{e, t, ω}^{S} + \sum_{m \in H} A_{e, m}^{E} Q_{m, t, ω} - \sum_{j \in EN} B_{e, j}^{E} Q_{e, j, t, ω}^{L} = Q_{e, t, ω}^{D} \\ \forall e, t, ω \end{matrix}$ (2b) $\begin{matrix} P_{e, j, t, ω}^{L} = {GL}_{e, j} (Δ V_{e, t, ω} - Δ V_{j, t, ω}) - \\ {BL}_{e, j} (θ_{e, t, ω} - θ_{j, t, ω}) \forall e, j, t, ω \end{matrix}$ (2c) $\begin{matrix} Q_{e, j, t, ω}^{L} = - {BL}_{e, j} (Δ V_{e, t, ω} - Δ V_{j, t, ω}) - \\ {GL}_{e, j} (θ_{e, t, ω} - θ_{j, t, ω}) \forall e, j, t, ω \end{matrix}$ (2d) $θ_{e, t, ω} = 0 \forall e = Slack bus, t, ω$ (2e) $\begin{matrix} G_{g, t, ω}^{S} + \sum_{m \in H} A_{g, m}^{G} G_{m, t, ω} - \sum_{j \in GN} B_{g, j}^{G} G_{g, j, t, ω}^{L} = \\ G_{g, t, ω}^{D} \forall g, t, ω \end{matrix}$ (2f) $G_{g, j, t, ω}^{L} = κ_{g, j} sign (π_{g}, π_{j}) δ_{g, j, t, ω} \forall g, j, t, ω$ (2g) $\begin{matrix} Δ δ_{g, j, t, l, ω} = \\ sign (π_{g}, π_{j}) \frac{m_{l}^{π}}{m_{l}^{δ}} (Δ π_{g, t, l, ω} - Δ π_{j, t, l, ω}) \\ \forall g, j, t, l, ω \end{matrix}$ (2h) $δ_{g, j, t, ω} = \sum_{l \in L} Δ δ_{g, j, t, l, ω} \forall g, j, t, ω$ (2i) $π_{g, t, ω} = {\underline{π}}_{g} + \sum_{l \in L} Δ π_{g, t, l, ω} \forall g, t, ω$ (2j) $\begin{matrix} H_{h, t, ω}^{S} + \sum_{m \in H} A_{h, m}^{H} H_{m, t, ω} - \sum_{j \in HN} B_{h, j}^{H} H_{h, j, t, ω}^{L} = \\ H_{h, t, ω}^{D} \forall h, t, ω \end{matrix}$ (2k) $H_{h, j, t, ω}^{L} = c_{h, j} {\dot{m}}_{h, j} (T_{h, t, ω} - T_{j, t, ω}) \forall h, j, t, ω$ (2l) ${\underline{V}}_{e} - 1 ⩽ Δ V_{e, t, ω} ⩽ {\bar{V}}_{e} - 1 \forall e, t, ω$ (2m) $\begin{matrix} P_{e, j, t, ω}^{L} cos (k \times Δ α) + Q_{e, j, t, ω}^{L} sin (k \times Δ α) ⩽ {\bar{S}}_{e, j}^{L} \\ \forall e, j, t, k, ω \end{matrix}$ (2n) $\begin{matrix} P_{e, t, ω}^{S} cos (k \times Δ α) + Q_{e, t, ω}^{S} sin (k \times Δ α) ⩽ {\bar{S}}_{e}^{S} \\ \forall e = Slack bus, t, k, ω \end{matrix}$ (2o) ${\underline{π}}_{g} ⩽ π_{g, t, ω} ⩽ {\bar{π}}_{g} \forall g, t, ω$ (2p) $- {\bar{G}}_{g, j}^{L} ⩽ G_{g, j, t, ω}^{L} ⩽ {\bar{G}}_{g, j}^{L} \forall g, j, t, ω$ (2q) $- {\bar{G}}_{g}^{S} ⩽ G_{g, t, ω}^{S} ⩽ {\bar{G}}_{g}^{S} \forall g = Slack bus, t, ω$ (2r) ${\underline{T}}_{h} ⩽ T_{h, t, ω} ⩽ {\bar{T}}_{h} \forall h, t, ω$ (2s) $- {\bar{H}}_{h, j}^{L} ⩽ H_{h, j, t, ω}^{L} ⩽ {\bar{H}}_{h, j}^{L} \forall h, j, t, ω$ (2t) $- {\bar{H}}_{h}^{S} ⩽ H_{h, t, ω}^{S} ⩽ {\bar{H}}_{h}^{S} \forall h = Slack bus, t, ω$ (2u) $\sum_{m \in H} {RU}_{m, t, ω}^{E} ⩾ {DU}_{t, ω}^{E} \forall t, ω$ (2v) $\sum_{m \in H} {RD}_{m, t, ω}^{E} ⩾ {DD}_{t, ω}^{E} \forall t, ω$ (2w) $\sum_{m \in H} {RU}_{m, t, ω}^{H} ⩾ {DU}_{t, ω}^{H} \forall t, ω$ (2x) $\sum_{m \in H} {RD}_{m, t, ω}^{H} ⩾ {DD}_{t, ω}^{H} \forall t, ω}$ (2y)

Hub constraints: In this paper, the hub contains a RES, a CHP, an electricity storage system, (ESS), a thermal storage system (TSS), an EVs parking lot and a boiler. Accordingly, the active power of the hub as equation (3a) is equal to the difference between CHP, RES, ESS and EVs generation active power and ESS, EVs and passive load demands. Based on (3b), the hub reactive load is provided by CHP generator and ESS and EVs chargers, and the excess reactive power of these sources is considered to hub reactive power that is injected into the electrical network. A hub is a gas consumer according to (3c) to supply CHP, boiler and passive gas load, and hub heating power as (3d) is equal to the difference between CHP, boiler and TSS generation heating power and passive heating load. In addition, the model of different sources and storages is expressed in (3e)-(3aa). The CHP model is presented in (3e)-(3 h) that refer respectively to CHP heating and gas power equations based on CHP active power, generation electrical and heating power limits in CHP [24]. Also, boiler constraints are formulated in (3i) and (3 j), where its input gas demand based on its output heating power is calculated by (3i), and its capacity limit to product heating power is expressed in (3 j). Moreover, TSS constraints are modeled in (3k)-(3o), which present respectively stored energy equation in TSS, TSS charge and discharge rate limits, initial stored energy value in TSS, and TSS stored energy limitation. The ESS model is demonstrated in (3p)-(3u), where equations (3p)-(3t) are the same as the TSS model and constraint (3u) denotes the ESS charger capability curve limit [20]. Finally, EVs parking lot formulation is given in (3v)-(3aa) [25 –27], where stored energy in the EVs battery at the parking lot is obtained from equation (3v). Also, charge and discharge rate limits of all EV batteries in the parking lot follow equations (3 w) and (3x), respectively. Initial/final stored energy in EVs at arrival/departure time is according to (3y)/(3z), and the capacity limit of all EV chargers in the parking lot is presented in (3aa). In this equation, the values of ${\bar{P}}^{EV, C}$ , ${\bar{P}}^{EV, D}$ and ${\bar{S}}^{EV}$ at time t are equal to $\sum_{v = 1}^{N_{t}} {CR}_{v}$ , $\sum_{v = 1}^{N_{t}} {DR}_{v}$ and $\sum_{v = 1}^{N_{t}} {CC}_{v}$ , respectively. CR/DR and CC are charge/discharge rate of EV battery and EV charger capacity, and N_t refers to the total EVs number in the parking lot at hour t. The values of E^A and E^D are obtained by equations $\sum_{v = 1}^{N_{A}} {SOC}_{v} {BC}_{v}$ and $\sum_{v = 1}^{N_{D}} {BC}_{v}$ , respectively. In these equations, SOC/BC is the state of charge (SOC)/battery capacity defined to EV battery, and N_A/N_D is the total number of EVs in the parking lot at arrival/departure time [23].

Reserve model of the hub is expressed in (3bb)-(3ff), where hub up/down reserve power in the electrical system is calculated by constraint (3bb)/(3cc) and in the district heating system, it is achieved by (3dd)/(3ee). Note the hub reserve power is positive according to (3ff), and the maximum/minimum capacity of the hub to inject/absorb different power depends on the maximum/minimum capacity of its local sources and storages, where this statement is created in (3bb)-(3ee). $\begin{matrix} Constraints of Energy Hub : \\ {\begin{matrix} P_{m, t, ω} = P_{m, t, ω}^{C} + P_{m, t, ω}^{R} + (P_{m, t, ω}^{ES, D} - P_{m, t, ω}^{ES, C}) \\ + (P_{m, t, ω}^{EV, D} - P_{m, t, ω}^{EV, C}) - P_{m, t, ω}^{D} \forall m, t, ω \end{matrix} \end{matrix}$ (3a) $\begin{matrix} Q_{m, t, ω} = Q_{m, t, ω}^{C} + Q_{m, t, ω}^{ES} + Q_{m, t, ω}^{EV} - Q_{m, t, ω}^{D} \\ \forall m, t, ω \end{matrix}$ (3b) $G_{m, t, ω} = - G_{m, t, ω}^{C} - G_{m, t, ω}^{B} - G_{m, t, ω}^{D} \forall m, t, ω$ (3c) $\begin{matrix} H_{m, t, ω} = H_{m, t, ω}^{C} + H_{m, t, ω}^{B} + (H_{m, t, ω}^{TS, D} - H_{m, t, ω}^{TS, C}) \\ - H_{m, t, ω}^{D} \forall m, t, ω \end{matrix}$ (3d) $H_{m, t, ω}^{C} = \frac{P_{m, t, ω}^{C} (1 - η_{m}^{T} - η_{m}^{L})}{η_{m}^{T}} . η_{m}^{H} \forall m, t, ω$ (3e) $G_{m, t, ω}^{C} = \frac{P_{m, t, ω}^{C}}{η_{m}^{T} {HV}_{m}^{C}} \forall m, t, ω$ (3f) $\begin{matrix} P_{m, t, ω}^{C} cos (k \times Δ α) + Q_{m, t, ω}^{C} sin (k \times Δ α) ⩽ {\bar{S}}_{m}^{C} \\ \forall m, t, k, ω \end{matrix}$ (3g) $0 ⩽ H_{m, t, ω}^{C} ⩽ {\bar{H}}_{m}^{C} \forall m, t, ω$ (3h) $G_{m, t, ω}^{B} = \frac{H_{m, t, ω}^{B}}{η_{m}^{B} {HV}_{m}^{B}} \forall m, t, ω$ (3i) $0 ⩽ H_{m, t, ω}^{B} ⩽ {\bar{H}}_{m}^{B} \forall m, t, ω$ (3j) $\begin{matrix} E_{m, t + 1, ω}^{TS} = E_{m, t, ω}^{TS} + η^{TS, C} H_{m, t, ω}^{TS, C} - \frac{1}{η^{TS, D}} H_{m, t, ω}^{TS, D} \\ \forall m, t, ω \end{matrix}$ (3k) $0 ⩽ H_{m, t, ω}^{TS, C} ⩽ {\bar{H}}_{m}^{TS, C} {ts}_{m, t} \forall m, t, ω$ (3l) $0 ⩽ H_{m, t, ω}^{TS, D} ⩽ {\bar{H}}_{m}^{TS, D} (1 - {ts}_{m, t}) \forall m, t, ω$ (3m) $E_{m, t, ω}^{TS} = E_{m}^{TS, I} \forall m, t = 1, ω$ (3n) ${\underline{E}}_{m}^{TS} ⩽ E_{m, t, ω}^{TS} ⩽ {\bar{E}}_{m}^{TS} \forall m, t, ω$ (3o) $\begin{matrix} E_{m, t + 1, ω}^{ES} = E_{m, t, ω}^{ES} + η^{ES, C} P_{m, t, ω}^{ES, C} - \frac{1}{η^{ES, D}} P_{m, t, ω}^{ES, D} \\ \forall m, t, ω \end{matrix}$ (3p) $0 ⩽ P_{m, t, ω}^{ES, C} ⩽ {\bar{P}}_{m}^{ES, C} {es}_{m, t} \forall m, t, ω$ (3q) $0 ⩽ P_{m, t, ω}^{ES, D} ⩽ {\bar{P}}_{m}^{ES, D} (1 - {es}_{m, t}) \forall m, t, ω$ (3r) $E_{m, t, ω}^{ES} = E_{m}^{ES, I} \forall m, t = 1, ω$ (3s) ${\underline{E}}_{m}^{ES} ⩽ E_{m, t, ω}^{ES} ⩽ {\bar{E}}_{m}^{ES} \forall m, t, ω$ (3t) $\begin{matrix} (P_{m, t, ω}^{ES, D} - P_{m, t, ω}^{ES, C}) cos (k \times Δ α) + \\ Q_{m, t, ω}^{ES} sin (k \times Δ α) ⩽ {\bar{S}}_{m}^{ES} \forall m, t, k, ω \end{matrix}$ (3u) $\begin{matrix} E_{m, t + 1, ω}^{EV} = E_{m, t, ω}^{EV} + η^{EV, C} P_{m, t, ω}^{EV, C} - \frac{1}{η^{EV, D}} P_{m, t, ω}^{EV, D} \\ \forall m, t, ω \end{matrix}$ (3v) $0 ⩽ P_{m, t, ω}^{EV, C} ⩽ {\bar{P}}_{m, t, ω}^{EV, C} {ev}_{m, t} \forall m, t, ω$ (3w) $0 ⩽ P_{m, t, ω}^{EV, D} ⩽ {\bar{P}}_{m, t, ω}^{EV, D} (1 - {ev}_{m, t}) \forall m, t, ω$ (3x) $E_{m, t, ω}^{EV} = E_{m, t, ω}^{A} \forall m, t = Arrival time, ω$ (3y) $E_{m, t, ω}^{EV} = E_{m, t, ω}^{D} \forall m, t = Departure time, ω$ (3z) $\begin{matrix} (P_{m, t, ω}^{EV, D} - P_{m, t, ω}^{EV, C}) cos (k \times Δ α) + \\ Q_{m, t, ω}^{EV} sin (k \times Δ α) ⩽ {\bar{S}}_{m, t, ω}^{EV} \forall m, t, k, ω \end{matrix}$ (3aa) $\begin{matrix} P_{m, t, ω} + {RU}_{m, t, ω}^{E} ⩽ {\bar{S}}_{m}^{C} + {\bar{P}}_{m}^{ES, D} + {\bar{P}}_{m, t, ω}^{EV, D} + P_{m, t, ω}^{R} \\ \forall m, t, ω \end{matrix}$ (3bb) $P_{m, t, ω} - {RD}_{m, t, ω}^{E} ⩾ - {\bar{P}}_{m}^{ES, C} - {\bar{P}}_{m, t, ω}^{EV, C} \forall m, t, ω$ (3cc) $\begin{matrix} H_{m, t, ω} + {RU}_{m, t, ω}^{H} ⩽ {\bar{H}}_{m}^{C} + {\bar{H}}_{m}^{B} + {\bar{H}}_{m}^{TS, D} \\ \forall m, t, ω \end{matrix}$ (3dd) $H_{m, t, ω} - {RD}_{m, t, ω}^{H} ⩾ - {\bar{H}}_{m}^{TS, C} \forall m, t, ω$ (3ee) $\begin{matrix} {RU}_{m, t, ω}^{E}, {RD}_{m, t, ω}^{E}, {RU}_{m, t, ω}^{H}, {RD}_{m, t, ω}^{H} ⩾ 0 \\ \forall m, t, ω \end{matrix}}$ (3ff) It should be noted that the objective function (1) is non-linear due to the absolute terms. To achieve a linear format for equation (1), the proposed problem is rewritten as follows: $\begin{matrix} max Profit = \\ \sum_{ω \in S} ρ_{ω} \sum_{t \in ST} \sum_{m \in H} [\begin{matrix} {\begin{matrix} λ_{t, ω}^{E} (P_{m, t, ω} + K_{Q} | Q_{m, t, ω} |) + \\ λ_{t, ω}^{G} G_{m, t, ω} + λ_{t, ω}^{H} H_{m, t, ω} \end{matrix}} + \\ {λ_{t, ω}^{UE} {RU}_{m, t, ω}^{E} + λ_{t, ω}^{UH} {RU}_{m, t, ω}^{H}} + \\ {λ_{t, ω}^{DE} {RD}_{m, t, ω}^{E} + λ_{t, ω}^{DH} {RD}_{m, t, ω}^{H}} \end{matrix}] \\ - \sum_{ω \in S} ρ_{ω} \sum_{t \in ST} \sum_{m \in H} VOLF [\begin{matrix} F_{m, t, ω}^{E +} + F_{m, t, ω}^{E -} + \\ F_{m, t, ω}^{H +} + F_{m, t, ω}^{E -} \end{matrix}] \end{matrix}$ (4a) Subject to: $\begin{matrix} F_{m, t, ω}^{E +} - F_{m, t, ω}^{E -} = P_{m, t, ω} - P_{m, t, 1} \\ \forall m, t, ω, F_{m, t, ω}^{E +}, F_{m, t, ω}^{E -} ⩾ 0 \end{matrix}$ (4b) $\begin{matrix} F_{m, t, ω}^{H +} - F_{m, t, ω}^{H -} = H_{m, t, ω} - H_{m, t, 1} \\ \forall m, t, ω, F_{m, t, ω}^{H +}, F_{m, t, ω}^{H -} ⩾ 0 \end{matrix}$ (4c) $Constraints (2) and (3)$ (4d)

According to (4b)/(4c)/(4d), the term P_ω - P₁/H_ω - H₁/Q can be equal to the difference between two positive variables of F^E + and F^E - (F^H + and F^H -) / (Q+ and Q^-); hence, the absolute value of P_ω - P₁/H_ω - H₁/Q is equal to the summation of F^E + and F^E - (F^H + and F^H -) / (Q+ and Q^-). Also, variables F+ and F^- (Q+ and Q^-) are defined as upward and downward flexibility power (capacitive and inductive reactive power) based on [23, 28]. Noted that, there is upward flexibility power if active (heating) power in scenario ω is greater than in scenario 1; otherwise, it is downward flexibility power.

3 Proposed HSRO model

In the proposed strategy, (1)-(3), load (P^D, Q^D, H^D and G^D), market price (λ^E, λ^G, λ^H, λ^UE, λ^UH, λ^DE and λ^DH), reserve demand (DU^E, DD^E, DU^H and DD^H), RES generation power (P^R), and EVs parking lot parameters ( ${\bar{P}}^{EV, C}$ , ${\bar{P}}^{EV, D}$ , ${\bar{S}}^{EV}$ , E^A, E^D) are uncertain parameters. Also, the robust capability of the hub to improve system flexibility is one of the main objectives in this paper. Hence, the hybrid stochastic/robust optimization (HSRO) method is used to model the proposed uncertainty parameters. Among the flexibility sources in the hub, i.e., EVs, ESS, TSS, CHP, boiler, the EVs parking lot contains uncertainty parameters. Therefore, robust programming is used to model parameters ${\bar{P}}^{EV, C}$ , ${\bar{P}}^{EV, D}$ , ${\bar{S}}^{EV}$ , E^A, E^D, where it can achieve the robust capability of the hubs in increasing system flexibility. In this paper, the bounded uncertainty-based robust optimization (BURO) is expressed to this issue, where the true value of the uncertain parameter in this method is equal to (1+γ) $\bar{u}$ or (1 - γ) $\bar{u}$ based on the location of this parameter in the optimization problem [20]. $\bar{u}$ is the forecasted value of the uncertain parameter and γ ≥ 0 is the uncertainty level [20]. Other uncertain parameters in this paper are modeled using scenario-based stochastic programming (SBSP). In this approach, the Monte Carlo Simulation (MCS) generates numerous scenario samples based on the normal probability distribution function (PDF) for uncertainties of load, energy price and reserve demand, and it uses the Weibull and Beta PDF for the uncertainty of RES generation power. In the next step, the SBM is utilized to decrease the number of scenarios and enhance the tractability of the proposed MILP operation tool [29]. Therefore, the proposed HSRO model is written as follows: $Objective function (4 a)$ (5a) Subject to: $0 ⩽ P_{m, t, ω}^{EV, C} ⩽ (1 - γ) {\bar{P}}_{m, t}^{EV, C} {ev}_{m, t} \forall m, t, ω$ (5b) $0 ⩽ P_{m, t, ω}^{EV, D} ⩽ (1 - γ) {\bar{P}}_{m, t}^{EV, D} (1 - {ev}_{m, t}) \forall m, t, ω$ (5c) $E_{m, t, ω}^{EV} = (1 + γ) E_{m, t}^{A} \forall m, t = Arrival time, ω$ (5d) $E_{m, t, ω}^{EV} = (1 + γ) E_{m, t}^{D} \forall m, t = Departure time, ω$ (5e) $\begin{matrix} (P_{m, t, ω}^{EV, D} - P_{m, t, ω}^{EV, C}) cos (k \times Δ α) + \\ Q_{m, t, ω}^{EV} sin (k \times Δ α) ⩽ (1 - γ) {\bar{S}}_{m, t}^{EV} \forall m, t, k, ω \end{matrix}$ (5f) $\begin{matrix} P_{m, t, ω} + {RU}_{m, t, ω}^{E} ⩽ {\bar{S}}_{m}^{C} + {\bar{P}}_{m}^{ES, D} + \\ {\bar{P}}_{m, t}^{EV, D} (1 - γ) + P_{m, t, ω}^{R} \forall m, t, ω \end{matrix}$ (5g) $\begin{matrix} P_{m, t, ω} - {RD}_{m, t, ω}^{E} ⩾ - {\bar{P}}_{m}^{ES, C} - {\bar{P}}_{m, t}^{EV, C} (1 - γ) \\ \forall m, t, ω \end{matrix}$ (5h) $Constraints (2), (3 a) - (3 v), and (3 dd) - (3 ff)$ (5i)

In problem (5), the objective function (5a) is the same as Equation (4a), constraints (5b)–(5h) refer to the robust or BURO model of (3 w)-(3cc) according to the true value of uncertainty of ${\bar{P}}^{EV, C}$ , ${\bar{P}}^{EV, D}$ , ${\bar{S}}^{EV}$ , E^A, E^D in the worst-case scenario. Note that in this scenario or model, the feasible region is less than that in other scenarios; thus, the true value of the EVs uncertain parameter according to (5b)-(5 h) is equal to its up-limit ((1+γ) $\bar{u}$ ) or down-limit ((1 - γ) $\bar{u}$ ) due to the linear format of constraints (3 w)-(3cc) [20]. This case depends on the min/max term of the objective function, the positive or negative coefficient of the uncertain parameter, the location of the uncertain parameter in the objective function or constraint. Accordingly, maximum profit will be obtained if the active power of the hub (P) includes high value, where this condition can occur if EVs demand reduces. Therefore, in the worst-case scenario, the true values of E^A and E^D ( ${\bar{P}}^{EV, C}$ , ${\bar{P}}^{EV, D}$ and ${\bar{S}}^{EV}$ ) are equal to its up-limit (down-limit) to achieve the small feasible region with respect to other scenarios. Moreover, EVs uncertain parameters are created in one scenario, hence, index ω is removed from these uncertainty parameters in constraints (5b)-(5 h). Simultaneously, the SBSP method is used in problem (5) to model the uncertainty of load, market price, reserve demand, and RES generation power. Finally, constraint (5i) is the same as equations of formulation (1)–(3) that do not include EVs uncertainty parameters.

4 Numerical results and discussion

4.1 Case study

The proposed strategy is implemented on the standard 9-bus electrical network, 4-node natural gas grid and 7-node district heating system as plotted in Fig. 2 [18]. The characteristics of the distribution lines and post, gas and heating pipelines and station, and peak electricity and heating load data are expressed in [18], where the gas load is considered zero. This system includes base power of 1 MW, the base voltage of 1 kV, the base pressure of 10 bar, and the base temperature of 100 C⁰, and also, the allowed range of bus voltage and node pressures and temperature are considered to be in the range [0.9, 1.1] p.u. [18]. The daily curves of the up and down reserve requirements or demand in the electrical [30] and district heating networks are shown in Figs. 3(a) and (3b), respectively. Also, the daily DA energy price curve related to proposed different energy networks is presented in [18] so that gas energy price is 12/18 $/MWh at period 1:00-4:00 and 23:00-24:00/5:00-22:00. The electricity energy price is equal to 17.6/26.4/33 $/MWh at hours 1:00-7:00/8:00-16:00 and 23:00-24:00/17:00-22:00, and heating energy price is 22/30 $/MWh at hours 1:00-4:00 and 16:00-24:00/5:00-15:00. K_Q is 0.08, and daily DA up and down reserve price curves in different networks are the same as those of the daily DA energy price curve in that network. In addition, the proposed system contains 7 hubs, where their locations are shown in Fig. 2. Hubs 1-3 and 5 include only RES, ESS and EVs parking lot, there are CHP, boiler and TSS in hub 4, and hubs 6 and 7 contain all proposed sources and storages in the problem (5). Data of different sources and storages are expressed in [18], and daily power rate curve of RES with types of wind system and photovoltaic is based on data in [31].

Fig. 2

Schematic of the proposed test case system [18].

Finally, the standard deviation for the uncertainty of load, reserve demand, the market price and RES generation power is considered 10%. In the proposed SBSP, the MCS generates 1000 scenarios for these uncertain parameters, and thus, the SBM obtain 30 scenario samples with high probability occurrence. Moreover, EVs uncertainty parameters according to problem (5) are modeled based on the BURO method.

4.2 Results

The proposed problem in (5) is simulated in GAMS 25.2, where it is solved by CPLEX [32]. Also, the total number of linearization segments in the conventional piecewise linearization method and the circular plan are 5 and 45, respectively [18]. According to [18], the maximum computational error for active and reactive power, voltage, gas power, gas pressure, and the district heating network variables is about 2.5%, 0.5%. 0.9%, 0.1% and zero, respectively.

A) CEM Capabilities for sources and FSs in the energy hub scheme: The potential of the CEM strategy that establishes coordination between sources and FSs in the energy hub framework by considering model (3) is investigated in this section. According to this purpose, the results of the uncoordinated energy management (UEM) method is presented in this section, where this method is repeated for each source or storage. In other words, it uses problem (5) considering only RES, CHP, boiler, ESS, TSS or EVs parking lot. Accordingly, the capabilities of these two methods in the proposed energy and reserve market based on the different uncertainty level of EVs parameters (γ) are reported in Table 2. According to this table, for γ=0, the total expected revenue of sources and storage in UEM strategy for DA energy and reserve markets is equal to 194.2 $ (79.8 + 40.3 + 21.2 + 57. 5 + 46.7 –51.3) and 107.6 $ (24.7 + 14.2 + 10.3 + 16.4 + 32.7 + 9.3), respectively. Yet, if these sources and storages are coordinated in a hub based on the CEM method, their revenue or hub revenue will increase respectively to 229.35$ and 131.6 $ in DA energy and reserve markets with a gain rate of 18.1% ((229.35 –194.2)/194.2) and 22.3% ((131.6 –107.6)/107.6) with respect to UEM strategy. There is this rule in the other value of γ, but note that increasing γ increases/decreases EVs energy cost/reserve revenue in the DA energy/reserve market because the EVs demand increases in this condition due to increasing EVs energy in arrival and departure hours and reducing EVs charge and discharge rates based on the problem (5). However, the revenue of other sources and storage in the proposed markets are independent of the uncertainty level of EVs parameter according to the formulation (5) and the results of Table 2. Hence, increasing γ will reduce the revenue of all sources and storages in both UEM and CEM methods.

Table 1
Expected hubs profit in the different DA markets according to different energy management methods in γ

γ 0 0.1 0.2

Energy management method DA market DA market DA market

Energy Reserve Energy Reserve Energy Reserve

UEM CHP 79.8 24.7 79.8 24.7 79.8 24.7

Boiler 40.3 14.2 40.3 14.2 40.3 14.2

TSS 21.2 10.3 21.2 10.3 21.2 10.3

RES 57.5 16.4 57.5 16.4 57.5 16.4

ESS 46.7 32.7 46.7 32.7 46.7 32.7

EVs –51.3 9.3 –58.4 7.5 –64.1 5.7

Sum 194.2 107.6 187.1 105.8 181.4 104

CEM 229.35 131.6 218.72 128.44 209.88 125.32

Gain rate (%) 18.1 22.3 16.9 21.4 15.7 20.5

γ	0	0.1	0.2
UEM	CHP	79.8	24.7	79.8	24.7	79.8	24.7
	Boiler	40.3	14.2	40.3	14.2	40.3	14.2
	TSS	21.2	10.3	21.2	10.3	21.2	10.3
	RES	57.5	16.4	57.5	16.4	57.5	16.4
	ESS	46.7	32.7	46.7	32.7	46.7	32.7
	EVs	–51.3	9.3	–58.4	7.5	–64.1	5.7
	Sum	194.2	107.6	187.1	105.8	181.4	104
CEM	229.35	131.6	218.72	128.44	209.88	125.32
Gain rate (%)	18.1	22.3	16.9	21.4	15.7	20.5

Table 2

Expected value of network indices in the different cases

Case	EEL (p.u)	GEL (p.u)	HEL (p.u)
Case 1: Power flow (γ=0, 0.1, 0.2)	5.1211	0	2.1156
Case 2: Power flow with RESs (γ=0, 0.1, 0.2)	5.4102	0	2.1156
Case 3: Hub including RES and EVs (γ=0)	4.3971	0	2.1156
Case 3: Hub including RES and EVs (γ=0.1)	4.4893	0	2.1156
Case 3: Hub including RES and EVs (γ=0.2)	4.5732	0	2.1156
Case 4: Hub including RES and ESS (γ=0, 0.1, 0.2)	4.2312	0	2.1156
Case 5: Hub including RES and CHP (γ=0, 0.1, 0.2)	4.0126	2.1183	1.5832
Case 6: Hub including RES, EVs, ESS and CHP (γ=0)	2.9511	2.1183	1.5832
Case 6: Hub including RES, EVs, ESS and CHP (γ=0.1)	3.0112	2.1294	1.5832
Case 6: Hub including RES, EVs, ESS and CHP (γ=0.2)	3.0698	2.1406	1.5832
Case 7: Hub including CHP (γ=0, 0.1, 0.2)	3.7805	2.1818	1.5312
Case 8: Hub including CHP and boiler (γ=0, 0.1, 0.2)	3.7805	2.9831	1.3703
Case 9: Hub including CHP and TSS (γ=0, 0.1, 0.2)	3.7805	2.1818	1.4708
Case 10: Hub including CHP, boiler and TSS (γ=0, 0.1, 0.2)	3.7805	2.9831	1.3176
Case 11: Proposed strategy (γ=0)	2.9511	2.9831	1.3176
Case 11: Proposed strategy (γ=0.1)	3.0112	2.9943	1.3176
Case 11: Proposed strategy (γ=0.2)	3.0698	3.0055	1.3176
Case	MVD (p.u)	MPD (p.u)	MTD (p.u)
Case 1: Power flow (γ=0, 0.1, 0.2)	0.113	0	0.111
Case 2: Power flow with RESs (γ=0, 0.1, 0.2)	0.128	0	0.111
Case 3: Hub including RES and EVs (γ=0)	0.073	0	0.111
Case 3: Hub including RES and EVs (γ=0.1)	0.075	0	0.111
Case 3: Hub including RES and EVs (γ=0.2)	0.077	0	0.111
Case 4: Hub including RES and ESS (γ=0, 0.1, 0.2)	0.069	0	0.111
Case 5: Hub including RES and CHP (γ=0, 0.1, 0.2)	0.064	0.011	0.082
Case 6: Hub including RES, EVs, ESS and CHP (γ=0)	0.057	0.011	0.082
Case 6: Hub including RES, EVs, ESS and CHP (γ=0.1)	0.059	0.011	0.082
Case 6: Hub including RES, EVs, ESS and CHP (γ=0.2)	0.061	0.011	0.082
Case 7: Hub including CHP (γ=0, 0.1, 0.2)	0.061	0.014	0.080
Case 8: Hub including CHP and boiler (γ=0, 0.1, 0.2)	0.061	0.038	0.070
Case 9: Hub including CHP and TSS (γ=0, 0.1, 0.2)	0.061	0.014	0.074
Case 10: Hub including CHP, boiler and TSS (γ=0, 0.1, 0.2)	0.061	0.038	0.068
Case 11: Proposed strategy (γ=0)	0.057	0.038	0.068
Case 11: Proposed strategy (γ=0.1)	0.059	0.038	0.068
Case 11: Proposed strategy (γ=0.2)	0.061	0.038	0.068

In addition, the daily power curve of hubs in DA energy and reserve markets is plotted in Figs. 4 and 5. Figure 4 shows the daily curve of hubs active, reactive, gas and heating power in the DA energy market. According to Fig. 4(a), the hub’s active power is negative in the period 1:00–8:00, thus, they are consumers due to charging all electrical storages in these hours arising from low energy price based on section 4.1. However, all hubs can sell energy at other hours to the DA energy market, because:

CHPs can inject high active power into the electrical network because electricity energy price is greater than gas energy price at these hours, i.e. 9:00-24:00, based on Section 4.1.

RESs inject high power into the electrical network especially during period 9:00-20:00 because RES with type of photovoltaic and wind system generates high active power in this period in comparison with other hours based on [31].

All electrical storages such as ESS and EVs parking lot can be discharged during the peak electrical load, 17:00-22:00, where this period includes high electricity energy price with respect to other hours according to Section 4.1.

Fig. 3

Daily curve of a) heating reserve requirements, b) electrical reserve requirement [30].

Fig. 4

Expected daily hubs curve in the a) active power, b) reactive power and, c) gas/heating power based on different uncertainty levels of EVs parameter (γ).

Figure 4(b) shows the daily reactive power curve of hubs so that the hubs inject high reactive power into the electrical network at period 1:00-8:00. Because, ESSs and EVs are charged at these hours based on Fig. 4(a), thus, the electrical energy demand will be high at this period. Hence, it is expected that the value of the voltage drop is high if the reactive power management is not used in this condition. Therefore, to cope with this issue, CHPs, ESSs and EVs will inject high reactive power into the proposed network to regulate bus voltage. However, all electrical sources and storages are trying to inject high active power to increase their revenue at other hours. Hence, a high percentage of their capacity is spent to generate active power and a low percentage is allocated to generate reactive power, thus, the hubs reactive power has a low value at period 9:00-24:00 compared to 1:00-8:00. Finally, the daily hub’s heating and gas power curve in the DA energy market is presented in Fig. 4(c). Accordingly, hubs include low/high heating power at hours 1:00-4:00 and 16:00-24:00/5:00-15:00, because:

CHPs and boilers inject high heating power into the district heating grid at all simulation hours because the heating energy price is greater than the gas energy price based on Section 4.1 at this period.

TSSs are charged/discharged at period 1:00-4:00 and 16:00-24:00/5:00-15:00 due to the low/high heating energy price at this period to increase the hubs revenue.

Based on Fig. 4(c), hubs gas power is close to 3.8 p.u. at all simulation time, where this statement demonstrates the hub’s capability to obtain a flat gas power profile for gas stations. Besides, the impacts of the increasing uncertainty level of EVs parameter are investigated in Fig. 4, so that increasing γ will e reduce the hubs generation active and reactive power based on Figs. 4(a) and 4(b) due to increasing/decreasing EVs energy demand/EVs generation capacity ( ${\bar{P}}^{EV, D}$ , ${\bar{S}}^{EV}$ ) in the worst-case scenario according to the formulation (5). However, hub’s heating and gas power is independent of γ based on Fig. 4(c) because EVs operate in the electrical network.

In Fig. 5, the hubs up and down reserve power in the electrical and heating reserve markets are presented. Accordingly, all hubs can provide reserve requirements in electrical and district heating networks because the daily curves in Fig. 5 are upper than the daily curves in Fig. 2. Moreover, the daily hubs reserve power curve for different values of γ is the same in the district heating system based on Fig. 5(a) because it is independent of EVs uncertainty. Nonetheless, increasing γ causes daily reserve power curve of hubs in the electrical network to be shifted down in comparison with γ=0 according to Fig. 5(b). Because the EVs’ capacity to generate reserve power is reduced in this condition due to increasing EVs demand. Finally, the daily curve of total hubs revenue in the DA energy and reserve markets for different uncertainty levels of EVs parameters is based on Fig. 6. Note that hubs revenue at each hour depends on the different power of the hub and market price according to the first part of the equation (1). Since the hourly changing value of hubs active power is high with respect to hubs reactive, heating and gas power based on Figs. 4 and 5, and also, the electrical energy price is high in comparison with other market prices at more simulation hours according to section 4.1, the format of the daily hubs revenue curve is close to daily hubs active power curve format. Moreover, all hubs include cost (hubs revenue is negative) at period 1:00-8:00 due to charging of ESSs and EVs; however, hubs can achieve revenue from selling energy and reserve power in the DA energy and reserve markets.

Fig. 5

Expected daily hubs up and down reserve power curve in a) electrical, and b) heating networks based on different uncertainty levels of EVs parameter (γ).

Fig. 6

Expected daily hubs revenue curve based on different uncertainty levels of EVs parameter (γ).

B) Hub’s impacts on different energy networks: Table 3 presents the impacts of hubs based on the CEM on the electricity, gas and district heating networks in 11 cases. Cases 1 and 2 refer to power flow analysis in these networks with and without RESs. The other cases consider the problem (5) while hub includes various resources and storages in each case as reported in Table 3. According to this table, the electrical energy loss (EEL) and maximum voltage deviation (MVD) include a high value in Case 2, where the value of these indices is more than Case 1 that considers power flow analysis. This statement demonstrates the RESs in Case 2 inject high active power into the electrical network so that the values of EEL and MVD are increased in comparison with Case 1. However, the values of these indices can be reduced in cases that contain flexible sources such as EVs parking lot (Case 3), ESS (Case 4), CHP (Case 5) or all these sources (cases 6 and 11). The best condition from the viewpoint of these indices occurs in cases 6 and 11 due to a minimum value of EEL and MDV compared to other cases. Regarding the district heating network indices, the heating energy loss (HEL) and maximum temperature deviation (MTD) contain high values in cases 1-4 that do not include local sources to generate heating power. However, these indices will be improved if these sources in the hub framework are added to the district heating grid based on the results of cases 5-11 in Table 3. Moreover, the gas system indices such as gas energy loss (GEL) and maximum pressure deviation (MPD) are zero in cases 1-4 due to zero gas demand. Despite that, by adding a hub with a CHP and a boiler that are gas consumers into the natural gas network in cases 5-11, the GEL and MPD will be increased with respect to the results of cases 1-4. In general, it is seen that Case 11 (the proposed strategy) has the best condition from the viewpoint of system operation indices, so that it can improve EEL, MVD, HEL and MTD about to 45.5% ((5.4102 –2.9511)/5.4102), 52.3% ((0.128 –0.061)/0.128), 37.7% ((2.1156 –1.3176)/2.1156) and 38.7% ((0.111 –0.068)/0.111), respectively. These benefits will be obtained while GEL and MPD are increased from zero to 2.9831 p.u. and 0. 038 p.u. The impacts of the uncertainty level of EVs parameter (γ) create only in cases 3, 6 and 11 that are contained EVs parking lot. Accordingly, increasing γ causes increases EEL and GEL as well as MVD in comparison with γ=0. Because, EVs demand will be increased if γ is a high value based on problem (5); hence, increasing EVs demand results in increasing EEL and MVD, and GEL due to increasing CHP generation active power in this condition.

C) Investigation of the flexibility of the energy hub: Fig. 7 shows the curves of the COLF and hub revenue in VOLF. Accordingly, the high revenue is obtained for all hubs from DA energy and reserve markets in the case of VOLF = 0. In this condition, the COLF will be a high value with respect to the other value of VOLF according to Fig. 7(a). In other words, the high flexibility (COLF = 0) will occur for hubs if VOLF is increased to 10⁶ $/MWh; hence, the total hubs revenue reduces about to 24.3% in VOLF = 10⁶ compared to VOLF = 0 according to Fig. 7(b). Note that the hubs revenue is decreased in high flexibility condition because the operation cost of flexibility sources such as ESS, TSS, CHP, boiler, and EVs parking lot is increased in these cases. In addition, increasing γ is related to increasing COLF and decreasing hubs revenue as well as increasing the optimal value of VOLF that is suitable to obtain zero COLF based on Fig. 7.

Fig. 7

Expected COLF and revenue curves in VOLF based on different uncertainty levels of EVs parameter (γ).

5 Conclusions

In this paper, the CEM strategy is implemented on different sources and storages in a flexible grid-connected energy hub framework, where the proposed hub can participate in DA energy and reserve markets to achieve high profit to these sources and storages. Hence, the proposed problem maximizes the difference between the expected hubs revenue in the proposed markets and the expected COLF while it is subject to LOPF formulation in different energy networks, and the hub model contains the renewable and flexible sources constraints and the hub reserve equations. In the next step, the HSRO method based on a hybrid approach of BURO and SBSP is used to model uncertainties of the load, reserve demand, market price, RES generation power and EVs demand. Finally, according to the numerical results obtained for different case studies, the use of CEM strategy for renewable and flexible sources in the hub scheme can increase hubs revenue by approximately 18.1% and 22.3% in DA energy and reserve markets concerning UEM strategy. Also, the presence of the hub in different energy networks improves EEL, MVD, HEL and MTD respectively about to 45.5%, 52.3%, 37.7% and 38.7% in comparison with power flow analysis in these networks, while GEL and MPD in the proposed case are increased compared to the power flow analysis in the natural gas network. Moreover, the high flexibility (COLF = 0) hub can be achieved for a high value of VOLF with reducing hub revenue about to 24.3% in this condition. Therefore, the benefits of the proposed strategy are: 1) obtaining a hub with high flexibility, 2) achieving a suitable local source as a hub in different energy networks to improve indices of these grids, and 3) improving the revenue of renewable and flexible sources in DA energy and reserve markets using a CEM strategy in the hub.

References

Thillainathan

and Dipti

, Multi-agent system for managingdistributed energy storage and electrical vehicles, IntelligentDecision Technologies 9 (2015), 181–190.

Ghorbani

, Afsharchi

and Derhami

, A multi-agentreinforcement learning algorithm with fuzzy approximation fordistributed stochastic unit commitment, Journal of Intelligent & Fuzzy Systems 37 (2019), 6613–6628.

Kim

, Cho

and Shin

, Advanced power distribution systemconfiguration for smart grid, IEEE Transactions on Smart Grid 4 (2013), 353–358.

Liu

, Zhang

, Dai

and Wu

, Intelligent modeling andoptimization for smart energy hub, IEEE Transactions onIndustrial Electronics 66 (2019), 9898–9908.

Mansouri

S.A.

, Ahmarinejad

, Ansarian

, Javadi

M.S.

and Catalao

J.P.S.

, Stochastic planning and operation of energy hubsconsidering demand response programs using Benders decompositionapproach, International Journal of Electrical Power & Energy Systems 120 (2020), 106030.

Gholinejad

H.R.

, Loni

, Adabi

and Marzband

, A hierarchicalenergy management system for multiple home energy hubs inneighborhood grids, Journal of Building Engineering 28(2020), 101028.

Guen

M.L.

, Mosca

, Perera

A.T.D.

, Coccolo

, Mohajeri

and Scartezzini

J.L.

, Achieving energy sustainability in future neighborhoods through building refurbishment and energy hub concept: a case study in Hemberg-Switzerland, Energy Procedia 122 (2017), 265–270.

Bostan

, Setayesh-Nazar

, Shafie-khah

and CatalãoOptimal

J.P.S.

, scheduling of distribution systems considering multipledownward energy hubs and demand response programs, Energy 190 (2020), 116349.

Nikmehr

, Distributed robust operational optimization of networked microgrids embedded interconnected energy hubs, Energy 199 (2020), 117440.

10.

Mirzapour-Kamanaj

, Majidi

, Zare

and Kazemzadeh

, Optimal strategic coordination of distribution networks and interconnected energy hubs: A linear multi-follower bi-level optimization model, International Journal of Electrical Power & EnergySystems 119 (2020), 105925.

11.

Wang

, Liu

and Liu

, Coordinating energy management for multiple energy hubs: From a transaction perspective, International Journal of Electrical Power & Energy Systems 121 (2020), 106060.

12.

Heidari

, Mortazavi

S.S.

and Bansal

R.C.

, Stochastic effects of ice storage on improvement of an energy hub optimal operation including demand response and renewable energies, Applied Energy 261 (2020), 114393.

13.

Heidari

, Mortazavi

S.S.

and Bansal

R.C.

, Equilibrium state of aprice-maker energy hub in a competitive market with priceuncertainties, IET Renewable Power Generation 14(2020), 976–985.

14.

, Liu

, Zhang

and Wang

, Real-time rolling horizonenergy management for the energy-hub-coordinated prosumer communityfrom a cooperative perspective, IEEE Transactions on PowerSystems 34 (2019), 1227–1242.

15.

Jadidbonab

, Dolatabadi

, Mohammadi-Ivatloo

, Abapour

and Asadi

, Risk-constrained energy management of PV integrated smartenergy hub in the presence of demand response program and compressedair energy storage, IET Renewable Power Generation 13(2019), 998–1008.

16.

, Wei

, Mei

, Hu

and Wu

, Participation of an energyhub in electricity and heat distribution markets: An MPEC approach, IEEE Transactions on Smart Grid 10 (2019), 3641–3653.

17.

Bahrami

, Toulabi

, Ranjbar

, Moeini-Aghtaie

and Ranjbar

A. M.

, A decentralized energy management framework for energy hubsin dynamic pricing markets, IEEE Transactions on Smart Grid 9 (2018), 6780–6792.

18.

Dini

, Pirouzi

, Norouzi

M.A.

and Lehtonen

, Grid-connected energy hubs in the coordinated multi-energy management based on day-ahead market framework, Energy 188 (2019), 1–12.

19.

Flexibility and Aggregation Requirements for their interaction in the market, Eurelectric, 2014.

20.

Bozorgavari

S.A.

, Aghaei

, Pirouzi

, Nikoobakht

, Farahmand

and Korpås

, Robust planning of distributed battery energy storage systems in flexible smart distribution networks: A comprehensive study, Renewable and Sustainable Energy Reviews 123 (2020), 109739.

21.

Norouzi

M.A.

, Aghaei

and Pirouzi

, Enhancing distribution network indices using electric spring under renewable generation permission, International Conference on Smart Energy Systems and Technologies (SEST) 2019, pp. 1–6.

22.

Shahbazi

, et al., Hybrid stochastic/robust optimization model for resilient architecture of distribution networks against extreme weather conditions, International Journal of Electrical Power & Energy Systems 126 (2021), 106576.

23.

Shahbazi

, et al., Effects of resilience-oriented design on distribution networks operation planning, Electric Power Systems Research 191 (2021), 106902.

24.

Kim

J.S.

and Edgar

T.F.

, Optimal scheduling of combined heat and power plants using mixed-integer nonlinear programming, Energy 77 (2014), 675–690.

25.

Pirouzi

, Latify

M.A.

and Yousefi

G.R.

, Conjugate active and reactive power management in a smart distribution network through electric vehicles: A mixed integer-linear programming model, Sustainable Energy, Grids and Networks 22 (2020), 100344.

26.

Pirouzi

, et al., Robust linear architecture foractive/reactive power scheduling of EV integrated smartdistribution networks, Electric Power Systems Research 155 (2018), 8–20.

27.

Pirouzi

, et al., Two alternative robust optimization models forflexible power management of electric vehicles in distributionnetworks, Energy 141 (2017), 635–651.

28.

Bozorgavari

S.A.

, Aghaei

, Pirouzi

, Vahidinasab

, Farahmand

and Korpås

, Two-stage hybrid stochastic/robust optimal coordination of distributed battery storage planning and flexible energy management in smart distribution network, Journal of Energy Storage 26 (2019), 100970.

29.

Norouzi

M.A.

, Aghaei

, Pirouzi

, Niknam

and Lehtonen

, Flexible operation of grid-connected microgrid using ES, IET Generation, Transmission & Distribution 14 (2020), 254–264.

30.

Nord Pool, https://www.nordpoolgroup.com/historical-market-data/

31.

Abrisham Foroushan Asl

, Gandomkar

and Nikoukar

, Optimal protection coordination in the micro-grid including inverter-based distributed generations and energy storage system with considering grid-connected and islanded modes, Electric Power Systems Research 184 (2020), 106317.

32.

Generalized Algebraic Modeling Systems (GAMS). [Online]. Available: http://www.gams.com

γ		0		0.1		0.2
Energy management method		DA market		DA market		DA market
		Energy	Reserve	Energy	Reserve	Energy	Reserve
UEM	CHP	79.8	24.7	79.8	24.7	79.8	24.7
	Boiler	40.3	14.2	40.3	14.2	40.3	14.2
	TSS	21.2	10.3	21.2	10.3	21.2	10.3
	RES	57.5	16.4	57.5	16.4	57.5	16.4
	ESS	46.7	32.7	46.7	32.7	46.7	32.7
	EVs	–51.3	9.3	–58.4	7.5	–64.1	5.7
	Sum	194.2	107.6	187.1	105.8	181.4	104
CEM	229.35	131.6	218.72	128.44	209.88	125.32
Gain rate (%)	18.1	22.3	16.9	21.4	15.7	20.5