Multiobjective Probabilistic Planning of Energy Hub With Hydrogen Storage Technologies Considering Demand Response Programs

Abstract

This paper develops a stochastic bi-objective energy management system (EMS) for an integrated energy hub (EH) comprising photovoltaic (PV) arrays, wind turbines (WTs), a dual-fuel boiler, combined heat and power (CHP) generation, electric vehicle (EV) charging infrastructure, and hydrogen storage systems, interconnected with the main grid. The proposed EMS framework simultaneously minimizes operational expenditures (OPEX) and carbon emissions while addressing uncertainties in renewable generation and load demand through probabilistic modeling and demand response programs (DRPs). A novel modified multi-objective grasshopper optimization algorithm (MMOGOA) with adaptive mutation operators is introduced to solve this complex optimization problem, demonstrating superior convergence characteristics and 7.2% lower OPEX compared to conventional MOEAs (Non-dominated Sorting Genetic Algorithm [NSGA-II] and MOPSO) in baseline scenarios. Comprehensive simulations reveal that demand response program (DRP) implementation achieves significant reductions (18.87% in costs and 14.62% in emissions), while uncertainty incorporation increases costs by 10% and emissions by 4.38%, with MMOGOA consistently maintaining performance dominance across all operational regimes. The results quantitatively highlight the importance of optimizing DRP participation and managing uncertainties to improve the efficiency and sustainability of energy management systems (EMSs).

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

1. Introduction

The energy hub (EH) stands as a pivotal concept in orchestrating the joint operation of diverse energy carrier infrastructures, notably encompassing natural gas and electricity. Its multifaceted capabilities extend to achieving a spectrum of objectives, ranging from minimizing operational costs and enhancing reliability and efficiency to curbing emissions [1–3]. The optimal functioning of multicarrier energy systems bears the potential to yield manifold benefits across technical, economic, and environmental domains. These advantages encompass heightened system reliability, reductions in operational costs, fuel consumption, and greenhouse gas emissions [4]. Over recent years, a myriad of methodologies have been proposed for the energy management system (EMS) governing EH operations.

In terms of methods and planning for the EMS in energy hubs (EHs), the literature can be categorized into two main groups of traditional methods and artificial intelligence techniques.

Traditional energy management methods include linear programing (LP) [5], nonlinear programing (NLP) [6], and rule-based dynamic programing [7]. Recent studies focus on optimizing EHs using mixed-integer linear programing (MILP). One approach combines energy trading with the national grid and battery storage, solving the problem with LP and MILP [8]. Another model uses MILP for energy consumption planning, integrating a demand response (DR) program to handle time-varying generation resources and reduce peak load [9]. A university campus EMS balances energy transactions using MILP and predicts energy output with a radial basis function neural network (RBF-NN) [8]. Battery management is highlighted as an economic solution to reduce operating costs. Other models aim to maximize energy exchange revenue through DR and use Lyapunov optimization for real-time management of thermal and electrical energy, achieving lower operational costs [10, 11].

Recent studies in EH research have focused on employing intelligent methods and advanced algorithms for optimal energy management. In [12], a reinforcement learning framework was developed for sustainable EH management in urban environments. Research [13] presented a distributed multi-objective control strategy considering energy conversion characteristics. In [14], a multi-objective model based on innovative optimization algorithms was designed for power systems and microgrids. Study [15] addressed economic and environmental aspects by proposing an integrated model for load management, pollution control, and pricing in large-scale EHs. In [16], distributed energy management solutions were developed for interconnected EHs, while research [17] employed a two-timescale reinforcement learning system to enhance operational flexibility. A deep reinforcement learning method was applied in [18] for real-time flexibility optimization in grid-connected smart building communities. Also, the metaheuristic methods, such as particle swarm optimization (PSO) [19], genetic algorithm (GA) [20], Whale Optimization Algorithm (WOA) [21], Artificial Bee Colony (ABC) [22] are used for EH energy management.

Demand response programs (DRPs) represent a pivotal method in smart grid energy management, guiding consumers toward optimal consumption patterns through dynamic electricity pricing or incentive policies. While DRPs yield substantial benefits for power systems and economic advantages for consumers, they introduce complexities in EH management systems [23]. For instance [24], proposed optimal planning for an EH integrating photovoltaic (PV) distributed generation, compressed air storage, and demand-responsive loads to minimize operating costs, employing metaheuristic algorithms in a risk-aware market framework. However, effective EH management must account for uncertainties, particularly in renewable generation and load estimation due to variable weather conditions—a critical yet often underemphasized factor in studies prioritizing power system stability and financial metrics [25]. These uncertainties lead to deviations from optimal plans, resulting in operational inefficiencies [26]. Weather variability not only impacts renewable sources (e.g., PV and wind output) but also alters load profiles, especially in household demand linked to heating, ventilation, and air conditioning (HVAC) systems. Despite this, most studies model wind and solar uncertainty using normal distributions without detailed parameter reporting [27], while overlooking weather-derived load forecasting improvements [28]. To address these gaps [29], leveraged machine learning tools for robust energy management of electrical and thermal resources, explicitly incorporating uncertainties in load predictions.

This study presents a probabilistic energy management model for an EHthat includes renewable energy sources, hydrogen energy storage systems. The model aims to optimize operations by minimizing costs and increasing flexibility, while also considering DRPs and uncertainties. A mutated multi-objective grasshopper optimization algorithm is proposed. In the modified version of the multi-objective grasshopper algorithm, a genetic mutation is introduced for weaker grasshoppers, enhancing the algorithm’s precision in identifying optimal points. The main contributions of this paper include:

Innovative EMS: An advanced EMS is developed, integrating renewable energy sources, electric vehicle (EV) parking, and hydrogen-based storage. It also considers DRPs and operational uncertainties, aiming to reduce costs and improve system adaptability.

Comprehensive cost-based and emission objective functions: A detailed set of objective functions is proposed to evaluate various operational scenarios. These scenarios are weighted according to their probability, providing a realistic representation of practical conditions.

Novel optimization algorithm: A mutated multi-objective grasshopper optimization algorithm is introduced to tackle EH management challenges. This enhanced algorithm effectively navigates a broader solution space and reduces the likelihood of convergence to local optima.

In continuation of the article, Section 2 presents the problem formulation, including a stochastic model for the EH, scenario selection using the roulette wheel mechanism (RWM), and the objective functions and constraints for optimization. Section 3 introduces the mutant multi-objective grasshopper optimization algorithm used to solve the EH energy management problem. Section 4 evaluates the simulation results, first by assessing modified multi-objective grasshopper optimization algorithm’s (MMOGOA) accuracy in comparison with other optimization methods, and then, by analyzing EH energy management under different scenarios—without uncertainties and demand response program (DRP), with DRP but without uncertainties, and with both DRP and uncertainties. The conclusion in Section 5 highlights key findings and suggests future research into resilience and reliability in EHs for more robust and efficient energy management.

2. Problem Formulation

2.1. Stochastic Model of Hub Energy

The EH under study, depicted in Figure 1, integrates renewable energy systems including PV panels and wind turbines (WTs). In addition, it incorporates a boiler, a combined heat and power (CHP) unit, thermal storage, and EV parking facilities for energy storage purposes. The EH is also equipped with an electrolyzer, enabling the use of surplus electricity to convert water into hydrogen. When needed, the stored hydrogen undergoes a methanation process to produce natural gas (CH₄), which can then be used as fuel in the CHP unit or the boiler.

[figure(s) omitted; refer to PDF]

