A hybrid robust‐stochastic approach for optimal scheduling of interconnected hydrogen‐based energy hubs

Danida Fellowship Centre, Grant/Award Number: 18‐M06‐AAU Abstract The energy hub (EH) concept is an efficient way to integrate various energy carriers. In addition, demand response programmes (DRPs) are complementary to improving an EH's efficiency and increase energy system flexibility. The hydrogen storage system, as a green energy carrier, has an essential role in balancing supply and demand precisely, similar to other storage systems. A hybrid robust‐stochastic approach is applied herein to address fluctuations in wind power generation, multiple demands, and electricity market price in a hydrogen‐based smart micro‐energy hub (SMEH) with multi‐energy storage systems. The proposed hybrid approach enables the operator to manage the existing uncertainties with more flexibility. Also, flexible electrical and thermal demands under an integrated demand response programme (IDRP) are implemented in the proposed SMEH. The optimal scheduling of the hydrogen‐based SMEH problem considering wind power generation and electricity market price fluctuations, as well as IDRP, is modelled via a mixed‐integer linear programming problem. Finally, the validity and applicability of the proposed model are verified through simulation and numerical results.


| INTRODUCTION
In conventional power systems, various energy carriers, for example, natural gas and electricity, were designed and operated separately. Therefore, natural gas energy is the input energy for conventional electric power units in order to supply only electrical demands, which is defined as a separately operated approach. However, current efforts aim to provide an infrastructure to design and operate multi-carrier energy together. In addition, the energy management of multiple energies has a significant role in research papers [1][2][3][4]. The synergy of various energy carriers under the energy hub (EH) concept is an efficient way to achieve the goal of integrating energy carriers. EH is an infrastructure that can help the power system operate more efficiently, cost-effectively, and reliably with lower pollutant emissions [5][6][7]. Furthermore, demand response programs (DRPs) can be applied to help load-shaping goals defined in demand-side management (DSM) strategies [8]. In this way, consumers' load patterns can be controlled, and consequently, the efficiency of EH and flexibility of the power system will be increased [9]. On the other hand, global concerns about environmental emission have compelled researchers to use green energy carriers such as hydrogen. Hydrogen gas is a clean, highly plentiful, and non-poisonous renewable fuel. By developing technology, electric power can be converted to hydrogen (P2H) in the low electricity price hours through electrolysation of water into oxygen and hydrogen and subsequently stored in a hydrogen storage tank to be used in industry, hydrogen-based vehicles, or converted to electricity (H2P) in the high electric price hours [10].
The first category in the literature is related to the various emerging technologies and their impacts on the EH systems. In [11], improving both the supplier-side and customer-side benefits by introducing the integrated demand response programme (IDR)-based non-cooperative game for both electricity and natural gas demands. The other management game approach for the EH is proposed in [10] to decrease the smart energy hub (S.E. Hub) energy cost and also the peak-toaverage ratio electrical load. The presented game in [12] is integrated demand-side management (IDSM) to handle the interaction among S.E. Hubs via the utilisation of a cloud computing structure. Various distributed energy sources such as CHP, wind turbine, energy storages, and DRP are applied in [13] to reduce loss factor and load factor by the presented new control strategy. Optimal daily scheduling management of EH integrated with the hydrogen system and DRP in [14] simultaneously minimises operating costs and emissions. The authors in [15] recommended optimal electric and heat management for the residential EH besides taking into account the inclusive DRP to verify the significant effects on the energy cost decrement. The other impacts of taking modern technologies into the EH are relevant to the short-term scheduling framework, which has been represented in [16]. Also, in [16], sensitivity analysis is rendered in determining the outcomes of considering smart grid appliances and input EH parameters on the expenditures. Some literature has been dedicated to different systems of uncertainty management based on risk-based approaches in addition to the presented works. For instance, in [17], a risk valuation structure with both deterministic and stochastic procedures for multi-objective decisions is taken into account for developing sustainable manufacturers' engineering projects. The authors in [18,19] have introduced the proper risk assessment methodology, that is, conditional value-at-risk (CVaR) for wind and PV-integrated smart multiple-carrier EH (SMCEH) in the presence of compressed air energy storage (CAES) and DRPs. Furthermore, in [20], the valuable outcomes and effects of applying CAES and P2G facilities to reduce the proposed EH system's operating cost besides considering various uncertain parameters of the system with utilising the CVaR risk measurement technique is validated. The EH management-based Information Gap Decision Theory (IGDT) model, along with consideration of DRP, has been accomplished in [21] to determine the robust decisions against the uncertainties of electricity price and wind power generation. A robust scheduling methodology is employed for the micro-EH in [22], which is integrated with technologies as an integrated demand response programme (IDR) and hydrogen storage system (HSS) to minimise total operating expenditures against variations in the electricity market. Another strong robustness-based game theory approach to determining the optimum scheduling problem of a multiple-EH system (MEHS) has been presented in [23]. Moreover, some literature has applied hybrid approaches to managing different uncertainties simultaneously. A hybrid stochastic/IGDT approach has been applied in [24] to consider wind generation and energy prices uncertainty. A hybrid interval/IGDT/stochastic model has been proposed in [25] to manage the uncertainties of the EH. To alleviate renewable uncertainties at the hub level, a hybrid box-polyhedral uncertainty set-based robust optimisation method has been proposed in [26] and also, power to gas technology along with a crossvector DR has been proposed that enables end-users to shift between electricity and natural gas. However, power to hydrogen technology, which has many economic benefits, has not been considered. Also, electricity price and load uncertainties have been ignored previously and these drawbacks are addressed herein.
According to the aforementioned works and to the best of the authors' knowledge, there is no focus on applying a hybrid stochastic/robust method in order to manage the different uncertainties in the EH. The proposed hybrid stochastic/robust procedure has been utilised in energy systems, for example, for scheduling the microgrid [27], etc. However, there is a research gap in utilising a hybrid stochastic/robust approach for multicarrier systems. Moreover, most articles have ignored utilising the power-to-hydrogen technology in the EH, which can bring many benefits to the EH operator. Herein, the day-ahead scheduling of hydrogen-based smart micro-energy hub (SMEH) is emphasised, as depicted in Figure 1, considering an incentive-based integrated demand response programme (IDRP). Using IDRP, both electrical and thermal loads of consumers can be controlled and shifted from on-peak periods to off-peak periods. HSS has an essential role in precise supplydemand balancing, similar to other storage systems. By applying HSS, the electrical power from a wind turbine unit can be converted to hydrogen at low electricity price periods and then reconverted to electricity at high electricity price periods. Finally, a hybrid stochastic/robust method is applied to manage fluctuations in wind power, multiple demands, and electricity market prices. The utilised hybrid approach enables the EH operator to make decisions with greater flexibility. Indeed, the EH operator can choose how much risk to take by adjusting the uncertainty budget parameter. Table 1 reviews the contribution of this study and compares it with previous works. In short, the main contributions are as follows: -Focussing on power-to-hydrogen technology and investigating its economic benefits. In addition to hydrogen storage, impact of multiple energy storage systems including electrical, gas, and thermal storage systems are evaluated.
-Various uncertain parameters including wind power, electrical and heat demands, and electricity market prices are considered and alleviated through a hybrid robust/stochastic approach. -To apply the demand-side management, an incentive-based demand response programme is applied for both electrical and thermal demands. Therefore, multi-carrier system endusers shift their loads and reduce the EH costs.
A brief description of HSS technology is given in Section 2. The problem formulation of the objective function, all constraints related to the technologies including HSS and IDRP, and also scenario-based stochastic objective function are introduced in Section 3. Simulation and analysis of numerical results are provided in Section 4. Finally, Section 5 concludes the paper.