Electric and thermal loads, as the primary sources of uncertainty in the EH, play a crucial role in the energy management program. Variations in these loads can arise from temporal fluctuations and weather conditions. Daily and weekly changes in electricity demand are largely influenced by consumer behavior patterns. To model the uncertainty in loads, a normal probability density function (PDF) can be employed as a suitable approach to represent load variations in the EH. The PDF for thermal and electrical load fluctuations is expressed by Equation (1), which applies equally to both types of loads [30].$$\begin{array}{}\text{(1)}& f{P}_L=\frac{1}{\sqrt{2\pi }\times {\sigma }_{PL}}\exp -\frac{{P_L}^{-}}{2\times {\sigma }_{PL}^2}.\end{array}$$

In Equation (1), p_L, μ_PL, and σ_PL represent the EH load, the average and standard deviation of heat or electrical load, respectively. These parameters are calculated based on the recorded load data from the EH. In addition to load uncertainties, the uncertainties associated with renewable energy sources, specifically solar and wind energy, must be accurately incorporated into the proposed EH model. The intermittent nature of wind turbine (WT) power generation introduces unique challenges for energy management within the EH. The Weibull PDF is a suitable mathematical model for representing wind speed distribution.$$\begin{array}{}\text{(2)}& f_{WT}\upsilon =\begin{array}{cc}\frac{k}{p}\times {\frac{\upsilon }{p}}^{k-1}\times \exp -{\frac{\upsilon }{p}}^k& \upsilon \ge 0\\ 0\hfill & \text{otherwise}\end{array}.\end{array}$$

In Equation (2), υ is a wind speed and also, p and k are scaling and shaping factors of the Weibull PDF. The power generation of the WT is calculated based on Equation (3).$$\begin{array}{}\text{(3)}& Pt=\begin{array}{cc}0\hfill & \upsilon t<{\upsilon }_i\text{and}\upsilon t>{\upsilon }_o\\ P_r\frac{\upsilon -{\upsilon }_i}{{\upsilon }_r-{\upsilon }_i}\hfill & {\upsilon }_i\le \upsilon t\le {\upsilon }_r\\ P_r\hfill & {\upsilon }_r\le \upsilon t\le {\upsilon }_o\end{array}.\end{array}$$

In Equation (3), υ is wind speed (m/s) and ${\upsilon }_i,{\upsilon }_o,\text{\hspace{0.17em}and}{\upsilon }_r$are the WT cut in, cutoff and nominal speed. $Pt$ is the output power generation of the WT and $P_r$ is the WT rated power. The power output of a PV array is directly influenced by the solar irradiance incident on its surface. Solar irradiance exhibits probabilistic characteristics, leading to variability and uncertainty in the power generation of the PV system. The beta PDF serves as an appropriate model for characterizing the fluctuations in solar irradiance.$$\begin{array}{}\text{(4)}& f_{pv}si=\begin{array}{cc}\frac{\Gamma \alpha +\beta }{\Gamma \alpha \Gamma \beta }{si}^{\alpha -1}{1-si}^{\beta -1}& \alpha >0,\beta >0\\ 0\hfill & o.w\end{array}.\end{array}$$Here, $\text{si}$ signifies the solar irradiation, Г refers to the Gamma function, and $\alpha \text{\hspace{0.17em}and}\beta$ are the defining parameters of the beta distribution. Equation (5) is used to determine the power generation capacity of the PV system.$$\begin{array}{}\text{(5)}& P_{PV}t=P_{PV}^R\times \frac{Sit}{1000}\times {\eta }_{PV},\end{array}$$where η_PV represents the efficiency of the PV system, and $P_{\text{PV}}^R$denotes the rated capacity of the PV system.

Modeling the input and output power in parking lots (PLs) necessitates consideration of various uncertainties, such as EV charging and discharging schedules, battery types and capacities, the fraction of EVs using the parking facility, driver behaviors, and the timing of EV arrivals and departures, along with the state of charge (SOC) of the batteries. To effectively capture the statistical characteristics of EV arrival and departure times, the Weibull and Lognormal probability density functions (PDFs) are deemed the most suitable models. These functions accurately represent the variability in EV behavior, facilitating a comprehensive analysis of charging trends. Equations (6) and (7) outline these PDFs for modeling EV dynamics.

Due to the use of renewable energy sources in the studied EH and the stochastic nature of their power generation, incorporating a battery storage system is essential. The battery capacities of EVs parked in PLs can be utilized for electricity storage. However, uncertainties regarding the presence of EVs in the PLs and their SOC pose significant challenges to using EVs for energy storage. The modeling of input/output power for PLs must consider various sources of uncertainty, including the scheduling of EV charging and discharging, the battery type and capacity in EVs, the proportion of EVs in PLs, driver behavior, EV arrival and departure times, and SOC of the batteries. The Weibull and Lognormal PDFs (PDFs) are considered the most suitable models for capturing the arrival and departure times of EVs. These PDFs provide a good fit to the statistical distribution of EV arrival and departure patterns, allowing for effective modeling and analysis of EV charging behavior. Equations (6) and (7) are employed to represent the PDFs for EV arrival and departure [31].$$\begin{array}{}\text{(6)}& fArr=\frac{k}{\Omega }{\frac{Arr}{\Omega }}^{k-1}{e}^{-{Arr/\Omega }^k},\end{array}$$$$\begin{array}{}\text{(7)}& fDep=\frac{1}{Dep\times \sigma \times \sqrt{2\pi }}{e}^{-{\ln Dep-\mu }^2/2{\sigma }^2}.\end{array}$$

2.2. Scenarios Selection by RWM

In stochastic problems, selecting appropriate scenarios is of critical importance. In this study, the RWM is employed to select various scenarios for the proposed objective functions. In the roulette wheel theory, each uncertainty factor is divided into sections based on its PDF, and each section is assigned a value. These values can be normalized by dividing each section’s value by the base value. Random selections are then made using the roulette wheel. This mechanism increases the likelihood of selecting more probable responses, where the probability in the PDF is higher compared to other points. However, less probable points are also, given a chance since the probability of the most likely outcomes is less than one. In the roulette wheel method, a circle, as depicted in Figure 2, is divided into sections proportional to the existing probabilities [32].

[figure(s) omitted; refer to PDF]

In Figure 2, the area of each section is determined based on the existing probabilities and are expressed as follows.$$\begin{array}{}\text{(8)}& Ai=\text{Prob}i\times \pi r^2.\end{array}$$

If a pointer is present and the circle is spun randomly, the probability of the pointer stopping at section $i$ is equal to $P_i$. To apply this method to the problem, the circle is assumed to be represented as a line segment ranging from 0 to 1, with probabilities distributed along it as shown in Figure 3. To sample from this distribution, a random number is generated uniformly between 0 and 1.

[figure(s) omitted; refer to PDF]

The uniformity ensures that the chance of selecting each section is proportional to its length. In other words, the generated random number falls into section $P_i$ with a probability equal to $P_i$. The generated random number must then be sequentially checked against the defined conditions to determine the section. The RWM efficiently generates scenarios from PDFs by probabilistically selecting high-likelihood events while retaining rare but critical cases. It offers computational scalability without losing key dynamics and preserves statistical accuracy, making it ideal for renewable energy studies requiring robust uncertainty modeling of wind/solar generation and load fluctuations.

2.3. Objective Functions and Constraints

Within the framework of probabilistic analysis, the objective functions in this study are expressed in the form $y=fZ$, where $Z$ represents the vector of the problem’s random variables. Accordingly, for the probabilistic analysis, the objective function is calculated using Equation (10).$$\begin{array}{}\text{(9)}& y=\sum _{i=1}^{N_s}\text{Prob}i\times fZ,\end{array}$$where $\text{Prob}i$ represents the probability of each scenario, derived from the PDFs. The first objective function for this paper is the operational cost function, which is defined as Equation (10).$$\begin{array}{}\text{(10)}& f_1=\mathrm{Min}\sum _{i=1}^{N_s}\pi i{E}_{\text{Cost}}i+H_{\text{Cost}}i+T_{\text{Cost}}i.\end{array}$$