| NOVEL HSS TECHNOLOGY
Metal hydride-based HSS is a modern technology, which has various characteristics. Hydrogen adsorption and desorption occur at a fixed amount of pressure; the high-and lowtemperature metal hydrides have analogous self-regulation; thus, solid hydrogen storage, compared to liquid and gas storage forms, is able to accumulate more hydrogen [28]. As a result, the selection of hydrogen-storing substances will specify the storage capacity, and the form of storage will affect the allocation. The generated hydrogen might be transmitted to the methane distribution network or the hydrogen market.
Reconverting to electricity through the utilisation of a fuel cell (FC) is another option for use of the produced hydrogen. It is expected that the hydrogen-based economy will be extended significantly, replacing the petroleum-based economy over the next few years [28]. The optimistic outlook for using hydrogen commercially, predicts that millions of businesses will be constructed worldwide, such as manufacturing, theoretical scientific research and development (R&D), and sales. The conceptual structure of HSS technology is presented in Figure 2.

| Wind generation model
Wind generation mainly depends on wind speed. Equation (1) shows the relationship between wind turbine output power and wind speed.
where P W T ;Max is the maximum power of the wind turbine; V W T t;s is the wind speed at time step t and scenario s; V CI ; V CO are the cut-in and cut-off speeds; and V R is the rated speed. Therefore, the uncertain parameter is wind speed, which is unpredictable. Some methods have been proposed to model the wind speed in the literature, for example, a data-driven model [29] and clustering model [30]. However, generally, Weibull distribution has been used to model the intermittency of wind speed. The Weibull probability density function (PDF) is given in (2), and due to its adaptability, it is also used herein to produce different scenarios.
In Equation (2),k is a shape factor and c is a scale factor.

| Multi-demand model
In order to model the loads' uncertainty, normal distribution with a mean of μ d and standard deviation of σ d is utilised to produce scenarios as follows:

| Stochastic approach
As already outlined, the hydrogen-based SMEH problem is a minimisation problem in which the objective function is related to the total operating costs. In the objective function (4), various terms are considered as follows: where t and s are indices of hours and scenarios, respectively, while N T , N s are the number of hours and scenarios; π s is the probability of each scenario; λ E t;s and λ G t are the electricity price and the natural gas price, respectively; . The objective function must also satisfy all the following equality and nonequality constraints.
The generation of power and heat are dependent on each other in the CHP unit, as presented in Equations (5) (6) (7). Indeed, these relations indicate an operation region for the CHP unit, namely a feasible operation region (FOR), that is shown in Figure 3. Equations (8) and (9) show the limitations of generating power and heat via the CHP unit, respectively. The maximum electrical and thermal power productions of units are stated as P CHP A and H CHP B , respectively. In the above equations, the parameter is a significant positive number. Also, in Equation (10), expressing the limitation of electric power generated by a gas turbine (GT) unit is indicated. In addition, the binary decision variable of the CHP unit for both electric and heat production at time t and scenario s is denoted by are the four-point heat power produced correlated with electric power in each vertex; P GT t;s is the power generated by the gas turbine unit at time t and scenario s; P GT ;max is the maximum level of generated power by the gas turbine unit; and I GT t;s is the binary variable decision of gas turbine unit. Ramp-up and ramp-down constraints of gas-based power units, that is, CHP and GT, are presented in Equations (11) and (12). This means that the generated power should be lower than the ramp rate capacities of gas-based units for the next scheduling hour. The binary variables, that is, Y i;t;s and Z i;t;s are defined as the start-up and shut-down indicators obtained based on the unit commitment (UC) programme of units, which corresponds with Equations (13) and (14).
where i is the index for thermal units; P i;t;s is the electric power produced by the gas-based power units, that is, CHP and GT; P min i is the minimum level of generated electricity for these two gas-based power units; I i;t;s is the binary decision variable of the mentioned two units similar to the definition of the CHP unit; Y i;t;s and Z i;t;s are also binary variables to show start-up and shutdown states of the thermal units; and R up i and R dn i are the up and down ramp rates of two gas-based power units.
Also, for gas-based power units, linearised minimum uptime and down-time relations are considered as Equations (15) (16) (17) (18) and (19) (20) (21) (22), respectively. These limitations determine that the starting up and shutting F I G U R E 3 Gas and electricity price profiles MANSOUR-SATLOO ET AL.
-245 down of the gas-based units must be done after passing their relative minimum up-time and down-time. Constraints (23)- (27) are utilised to determine how much natural gas is consumed by gas-based power units. The start-up and shut-down natural gas consumption are located in Equations (23) (24) (25) (26) by SU i,t,s and SD i,t,s , respectively. Total gas consumed, that is, GC i,t,s , is denoted by Equation (27), in which the heat rate value (HRV) parameter is the conversion coefficient of the electric unit to the gas unit. The efficiency of power produced from each unit affects total gas consumption, which is indicated in Equation (27).
where sug i is the required natural gas at the start-up and shutdown period of unit i while, SU i,t,s and SD i,t,s are the natural gas consumption in the start-up and shut-down hours of unit i at time t and scenario s; GC i,t,s is the total gas consumption of unit i at time t and scenario s; η CHP and η GT are the efficiencies of CHP unit and gas turbine unit, respectively; and HRV i is the heat rate value of thermal units.
Heat power produced by the boiler unit is stated in Equation (28) denoted by H B t;s , which is limited through the maximum and minimum heat power production. Also, the boiler unit's gas consumption, that is, GB B t;s is defined as in Equation (29)  In the following constraints, electrical energy storage (ESS) must be charged in low-price hours while being operated in discharge mode in high-price hours to reduce the total imported electric power from the main grid. Therefore, constraints (30) and (31) illustrate the charging and discharging electric power limited by their maximum and minimum charging and discharging power values. This action must not occur at the same time, as illustrated in Equation (32) are the binary decision variables for charging and discharging power by TSS unit at time t and scenario s; η chr;T S and η dischr;T S are the efficiencies of charging and discharging heat power modes; E min, TS and E max, TS are the minimum and maximum thermal energy stored level; η TS is the energy dissipation efficiency; and E T S t;s is the level of stored thermal energy in TSS unit at time t and scenario s.
Also, in order to consider the gas storage system (GSS), constraints (42)-(47) are outstanding and determine the charging and discharging power of GSS for each hour and scenario [22,31]. are the charging and discharging gas power of GSS unit at time t and scenario s; G min,chr,GS and G max,chr,GS are the minimum and maximum state of charging power and also G min, dischr,GS and G max, dischr,GS are the minimum and maximum state of discharging gas power by GSS unit; I chr;GS t;s and I dischr;GS t;s are the binary decision variables for charging and discharging power by GSS unit at time t and scenario s; η chr, GS and η dischr,GS are the efficiencies of charging and discharging gas power modes; E min, GS and E max, GS are the minimum and maximum gas energy stored level; Δt is the interval time for calculating gas energy level of GSS unit; E GS t;s is the level of stored gas energy in GSS unit at time t and scenario s.
In HSS technology, the converted electricity to hydrogen, denoted as P P2H t;s , is the charging equivalent power while the hydrogen converted to electricity, that is, P H2P t;s , is the discharging equivalent electric power, which is restricted by Equations (48) and (49), respectively. Also, the other constraints (i.e., (50)-(53)) are precisely the same as the previously defined energy storage systems [10,22]. It should be noted that H2P conversion via a fuel cell is a clean process with zero carbon emission; however, supplying hydrogen to the gas network and then to the CHP unit to generate power causes carbon emission, which has environmental impacts. Also, supplying hydrogen to the gas network is possible up to 2%-7%, which presents the P2H conversion with a limitation [32]. On the other hand, if wind energy converts to hydrogen then to methane, the system efficiency decreases. According to these statements, P2H and H2P via a fuel cell have been proposed as a flexible and clean solution.
where P P2H t;s and P H2P t;s are the charging and discharging electric power of HSS unit at time t and scenario s; P min,P2H and P max,P2H are the minimum and maximum state of charging power and also P min,H2P and P max,H2P are the minimum and maximum state of discharging electric power by HSS unit; I P2H t;s and I H2P t;s are the binary decision variables for charging and discharging power by HSS unit at time t and scenario s; η P2H and η H2P are the efficiencies of charging and discharging electric power modes; E min, hyd and E max, hyd are the minimum and maximum energy stored level; and E hyd t;s is the level of stored energy in HSS unit at time t and scenario s.
-247 Therefore, Equations (54) and (55) introduce a specified amount of electrical demand, which can be curtailed from peak hours and shifted into the off-peak hours. Likewise, for the thermal demand, the acceptable range of curtailing and shifting load are determined in Equations (56) and (57). Predefined parameters DRE and DRT are considered as the maximum allowed participation factors for electrical and thermal demands, respectively. Therefore, the amount of transferred loads to off-peak times should be equal to the curtailed loads in peak times, which is expressed by Equations (58)