In Equation (10), $E_{\text{Cost}}$ represents the cost of requested electrical energy, $H_{\text{Cost}}$ represents the cost of requested thermal energy, and $T_{\text{Cost}}$ is the cost of converting electrical energy to hydrogen cost. Additionally, $\text{Prob}i$ denotes the probability of scenario $i$, and $N_s$ is the total number of scenarios considered for the study. The inclusion of the probability coefficient ensures that more probable scenarios have a greater impact on the final value of the objective function compared to less likely scenarios. The cost of requested electrical energy is calculated as shown in Equation (11).$$\begin{array}{}\text{(11)}& E_{\text{Cost}}=\sum _{t=1}^{24}{C}_{\text{NG}}t+C_{\text{PV}}t+C_{\text{WT}}t+C_{\text{CHP-E}}t+C_{\text{DRP-E}}t+C_{\text{PL}}t.\end{array}$$where $C_{\text{NG}},C_{\text{PV}},C_{\text{WT}},C_{\text{CHP-E}},C_{\text{PL}},\text{\hspace{0.17em}and}{C}_{\text{DRP-E}}$ represent the costs of power purchased from the grid, the costs of PV and wind distributed generation, the cost of electrical power generated by the CHP unit, DRP cost, and EV parking lot (PL) costs, respectively. The cost of purchasing energy from the grid is calculated using Equation (12).$$\begin{array}{}\text{(12)}& C_{\text{NG}}t=\varphi t\times P_{\text{NG}}t.\end{array}$$where $P_{\text{NG}}t$ represents the power purchased during time interval t,$\phi t$ denotes the cost per kilowatt of energy exchanged with the grid. The costs of power generated by renewable energy sources and the electrical power produced by the CHP unit are calculated using Equations (13)–(15), respectively.$$\begin{array}{}\text{(13)}& C_{\text{PV}}t={\sigma }_{\text{PV}}\times P_{\text{PV}}t,\end{array}$$$$\begin{array}{}\text{(14)}& C_{\text{WT}}t={\sigma }_{\text{WT}}\times P_{\text{WT}}t,\end{array}$$$$\begin{array}{}\text{(15)}& C_{\text{CHP-E}}t={\sigma }_{\text{CHP}}\times P_{\text{CHP}}t.\end{array}$$where $P_{\text{PV}}t,P_{\text{WT}}t,\text{\hspace{0.17em}and}{P}_{\text{CHP-E}}t$ represent the power generated by PV systems, WT, and CHP units, respectively. Similarly, ${\sigma }_{\text{PV}},{\sigma }_{\text{WT}},\text{\hspace{0.17em}and}{\sigma }_{\text{CHP}}$ denote the cost coefficients associated with these distributed generation resources. The cost of responsive electrical loads participating in the demand response program can be calculated using Equation (16).$$\begin{array}{}\text{(16)}& C_{\text{FL}-\mathrm{E}}t=\zeta t\times P_{\text{DRP}}t.\end{array}$$where $P_{\text{DRP}}t$ represents the amount of reduced electricity consumption, and $\xi t$ is the reward for load reduction at time t. Considering DRP and reducing energy consumption during peak hours can significantly contribute to lowering operational costs. Finally, the cost associated with PLs is determined using Equation (17).$$\begin{array}{}\text{(17)}& C_{\text{PL}}t=\vartheta t\times P_{\text{Ch}}t+P_{\text{dCh}}t.\end{array}$$In the above equation, $P_{\text{Ch}}t\text{and}{P}_{\text{dCh}}t$represent the input and output power of the PL, respectively, while $\vartheta t$ denotes the cost per kilowatt of electricity exchanged with the PL during hour t. The cost of the requested heat energy is calculated using Equation (18):$$\begin{array}{}\text{(18)}& H_{\text{Cost}}=\sum _{t=1}^{24}{C}_{\text{HO}}t+C_{\text{CHP-H}}t+C_{\text{DRP-H}}t+C_{\text{H-Storage}}t,\end{array}$$where $C_{\text{CHP-H}}$ represents the thermal cost of the CHP unit, $C_{\text{HO}}$ refers to the cost of the heat-only (HO) unit or boiler, $C_{\text{DRP-H}}$denotes the cost of thermal responsive loads participating in the DRP, and $C_{\text{H-Storage}}$ indicates the cost of thermal energy storage. The cost of the heating-only unit $C_{\text{HO}}$ is calculated using Equation (19).$$\begin{array}{}\text{(19)}& C_{\text{HO}}t=C_{\text{HO-F}}t+C_{\text{HO-SU/SD}}t,\end{array}$$where $C_{\text{HO-F}}$ represents the fuel cost, and $C_{\text{HO-SU/SD}}$ denotes the startup/shutdown cost. The fuel cost of the HO unit $C_{\text{HO-F}}$ is calculated using Equation (20).$$\begin{array}{}\text{(20)}& C_{\text{HO-F}}t={\alpha }_{\text{HO}}{H}_{\text{HO}}^{2}t+{\beta }_{\text{HO}}{H}_{\text{HO}}t+{\delta }_{\text{HO}},\end{array}$$where $H_{\text{HO}}$ represents the heat generated by the boiler, and ${\alpha }_{\text{HO}},{\beta }_{\text{HO}},\text{\hspace{0.17em}and}{\delta }_{\text{HO}}$ are the cost coefficients of the boiler. The startup/shutdown cost is calculated using Equation (21).$$\begin{array}{}\text{(21)}& C_{\text{HO-SU/SD}}t=SUC\times Ut-Ut-1,\end{array}$$where SUC represents the startup/shutdown cost, which is calculated for cold and hot startups using Equation (22).$$\begin{array}{}\text{(22)}& SUC=\begin{array}{c}HSC{T}_{\text{off}}\le CST+MDT\\ CSC{T}_{\text{off}}>CST+MDT\end{array},\end{array}$$where HSC and CSC represent the hot startup cost and cold startup cost, respectively. Additionally, CST is the time required for a cold startup, and MDT is the minimum downtime of the boiler. The heat generation cost by the CHP unit is expressed as Equation (23).$$\begin{array}{}\text{(23)}& C_{\text{CHP-H}}t={\tau }_{\text{CHP}}\times H_{\text{CHP}}t.\end{array}$$Here, $H_{\text{CHP}}$ represents the heat generated by the CHP unit, and ${\tau }_{\text{CHP}}$ is the cost per kilowatt of heat produced by this unit. Similar to electrical loads, thermal loads are categorized into flexible loads and other loads. The cost of participation of the loads in the DRP can be calculated using Equation (24).$$\begin{array}{}\text{(24)}& C_{\text{DRP-H}}t=\lambda t\times H_{\text{DRP}}t,\end{array}$$where H_DRP(t) represents the amount of reduced thermal energy, and λ(t) denotes the reward for consumption reduction at time t. Finally, the cost of the thermal storage system is calculated using Equation.$$\begin{array}{}\text{(25)}& C_{\text{H-Storage}}t=\epsilon t\times H_{\text{Ch}}t+H_{\text{dCh}}t.\end{array}$$

The cost of converting electrical energy into hydrogen, storing it, and producing natural gas is defined by Equation (26).$$\begin{array}{}\text{(26)}& T_{\text{Cost}}=\sum _{t=1}^{24}{\varphi }_T\times P_{T}t+{\theta }_{\text{Storage-G}}\times G_{\text{Storage-G}}t,\end{array}$$where ${\varphi }_T$ refers to the cost associated with the electrolyzer, which is responsible for converting electrical energy into hydrogen gas, while ${\theta }_{\text{G-Storage}}$ represents the cost of storing the produced gas. Furthermore, the emission function is introduced as the second objective in the analysis.$$\begin{array}{}\text{(27)}& f_2=\mathrm{Min}\sum _{i=1}^{N_k}\pi i\sum _{t=1}^{24}{P}_{\text{CHP}}t,i+H_{\text{CHP}}t,i\times E_{\text{CHP}}+H_{\text{HO}}t,i\times E_{\text{HO}}+P_{\text{NG}}t,i\times E_{\text{NG}},\end{array}$$Here, $E_{\text{CHP}},E_{\text{HO}},\text{and}{E}_{\text{NG}}$ represent the emissions from the CHP unit, the HO unit, and the grid network. These emissions include carbon dioxide ${\text{CO}}_2$, sulfur dioxide ${\text{SO}}_2$, and nitrogen oxides ${\text{NO}}_x$ released from each of the units. The optimization algorithm aims to maximize the use of low-emission renewable energy sources during most hours. Reducing production costs generally leads to an increase in environmental emissions, and vise versa. This implies that it is not possible to minimize both costs and emissions simultaneously. The main constraints for the EH operation are as follows [33]:

Electric power balance: The amount of generated and requested electric power must be equal in each hour.$$\begin{array}{}\text{(28)}& P_{\text{NG}}t+P_{\text{CHP}}t+P_{\text{WT}}t+P_{\text{PV}}t+P_{\text{dCh}}t=P_{\text{Load}}t+P_{\text{Loss}}t+P_{\text{Ch}}t+P_{T}t-P_{\text{DRP}}t.\end{array}$$

Thermal power balance constraint: The amount of generated and requested thermal power must be equal in each hour.$$\begin{array}{}\text{(29)}& H_{\text{CHP}}t+H_{\text{HO}}t=H_{\text{Load}}t+H_{\text{Loss}}t-H_{\text{DRP}}t.\end{array}$$

Power grid exchange power constraint: The limitation for the electric power exchanged with the power grid is as follows.$$\begin{array}{}\text{(30)}& P_{\text{NG}}^{\mathrm{Min}}\le P_{\text{NG}}t\le P_{\text{NG}}^{\mathrm{Max}}.\end{array}$$

CHP operation constraint: The operational constraints of the CHP unit are depicted in Figure 4.

[figure(s) omitted; refer to PDF]

Equations (31)–(39) express the formulation of CHP feasible areas$$\begin{array}{}\text{(31)}& P_{j,t}^{\text{CHP}}-P_{j,B}^{\text{CHP}}-\frac{P_{j,B}^{\text{CHP}}-P_{j,C}^{\text{CHP}}}{H_{j,B}^{\text{CHP}}-H_{j,C}^{\text{CHP}}}\times P_{j,t}^{\text{CHP}}-P_{j,B}^{\text{CHP}}\le 0,\end{array}$$$$\begin{array}{}\text{(32)}& P_{j,t}^{\text{CHP}}-P_{j,C}^{\text{CHP}}-\frac{P_{j,C}^{\text{CHP}}-P_{j,D}^{\text{CHP}}}{H_{j,C}^{\text{CHP}}-H_{j,D}^{\text{CHP}}}\times H_{j,t}^{\text{CHP}}-H_{j,C}^{\text{CHP}}\le 0,\end{array}$$$$\begin{array}{}\text{(33)}& P_{j,t}^{\text{CHP}}-P_{j,F}^{\text{CHP}}-\frac{P_{j,E}^{\text{CHP}}-P_{j,F}^{\text{CHP}}}{H_{j,E}^{\text{CHP}}-H_{j,F}^{\text{CHP}}}\times H_{j,t}^{\text{CHP}}-H_{j,F}^{\text{CHP}}\ge -1-X_{j,t,1}^{\text{CHP}}\times M,\end{array}$$$$\begin{array}{}\text{(34)}& P_{j,t}^{\text{CHP}}-P_{j,D}^{\text{CHP}}-\frac{P_{j,D}^{\text{CHP}}-P_{j,E}^{\text{CHP}}}{H_{j,D}^{\text{CHP}}-H_{j,E}^{\text{CHP}}}\times H_{j,t}^{\text{CHP}}-H_{j,D}^{\text{CHP}}\ge -1-X_{j,t,2}^{\text{CHP}}\times M,\end{array}$$$$\begin{array}{}\text{(35)}& 0\le P_{j,t}^{\text{CHP}}\le P_{j,A}^{\text{CHP}}\times V_{j,t}^{\text{CHP}},\end{array}$$$$\begin{array}{}\text{(36)}& 0\le H_{j,t}^{\text{CHP}}\le H_{j,B}^{\text{CHP}}\times V_{j,t}^{\text{CHP}},\end{array}$$$$\begin{array}{}\text{(37)}& X_{j,t,2}^{\text{CHP}}+X_{j,t,2}^{\text{CHP}}=V_{j,t}^{\text{CHP}},\end{array}$$$$\begin{array}{}\text{(38)}& H_{j,t}^{\text{CHP}}+H_{j,F}^{\text{CHP}}\le 1-X_{j,t,1}^{\text{CHP}}\times M,\end{array}$$$$\begin{array}{}\text{(39)}& H_{j,t}^{\text{CHP}}+H_{j,F}^{\text{CHP}}\ge -1-X_{j,t,2}^{\text{CHP}}\times M,\end{array}$$where $P_{j,t}^{\text{CHP}}$ is the unit electrical output and $H_{j,t}^{\text{CHP}}$ is the unit heat power output at t.

HO operation constraint: The heat production capacity of HO units is limited.$$\begin{array}{}\text{(40)}& H_{\text{HO}}^{\mathrm{Min}}\le H_{\text{HO}}^{k}t\le H_{\text{HO}}^{\mathrm{Max}}.\end{array}$$

EV PL constraint: The charging and discharging rates, as well as the charging and discharging levels of EVs in the PLs, are limited.$$\begin{array}{}\text{(41)}& P_{\text{dCh}}t\le P_{\text{dCh}}^{\mathrm{Max}},\end{array}$$$$\begin{array}{}\text{(42)}& P_{\text{Ch}}t\le P_{\text{Ch}}^{\mathrm{Max}},\end{array}$$$$\begin{array}{}\text{(43)}& E_{\text{Ch}}t\le {\psi }_k^{\text{charge-max}},\end{array}$$$$\begin{array}{}\text{(44)}& E_{\text{dCh}}t\ge {\psi }_k^{\text{charge-min}}.\end{array}$$

3. Mutant Multiobjective Grasshopper Algorithm

Nature-inspired algorithms logically divide the search process into two primary phases: exploration and exploitation. During the exploration phase, search agents are encouraged to move abruptly and explore the search space extensively, whereas, in the exploitation phase, they focus on making localized movements to refine potential solutions. The MOGOA algorithm is inspired by the natural foraging behavior of grasshoppers and strives to balance exploration and exploitation in multi-objective optimization problems. This algorithm is designed to identify a set of Pareto-optimal solutions for multi-objective optimization challenges. To simulate the collective behavior of grasshoppers, it is formulated as follows [34]:$$\begin{array}{}\text{(45)}& X_i=S_i+G_i+A_i,\end{array}$$where $X_i$ represents the position of the $i^{\text{th}}$ grasshopper, $S_i$ denotes the social interaction effect, $G_i$ signifies the gravitational force acting on $i^{\text{th}}$ grasshopper, and $A_i$ represents the advection caused by wind. To reduce randomness in behavior, this equation can be expressed as:$$\begin{array}{}\text{(46)}& X_i={\mathrm{rand}}_1{S}_i+{\mathrm{rand}}_2{G}_i+{\mathrm{rand}}_3{A}_i,\end{array}$$where ${\mathrm{rand}}_1,{\text{\hspace{0.17em}rand}}_2,{\text{\hspace{0.17em}and\hspace{0.17em}rand}}_3$ are random numbers in the range [0, 1].$$\begin{array}{}\text{(47)}& S_i=\sum _{\begin{array}{c}i=1\\ j\ne i\end{array}}^{N}s{d}_{ij}{d}_{ij}^{\ast },\end{array}$$where $d_{ij}$represents the distance between the $i^{\text{th}}$ and $j^{\text{th}}$ grasshoppers, and $d_{ij}^{\ast }$ is the unit vector pointing from the $i^{\text{th}}$ grasshopper to the $j^{\text{th}}$ grasshopper [34].$$\begin{array}{}\text{(48)}& d_{ij}=x_i-x_j,\end{array}$$$$\begin{array}{}\text{(49)}& d_{ij}^{\ast }=\frac{x_i-x_j}{d_{ij}}.\end{array}$$