| Hybrid stochastic/robust approach
In the previous section, the probabilistic formulation of the EH was introduced, in which the wind power uncertainty was alleviated; however, probabilistic-based methods are not powerful enough to tackle the electricity market uncertainties. Therefore, in this section, a hybrid stochastic/robust methodology is presented in order to deal with electricity price uncertainty. The EH operator can make a decision to be optimistic, deterministic, or pessimistic using an integer variable Γ s into the interval of ½0 : N k �. The robust part of the objective function is as follows: The above equation is a min-max model composed of two parts. The outer part minimises the EH's operation cost and the inner part deals with the electricity market uncertainty. The inner part of Equation (65), that is, the max term, can be stated as follows: To convert the min-max problem into a minimisation problem, the strong duality theory can be used.
where α s and β t;s are dual variables; y t;s is an auxiliary variable. Therefore, a hybrid stochastic/robust objective function is as follows:

| SIMULATION RESULTS
The proposed stochastic mixed-integer linear programming (MILP) methodology is solved through the CPLEX 12.9 solver in GAMS 27.3 software on a PC with an Intel Core i7 processor with 2.50 GHz and 8 GB of RAM. The Monte-Carlo simulation (MCS) approach is employed herein to generate 1000 scenarios to evaluate multiple demands and wind power generation with normal PDF and Weibull PDF, respectively. To decrease the number of scenarios, the fast-forward selection algorithm is used in Matlab R2019b; therefore, 10 scenarios with the incidence possibility of each scenario are produced.
Fast forward selection is an algorithm based on probability distance, which is utilised to reduce the number of scenarios in a way that reduced scenarios can represent the key feature of the original scenarios [33]. Finally, the reduced scenarios are illustrated in Table 2. It should be mentioned that simulation results are demonstrated only for specific scenario number 3. The electricity and natural gas hourly prices are depicted in Figure 3; also, wind power production and gas demand are illustrated in Figure 4. Figures 5 and 6 specify the electric power and heat power dispatches during scheduling time. The power dispatches of CHP and GT units are demonstrated in Figure 5 in which, due to low-electricity prices between hours 1-4, the CHP unit does not prefer to commit on. In two hours (i.e., hours 5 and 6), the CHP unit is in the maximum electric power mode and then its heat production is equal to zero. In the remaining hours of Figure 5, besides the CHP generated power, the GT is only committed in peak times (i.e., hours [16][17][18][19][20][21][22] to meet the electrical demand. By explanations for Figure 5, the heat power generated by CHP and gas boiler units are illustrated in Figure 6. According to Figure 6, the gas boiler contributes to supplying of thermal demand during the scheduling horizon and performs a considerable role, and the CHP unit is online after t=7. Electrical and thermal demands are presented before and after applying the proposed IDRP in Figures 7 and 8, respectively. As shown in Figure 7, the electrical demand is curtailed from peak hours (i.e., hours [12][13][14][15][16][17][18][19][20][21][22] and shifted to off-peak hours (i.e., hours 1-10 and 22-24). By referring to Figure 8, it can be seen that the thermal demand has the same behaviour in which the curtailed load has occurred in peak periods (i.e., hours [13][14][15][16][17][18][19][20][21][22] and increased demand takes place in off-peak scheduling periods (i.e., hours 1-12 and 23-24). The total gas and electricity purchased from the main upstream grid are reduced by these descriptions compared to the system which does not have any storages or any utilised IDRPs. Therefore, the obtained related result is shown in Figure 9. As is clear from Figure 9, in off-peak times (i.e., hours 5-7), the electricity price is low. Hence, the EH operator prefers to buy the required electric power instead of committing the units, which have a high-cost coefficient in power production. However, in peak hours (i.e., hours [19][20] and the periods with additional power, the rest of the electric power is sold to the main grid to minimise operating costs. Likewise, due to gas consumption in gas-based power units (i.e., CHP, gas turbine, and boiler), the imported gas from the main gas network is varied along the scheduling period.
The SoC of ESS, TSS, and GSS technologies are presented in Figure 10, which illustrates that the storages are charged and discharged in low-and high-price hours, respectively.

F I G U R E 6
Heat power produced by CHP and gas boiler units F I G U R E 7 Electrical demand before and after applying IDRP at scenario number 3 As depicted in Figure 10, TSS is discharged earlier than the discharging of GSS technology; however, ESS discharged electric power only in peak hours, which is later than the discharge starting time of the two other energy storages. It should be mentioned that the energy stage of all storages is equal in the first and final hours, as is evident in Figure 10. Figure 11 depicts that HSS technology must be charged in low-price hours to energise the SMEH system later in the day, which reduces the total operation cost and power purchased cost from the main grid. In addition, the SoC of HSS is presented in Figure 11, which verifies that the charging and discharging of power have taken place.
A comparison of all cases is provided in Table 3. As reported in this table, in the presence of all considered energy storage facilities, but without any utilised DRPs, the total operation cost is only decreased by 2%; however, employing IDRP results in a high reduction of 11%. The share of ESS, TSS, GSS, and HSS in reducing the total operation cost are $7.35, $4.01, $2.04, and $6.81, respectively. Furthermore, without utilised any DRPs, the reduction of exchanged power F I G U R E 8 Thermal demand before and after applying integrated demand response programme at scenario number 3 F I G U R E 9 Imported gas and exchanged electric power at scenario number 3 F I G U R E 1 0 State of charge of electrical energy storage, thermal storage system and gas storage system technologies cost with the main grid is about 2.5%; however, applying IDRP leads to a reduction of approximately 15.6%. The cost of natural gas purchased from the main grid is almost the same as the case of any energy storage, and DRPs are considered.
This part of the results section presents a hybrid stochastic/robust scheduling of the EH. To this end, the optimisation problem of (75)-(81) was solved, and the results are introduced. Figure 12 indicates the exchanged power with the main grid through different uncertainty budget amounts. Whatever the uncertainty budget increases, the EH operator becomes pessimistic and tries to make a decision with more robustness; therefore, the cost of the exchanged power with the main grid increases, as shown in Table 4. In addition, Table 4 shows the uncertainty budget impacts on the total operation cost. It is obvious that in the higher robustness levels, the EH operator bears high operation costs. Finally, Table 5 analyses the price deviation and uncertainty budget impact on the expected total operation cost. As can be seen, more price deviation in the high uncertainty budgets applies at higher operation costs. Therefore, more robustness levels mean more operating costs. The optimal scheduling of the hydrogen-based SMEH problem was analysed herein by considering fluctuations of wind power production, multiple demands, and electricity price, which is modelled as a MILP problem. A hybrid robust-stochastic model was introduced to address the uncertainties of wind power, various demands, and electricity price. The proposed model enabled the operator to use the advantages of both stochastic and robust optimisation approaches to manage the uncertainties. By considering HSS technology, the total operating cost reduction and the cost of exchanged power were obtained. Therefore, the EH system, which does not have any HSS, is equipped with the proposed HSS facility and the level of decreased total operation cost and the cost of electricity purchased from the main grid are nearly 1% and 1.1%, respectively. In addition, employing IDRP in the proposed SMEH system caused shifting electrical and thermal demands from peak periods into off-peak periods, which was conducted in 11% and 15.6% decrements in the total operating expenditure and cost of the exchanged power with the main grid, respectively. The other outstanding result is relevant to employing the proposed hybrid stochastic/robust to the SMEH system. In the highest price deviation and uncertainty budget levels, the considered system is shown to be a robust approach against the environment with various high uncertainties. Finally, simulation and numerical results indicated that the proposed model is optimal, economical, and applicable.