The function s is used to define the strength of social forces which determines the intensity of social interactions, is computed as follows:$$\begin{array}{}\text{(50)}& sr=f{e}^{\frac{-r}{l}}-e^{-r},\end{array}$$where f represents the attraction intensity, and l is the characteristic length scale of the attraction. Additionally, the gravitational component G is calculated as:$$\begin{array}{}\text{(51)}& G_i=-g{e}_g^{\ast },\end{array}$$where g is the gravitational constant, and $e_g^{\ast }$ is a unit vector pointing toward the Earth’s center. The advection component A is computed as:$$\begin{array}{}\text{(52)}& A_i=u{e}_w^{\ast },\end{array}$$where u is a drift constant, and $e_w^{\ast }$ is a unit vector in the direction of the wind. Since nymph grasshoppers lack wings, their movement is heavily influenced by the wind’s direction. Substituting S, G, and A into Equation (39), the expanded equation becomes:$$\begin{array}{}\text{(53)}& X_i=\sum _{\begin{array}{c}i=1\\ j\ne i\end{array}}^{N}s{x}_i-x_j\frac{x_i-x_j}{d_{ij}}-g{e}_g^{\ast }+u{e}_w^{\ast },\end{array}$$where N is the total number of grasshoppers. Since nymph grasshoppers reside on the ground, their positions should not fall below a certain threshold. However, this equation is not applied directly in group simulation or optimization algorithms because it prevents the algorithm from fully exploring and exploiting the search space around a solution. In fact, the model assumes a group in free space. This mathematical model cannot be directly used for solving optimization problems, primarily because grasshoppers rapidly settle in a comfortable region, and the swarm does not converge to a specific point. In the proposed version of the MOGOA algorithm, genetic mutations are introduced for a subset of grasshoppers to enhance diversity and prevent premature convergence. A percentage of the grasshoppers $M_{\text{f}}$ that perform poorly compared to the rest of the population are selected for genetic mutations. These mutations modify their positions in the search space to increase the likelihood of escaping local optima and exploring new regions.$$\begin{array}{}\text{(54)}& X_i=X_i+randn\times \sigma ,\end{array}$$where randn is a random number sampled from a uniform distribution in the range [−1, 1] and σ is a scaling factor controlling the mutation’s intensity. The mutation diversifies the search process and enables the algorithm to explore areas that might not be reachable through standard MOGOA operations. The flowchart of the proposed MMOGOA is shown in Figure 5.

[figure(s) omitted; refer to PDF]

4. Evaluation of Simulation Results

The simulation results are analyzed in two main parts. The first part assesses the effectiveness of the MMOGOA algorithm by applying it to benchmark functions and comparing its performance with that of other multi-objective optimization algorithms. The second part focuses on the energy management problem within the EH, which is further subdivided into three sections. In this part, the MMOGOA algorithm is utilized along with several other optimization techniques to solve the energy management issue in the EH. The simulations were carried out using MATLAB 2020a on an ACER laptop, equipped with an Intel i9 13900HX processor, 16 GB of RAM, 1TB storage, and an RTX 4060 graphics card.

4.1. The MMOGOA Accuracy Evaluation

The initial phase of the simulation involved assessing the performance of the proposed MMOGOA algorithm through the optimization of ZDT benchmark functions. The results were then compared against those achieved by several other algorithms, including the Multi-Objective Exponential Distribution Optimizer (MOEDO) [35], the Multi-Granularity Clustering-Based Evolutionary Algorithm (MGCBEA) [36], the Multistage Competitive Swarm Optimization Algorithm (MSCSOA) [37], and the Multi-Strategy Multi-Objective Differential Evolutionary Algorithm (MSMODE) [38]. After performing 20 independent runs to optimally solve the ZDT benchmark problems, various performance metrics, including generation distance (GD), spacing (SP), and inverse generation distance (IGD), to thoroughly assess the algorithm’s effectiveness. The optimization results are shown in Figure 6.

[figure(s) omitted; refer to PDF]

Smaller values for the calculated criteria across multiple iterations indicate improved performance and greater precision of the optimization method. As shown in Figure 5, the proposed MMOGOA algorithm consistently produces lower average values for these criteria compared to other algorithms across most benchmark functions. Furthermore, reduced standard deviation values suggest that the algorithm reliably identifies optimal solutions that are closely clustered together, rather than, dispersing randomly around the optimal point

4.2. The EH Energy Management

Figure 7 presents the hourly price variations for electricity and natural gas in the energy market [39].

[figure(s) omitted; refer to PDF]

Using the RWM and PDFs, ten distinct scenarios have been generated for the electrical and thermal loads within the microgrid, as well as for the power outputs of renewable energy sources. The corresponding values for these scenarios are depicted in Figure 8.

[figure(s) omitted; refer to PDF]

The arrival pattern of electric vehicles (EVs) to the PL is modeled using a Weibull PDF with parameters c = 8.4 and k = 2.8. Conversely, the departure pattern is characterized by a lognormal PDF with a mean of 2.4 and a standard deviation of 0.4. Figure 9 illustrates the number of EVs arriving at and departing from the PL based on these probability distributions.

[figure(s) omitted; refer to PDF]

The energy management of the EH is divided into three stages. In the first stage, energy management is carried out without considering uncertainties and the DRP. In the second stage, uncertainties are incorporated, but the DRP is still ignored. Finally, in the third stage, both uncertainties and the DRP are accounted for in the EH energy management. To assess the performance of the MMOGOA algorithm, its optimization results were compared with those obtained from well-established methods, including the Non-dominated Sorting Genetic Algorithm (NSGA-II), Multi-Objective Gray Wolf Optimization (MOGWO), and the Multi-Objective Whale Optimization Algorithm (MOWOA). The optimization algorithms parameters are included in Table 1.

Table 1

The optimization algorithms parameters.

MOPSO	MMOGOA	MOWOA	NSGA-II	MOPSO	MMOGOA
$\text{Pop}.=100$	$\text{Max\hspace{0.17em}Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$	$\mathrm{Max}\text{Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$
$\text{Max\hspace{0.17em}Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$	$\text{Max\hspace{0.17em}Iter}.=50$
$\text{Inertia\hspace{0.17em}weight},w=0.4\to 0.9$	$c_{\max }=1.0,c_{\min }=0.00001$	$\text{Spiral}\text{constant},b=1$	$\text{Crossover}\text{prob}.,p_a=0.9$	$\text{Inertia\hspace{0.17em}weight},w=0.4\to 0.9$	${\text{c}}_{\max }=1.0,{\text{c}}_{\min }=0.00001$
${\text{c}}_1=1.5,{\text{\hspace{0.17em}c}}_2=2.0$	$f=0.5,l=1.5$	$a=2\to 0\text{linear\hspace{0.17em}decrease}$	$\text{Mutation\hspace{0.17em}prob}.,p_m=0.2$	${\text{c}}_1=1.5,{\text{\hspace{0.17em}c}}_2=2.0$	$f=0.5,l=1.5$
$\text{Repository\hspace{0.17em}size}=100$	$\text{Mutation\hspace{0.17em}prob}.=0.1$	$\text{Prob.\hspace{0.17em}threshold},p=0.5$	${\eta }_a=20,{\eta }_m=20$	$\text{Rep\hspace{0.17em}size}=100$	$\text{Mutation\hspace{0.17em}prob}.=0.1$
$\text{Grid\hspace{0.17em}size}=30$	$\text{Grid\hspace{0.17em}size}=30$	$\text{Crossover}\text{rate}=0.8$	$\text{Tournament}\text{size}=2$	$\text{Grid\hspace{0.17em}divisions}=30$	$\text{Grid\hspace{0.17em}size}=30$

• Without considering uncertainties and DRP

In the first stage of the simulation analysis, the energy management of the EH was carried out without incorporating uncertainties or DRPs. The Pareto fronts for this phase of the simulations, shown in Figure 10.

[figure(s) omitted; refer to PDF]

The final selection of the best result is made using a fuzzy method. The results are accumulated in Table 2.

Table 2

The best results in the first section of the simulation.

Optimization algorithm	$E_{\text{Cost}}¢$	$H_{\text{Cost}}¢$	$T_{\text{Cost}}¢$	$\text{Total\hspace{0.17em}cost}¢$	$\text{Emission}\text{kg}$
MMOGOA	1133	767	325	2225	15.32
MOPSO	1156	808	283	2247	16.83
NSGA-II	1152	794	295	2241	17.36
MOWOA	1146	781	311	2238	16.41

The proposed MMOGOA algorithm achieves an EH operation cost of 2225 cents and emissions of 15.32 kg, showcasing superior performance compared to other optimization methods. For reference, the costs for MOPSO, NSGA-II, and MOWOA are 2247, 2241, and 2238 cents, with corresponding emission levels of 16.83, 17.36, and 16.41 p.u. respectively. Figure 11 provides an overview of the hourly charge variations across hydrogen storage, EV parking, and heat storage, demonstrating effective management of energy resources over time.

[figure(s) omitted; refer to PDF]

To achieve cost-effective energy management, the operator employs a strategic approach to energy storage and conversion. During low-demand periods with reduced electricity prices, excess power is acquired from the grid and either converted to natural gas or stored in EV batteries in the PL. This stored natural gas is reserved for future use. When natural gas prices rise, the stored reserves are utilized within the EH to produce both electricity and heat, meeting system demands efficiently and cost-effectively. This dynamic process leverages price fluctuations to optimize energy usage and reduce overall expenses. The operational analysis shows a significant increase in electricity drawn from the grid between 1:00 AM and 7:00 AM, primarily for EV charging and hydrogen production. During high electricity price hours, reliance on grid electricity decreases as EV batteries discharge and stored gas powers the CHP system. A portion of the stored gas is directed toward heating, while the excess is converted back into electricity. Additionally, approximately 11.32% of the EHs gas demand is supplied by stored hydrogen, underlining its importance in reducing reliance on external resources.

• Considering DRP without uncertainties

In the second phase of the simulation, energy management of the EH was optimized by incorporating DRPs into the framework, while uncertainties in system parameters were deliberately excluded to simplify the analysis. This phase assumed that approximately 30% of the total electrical and thermal loads were flexible and could participate in the DRP, enabling the system to shift or reduce energy demand based on pricing incentives or operational needs. By integrating DRPs, the EH could better align its load profile with available energy resources, improving efficiency and reducing costs. The results of this simulation are visually represented in Figure 12, where the pareto fronts obtained from various optimization algorithms demonstrate the trade-offs achieved between multiple objectives. These fronts highlight the system’s ability to balance competing factors, such as minimizing operational costs and emissions, while maintaining reliable energy supply to meet demand. This analysis emphasizes the role of DRPs in enhancing system flexibility and overall performance.

[figure(s) omitted; refer to PDF]

The optimal outcomes for this phase of the simulations are summarized in Table 3, where the best results from the various optimization algorithms are presented for comparison.

Table 3

The best results in the second section of the simulation.

Optimization algorithm	$E_{\text{Cost}}¢$	$H_{\text{Cost}}¢$	$T_{\text{Cost}}¢$	$\text{Total\hspace{0.17em}cost}¢$	$\text{Emission}\text{kg}$
MMOGOA	986	534	285	1805	13.08
MOPSO	1013	578	241	1832	13.23
NSGA-II	983	591	264	1838	13.31
MOWOA	1006	529	276	1811	13.16

The simulation results in this section highlight the positive impact of implementing a DRP within the EMS. The inclusion of DRP led to an 18.87% reduction in operating costs and a 14.62% decrease in emissions at the EH. Among the algorithms tested, the MMOGOA-based EMS delivered the most cost-effective and environmentally friendly outcomes. After optimization with the MMOGOA algorithm, the total operating cost was reduced to 1805 cents. In comparison, the costs using EMS with the MOPSO, NSGAII, and MOWOA algorithms were slightly higher, at 1832, 1838, and 1811 cents, respectively. Similarly, emissions were minimized with the MMOGOA-based EMS, yielding a value of 13.08 kg. By contrast, the emissions with EMSs optimized by the MOPSO, NSGAII, and MOWOA algorithms were 13.23, 13.31, and 13.16 kg, respectively. Figure 13 illustrates the hourly charge levels of the storage systems for this part of the simulations, showcasing how the DRP influenced energy storage and utilization patterns throughout the day.

[figure(s) omitted; refer to PDF]

Figure 10 illustrates that incorporating the DRP into the EMS effectively reduces the consumption of electrical and heat loads in the EH. This decreased demand directly lowers the reliance on hydrogen storage, resulting in a significant decline in the total energy conversion cost (T_cost). By optimizing energy utilization through DRP, the EH achieves a more cost-efficient operation while simultaneously minimizing unnecessary energy transformations. This highlights the potential of DRP to enhance both economic and operational performance within the EH. Compared to scenarios without DRP, both electricity demand during peak hours and natural gas consumption are significantly lower, leading to substantial cost savings. Additionally, the hydrogen storage system contributes approximately 4.83% to the total gas supply within the EH, enhancing its operational efficiency.

A comprehensive sensitivity analysis was conducted to examine the influence of DRP on operational costs and emissions. The findings, depicted in Figure 14, reveal a clear relationship between the percentage of load participation in DRPs and the corresponding reductions in EH costs and emissions. These results underscore the critical role of DRP in optimizing energy management and achieving cost-effective, environmentally sustainable operations in the EH.

[figure(s) omitted; refer to PDF]

Figure 11, indicates that as the percentage of load participation in the DRP increases, the reduction in cost and emissions steadily rises; however, this increase gradually follows a diminishing rate. In the initial stages of participation (0%–10%), the slopes of both curves are steeper, reflecting a more significant impact of the DRP on reducing costs and emissions. Nevertheless, the cost reduction curve (blue) is higher than the emissions reduction curve (red), highlighting the greater effectiveness of the DRP in reducing costs. In the range of 20% to 30% participation, the curves approach saturation, indicating that further increases in load participation have a diminished effect on cost and emissions reduction. Ultimately, the results suggest that the DRP is an effective tool for reducing costs and emissions, but optimizing the level of load participation is essential to achieve maximum efficiency.

• Considering DRP and uncertainties

In the third phase of simulation results analysis, the energy management of the EH integrates both uncertainties and DRP. This approach aims to reflect a more realistic operational environment, incorporating the inherent variability of energy resources and loads alongside the flexibility introduced by DRP. The outcomes of this advanced simulation phase are illustrated in Figure 15, showcasing the pareto fronts derived from the multi-objective optimization algorithms. These Pareto fronts highlight the trade-offs between competing objectives, such as cost minimization and emission reduction, under the combined influence of uncertainties and DRP. The results emphasize the robustness and adaptability of the proposed optimization strategies in achieving an efficient and balanced energy management solution.

[figure(s) omitted; refer to PDF]

Table 4 provides a detailed summary of the optimal solutions selected from the Pareto fronts generated by the multi-objective optimization algorithms. These solutions are accompanied by a breakdown of the EH operating cost components, offering insights into the cost distribution across various elements, such as electricity procurement, gas usage, and energy conversion processes. This detailed analysis highlights the efficiency of the proposed optimization methods in minimizing overall operational costs while addressing multiple objectives effectively.

Table 4

The best results in the third section of the simulation.

Optimization algorithm	$E_{\text{Cost}}¢$	$H_{\text{Cost}}¢$	$T_{\text{Cost}}¢$	$\text{Total\hspace{0.17em}cost}¢$	$\text{Emission}\text{kg}$
MMOGOA	975	723	314	2012	13.68
MOPSO	1041	708	294	2041	14.18
NSGA-II	1014	754	267	2035	14.07
MOWOA	985	764	277	2026	13.84

Considering uncertainties in the EH EMS significantly impacted the performance metrics of the EH, leading to an observable rise in both operational costs and emissions. Specifically, the integration of uncertainties caused operating costs to escalate by 10%, while emissions experienced a 4.38% increase compared to scenarios where uncertainties were ignored. Despite these challenges, the proposed MMOGOA algorithm proved to be the most effective among the tested optimization methods. Under the MMOGOA framework, the EH achieved the lowest operational cost of 2012 cents and an emission level of 13.68 kg. In comparison, the MOPSO algorithm yielded an operating cost of 2041 cents with emissions of 14.18 kg, while the NSGAII algorithm recorded costs of 2035 cents and emissions of 14.07 kg. Similarly, the MOWOA algorithm, though slightly better than MOPSO and NSGAII, still lagged behind MMOGOA with an operating cost of 2026 cents and emissions of 13.87 kg. The increased costs and emissions are attributed to the system’s need to account for unpredictable variations in energy demand and supply, which complicates energy allocation and storage management. Figure 16 highlights the hourly charging patterns of storage units, showcasing how uncertainties influence the utilization of these resources. This demonstrates the robustness of the MMOGOA algorithm in optimizing the EHs performance despite the added complexity of uncertainty, making it a reliable choice for real-world EMS applications.

[figure(s) omitted; refer to PDF]

The results from this phase of the simulations reveal that the storage capacity was used more efficiently compared to previous scenarios. This improvement is largely due to the reduced influence of uncertainties on the energy management of the EH. A key factor contributing to the lower operating costs and emissions observed with the EMS model optimized by the MMOGOA algorithm is its effective management of storage charge levels. As illustrated in Figure 13, there was a significant reduction in the amount of electricity and natural gas purchased from the grid, indicating more efficient utilization of on-site resources. In addition, to further minimize the effects of uncertainties and thereby decrease costs and emissions, hydrogen storage technology was more heavily employed. Also, the hydrogen storage system provided around 4.23% of the total natural gas demand in the EH, showcasing its critical role in balancing energy supply and demand while reducing reliance on external energy sources. This highlights the effectiveness of the MMOGOA-optimized EMS in enhancing energy storage management and promoting a more sustainable and cost-effective EH operation.

5. Conclusion and Future Study

This paper proposes a probabilistic bi-objective EMS model for an EH integrating renewable energy sources (PV systems and WT) with the main grid. The hub configuration comprises a boiler, CHP unit, thermal/EV parking facilities, and hydrogen storage technology. The EMS model simultaneously minimizes operational costs and emissions while incorporating DRPs and addressing uncertainties in energy demand and generation. For optimization, we develop a MMOGOA. The algorithm’s effectiveness is first verified through ZDT benchmark functions, demonstrating superior performance compared to conventional multi-objective optimization methods. The EMS model is evaluated through three scenarios: (1) baseline operation without uncertainties or DRPs, (2) system with DRPs, and (3) system with both DRPs and uncertainties. Simulation results indicate MMOGOA’s consistent superiority across all scenarios. The baseline scenario achieved minimal operating costs and emissions. DRP implementation yielded 18.87% cost reduction and 14.62% emission decrease. While uncertainty introduction increased costs by 10% and emissions by 4.38%, MMOGOA maintained optimal performance. These findings highlight the importance of DRP optimization and uncertainty management in sustainable energy systems. Future research is suggested to explore the resilience and reliability of the EH to further improve its robustness and ensure efficient operation in the face of varying environmental and operational conditions.

Author Contributions

Shahriar Karimian: formal analysis, software, writing – original draft. Majid Moazzami: writing – review and editing, methodology, investigation. Bahador Fani: writing – review and editing, supervision, methodology, investigation. Ghazanfar Shahgholian: writing – review and editing, supervision.

Acknowledgments

The authors have nothing to report.

References

[1] S. Zhong, X. Wang, H. Wu, Y. He, B. Xu, M. Ding, "Energy Hub Management for Integrated Energy Systems: A Multi-Objective Optimization Control Strategy Based on Distributed Output and Energy Conversion Characteristics," Energy, vol. 306, 2024.

[2] Z. Yadollahi, R. Gharibi, R. Dashti, A. T. Jahromi, "Optimal Energy Management of Energy Hub: A Reinforcement Learning Approach," Sustainable Cities and Society, vol. 102, 2024.

[3] K. Esapour, F. Moazzen, M. Karimi, M. Dabbaghjamanesh, A. Kavousi-Fard, "A Novel Energy Management Framework Incorporating Multi-Carrier Energy Hub for Smart City," IET Generation, Transmission and Distribution, vol. 17 no. 3, pp. 655-666, 2023.

[4] R. Parekh, B. Sedhom, S. Padmanaban, A. A. Eladl, "A Review of IoT-Enabled Smart Energy Hub Systems: Rising, Applications, Challenges, and Future Prospects," Renewable Energy Focus, vol. 51, 2024.

[5] S. M. Ezzati, F. Faghihi, H. M. Shourkaei, S. B. Mozafari, S. Soleymani, "Optimum Operation of Multi-Energy Carriers in the Context of an Energy Hub considering a Wind Generator Based on Linear Programming," Journal of Renewable and Sustainable Energy, vol. 10 no. 1, 2018.

[6] H. Ahmadisedigh, L. Gosselin, "Combined Heating and Cooling Networks With Part-Load Efficiency Curves: Optimization Based on Energy Hub Concept," Applied Energy, vol. 307, 2022.

[7] J. Bai, W. Wei, L. Chen, S. Mei, "Rolling-Horizon Dispatch of Advanced Adiabatic Compressed Air Energy Storage Based Energy Hub via Data-Driven Stochastic Dynamic Programming," Energy Conversion and Management, vol. 243, 2021.

[8] S. Rech, S. Casarin, C. S. Silva, A. Lazzaretto, "University Campus and Surrounding Residential Complexes as Energy-Hub: A MILP Optimization Approach for a Smart Exchange of Solar Energy," Energies, vol. 13 no. 11, 2020.

[9] A. Moser, D. Muschick, M. Gölles, "A MILP-Based Modular Energy Management System for Urban Multi-Energy Systems: Performance and Sensitivity Analysis," Applied Energy, vol. 261, 2020.

[10] A. Mansour-Saatloo, M. Agabalaye-Rahvar, M. A. Mirzaei, B. Mohammadi-Ivatloo, M. Abapour, K. Zare, "Robust Scheduling of Hydrogen Based Smart Micro Energy Hub With Integrated Demand Response," Journal of Cleaner Production, vol. 267, 2020.

[11] P. Li, W. Sheng, Q. Duan, Z. Li, C. Zhu, X. Zhang, "A Lyapunov Optimization-Based Energy Management Strategy for Energy Hub With Energy Router," IEEE Transactions on Smart Grid, vol. 11 no. 6, pp. 4860-4870, 2020.

[12] H. A. Honarmand, A. G. Shamim, H. Meyar-Naimi, "A Robust Optimization Framework for Energy Hub Operation considering Different Time Resolutions: A Real Case Study," Sustainable Energy, Grids and Networks, vol. 28, 2021.

[13] M. Ghorab, "Energy Hubs Optimization for Smart Energy Network System to Minimize Economic and Environmental Impact at Canadian Community," Applied Thermal Engineering, vol. 151, pp. 214-230, 2019.

[14] U. Güvenç, B. Özkaya, H. Bakir, S. Duman, O. Bingöl, "Energy Hub Economic Dispatch by Symbiotic Organisms Search Algorithm," pp. 375-385, .

[15] T. Liu, D. Zhang, H. Dai, T. Wu, "Intelligent Modeling and Optimization for Smart Energy Hub," IEEE Transactions on Industrial Electronics, vol. 66 no. 12, pp. 9898-9908, 2019.

[16] S. M. Moghaddas-Tafreshi, M. Jafari, S. Mohseni, S. Kelly, "Optimal Operation of an Energy Hub considering the Uncertainty Associated With the Power Consumption of Plug-in Hybrid Electric Vehicles Using Information Gap Decision Theory," International Journal of Electrical Power and Energy Systems, vol. 112, pp. 92-108, 2019.

[17] J. Zou, X. Yang, Z. Liu, J. Liu, L. Zhang, J. Zheng, "Multiobjective Bilevel Optimization Algorithm Based on Preference Selection to Solve Energy Hub System Planning Problems," Energy, vol. 232, 2021.

[18] X. Zeng, Z. Wu, Y. Wang, D. Song, G. Li, "Multi-Mode Energy Management Strategy for Hydraulic Hub-Motor Auxiliary System Based on Improved Global Optimization Algorithm," Science China Technological Sciences, vol. 63, pp. 2082-2097, 2020.

[19] H. Sadeghi, M. Rashidinejad, M. Moeini-Aghtaie, A. Abdollahi, "The Energy Hub: An Extensive Survey on the State-of-the-Art," Applied Thermal Engineering, vol. 161, 2019.

[20] D. Rakipour, H. Barati, "Probabilistic Optimization in Operation of Energy Hub With Participation of Renewable Energy Resources and Demand Response," Energy, vol. 173, pp. 384-399, 2019.

[21] X. Wang, Y. Liu, C. Liu, J. Liu, "Coordinating Energy Management for Multiple Energy Hubs: From a Transaction Perspective," International Journal of Electrical Power and Energy Systems, vol. 121, 2020.

[22] W. Gan, M. Yan, W. Yao, J. Wen, "Peer to Peer Transactive Energy for Multiple Energy Hub With the Penetration of High-Level Renewable Energy," Applied Energy, vol. 295, 2021.

[23] M. A. Mirzaei, M. Zare Oskouei, B. Mohammadi-Ivatloo, "Integrated Energy Hub System Based on Power-to-Gas and Compressed Air Energy Storage Technologies in the Presence of Multiple Shiftable Loads," IET Generation, Transmission and Distribution, vol. 14 no. 13, pp. 2510-2519, 2020.

[24] M. Jadidbonab, A. Dolatabadi, B. Mohammadi-Ivatloo, M. Abapour, S. Asadi, "Risk-Constrained Energy Management of PV Integrated Smart Energy Hub in the Presence of Demand Response Program and Compressed Air Energy Storage," IET Renewable Power Generation, vol. 13 no. 6, pp. 998-1008, 2019.

[25] M. Aghamohamadi, N. Amjady, A. Attarha, "A Linearized Energy Hub Operation Model at the Presence of Uncertainties: An Adaptive Robust Solution Approach," International Transactions on Electrical Energy Systems, vol. 30 no. 3, 2020.

[26] Y. Mu, C. Wang, M. Sun, W. He, W. Wei, "CVaR-Based Operation Optimization Method of Community Integrated Energy System Considering the Uncertainty of Integrated Demand Response," Energy Reports, vol. 8, pp. 1216-1223, 2022.

[27] J. Faraji, H. Hashemi-Dezaki, A. Ketabi, "Stochastic Operation and Scheduling of Energy Hub considering Renewable Energy Sources’ Uncertainty and N-1 Contingency," Sustainable Cities and Society, vol. 65, 2021.

[28] M. Jiang, T. Liu, H. H. Goh, "Estimation of Operation Cost of Residential Multiple Energy System considering Uncertainty of Loads and Renewable Energies," IEEE Access, vol. 9, pp. 4874-4885, 2020.

[29] Q. Cai, X. Luo, P. Wang, C. Gao, P. Zhao, "Hybrid Model-Driven and Data-Driven Control Method Based on Machine Learning Algorithm in Energy Hub and Application," Applied Energy, vol. 305, 2022.

[30] S. Habib, M. A. Kamarposhti, H. Shokouhandeh, I. Colak, "Economic Dispatch Optimization considering Operation Cost and Environmental Constraints Using the HBMO Method," Energy Reports, vol. 10, pp. 1718-1725, 2023.

[31] H. Shokouhandeh, M. Ahmadi Kamarposhti, F. Asghari, I. Colak, K. Eguchi, "Distributed Generation Management in Smart Grid With the Participation of Electric Vehicles with Respect to the Vehicle Owners’ Opinion by Using the Imperialist Competitive Algorithm," Sustainability, vol. 14 no. 8, 2022.

[32] M. Ghadi-Sahebi, R. Ebrahimi, V. Parvin-darabad, "Optimal Probabilistic Operation Management of Smart Parking Lot and Renewable Sources in Microgrid to Reduce Cost and Improve System Reliability considering Demand Response Program," International Transactions on Electrical Energy Systems, vol. 31 no. 12, 2021.

[33] S. Zhan, Y. Zhou, D. Feng, H. Wang, C. Fang, "Energy Optimization Dispatch Based on Two-Stage and Multi-Objective Comparative Analysis for PV and BESS Integrated Fast Charging Stations With V2G," IET Renewable Power Generation, vol. 18, 2024.

[34] B. Zhou, M. Wen, "A Multipopulation Grasshopper Algorithm for Multiobjective Material Distribution Scheduling Problem of Assembly Line With Flexible Line Inventory," Robotic Intelligence and Automation, vol. 44 no. 5, pp. 724-745, 2024.

[35] K. Kalita, J. V. N. Ramesh, L. Cepova, S. B. Pandya, P. Jangir, L. Abualigah, "Multi-Objective Exponential Distribution Optimizer (MOEDO): A Novel Math-Inspired Multi-Objective Algorithm for Global Optimization and Real-World Engineering Design Problems," Scientific Reports, vol. 14 no. 1,DOI: 10.1038/s41598-024-52083-7, 2024.

[36] Y. Tian, S. Shao, G. Xie, X. Zhang, "A Multi-Granularity Clustering Based Evolutionary Algorithm for Large-Scale Sparse Multi-Objective Optimization," Swarm and Evolutionary Computation, vol. 84, 2024.

[37] Q. Shang, M. Tan, R. Hu, Y. Huang, B. Qian, L. Feng, "A Multi-Stage Competitive Swarm Optimization Algorithm for Solving Large-Scale Multi-Objective Optimization Problems," Expert Systems with Applications, vol. 260, 2025.

[38] Y. Han, H. Peng, C. Mei, "Multi-Strategy Multi-Objective Differential Evolutionary Algorithm With Reinforcement Learning," Knowledge-Based Systems, vol. 277, 2023.

[39] M. Akbarpour, A. Abdollahi, M. Abdali, M. Vosoogh, "Economic-Environmental Optimization of an Energy Hub with DRPs and Power-to-Gas Technology under Uncertainties With Introducing a Novel Optimization Algorithm," Journal of Energy Management and Technology, vol. 8 no. 3, pp. 224-238, 2024.