Operation and Planning of Energy Hubs Under Uncertainty—A Review of Mathematical Optimization Approaches

Co-designing energy systems across multiple energy carriers is increasingly attracting attention of researchers and policy makers, since it is a prominent means of increasing the overall efficiency of the energy sector. Special attention is attributed to the so-called energy hubs, i.e., clusters of energy communities featuring electricity, gas, heat, hydrogen, and also water generation and consumption facilities. Managing an energy hub entails dealing with multiple sources of uncertainty, such as renewable generation, energy demands, wholesale market prices, etc. Such uncertainties call for sophisticated decision-making techniques, with mathematical optimization being the predominant family of decision-making methods proposed in the literature of recent years. In this paper, we summarize, review, and categorize research studies that have applied mathematical optimization approaches towards making operational and planning decisions for energy hubs. Relevant methods include robust optimization, information gap decision theory, stochastic programming, and chance-constrained optimization. The results of the review indicate the increasing adoption of robust and, more recently, hybrid methods to deal with the multi-dimensional uncertainties of energy hubs.

The associate editor coordinating the review of this manuscript and approving it for publication was Easter Selvan Suviseshamuthu .

I. INTRODUCTION
The joint co-optimization of different energy carriers, traditionally operated and planned separately, forms a prominent endeavor towards efficient energy and infrastructure utilization [1], since co-optimizing decisions across multiple energy carriers, features a synergistic effect [2]. Hence, the Energy Hub (EH), an energy management unit that encompasses multiple energy carriers and technologies, is envisioned as a new concept that enhances the integrated system's efficiency [3]. The operational decisions of an EH concern the quantities dispatched for each of the EH's energy carriers (and the necessary conversions across carriers) in operational time; the operational problem generally refers to satisfying the hub's energy demands (respecting energy flow constraints) in the most economically efficient way. The planning decisions on the other hand, concern the decisions of investing in infrastructure expansions for the EH relating to techno-economic aspects of the EH's future operation.
The spectrum of relevant technologies includes combined heat and power (CHP) technology, boiler, heat exchangers, different types of energy storage systems, power-to-X technologies, and more [4]. An architectural overview of an EH, and its featured devices, inputs, and outputs is presented in Fig 1. A. LITERATURE REVIEW Some of the well-known methods for decision making under uncertainty in energy systems, including probabilistic methods, robust optimization, information gap decision theory, possibilistic methods, and interval based analysis, have been categorized and explained in [5]. A simple formulation for each of these methods, as well as some energy related applications for each of them have been provided in this reference. In [6] probabilistic methods, possibilistic methods, and their combinations for the situations when the problem includes both probabilistic and possibilistics uncertain variables have been presented in details. In this reference numerical probabilistic approaches, including sequential, non-sequential and pseudo-sequential Monte Carlo simulation; as well as, analytical probabilistic approaches including convolution method, cumulant method, Gram-Charlier method, Edgeworth expansion, Taylor series, Cornish-Fisher expansion, and point estimation method have been simply formulated. A comprehensive review on optimization of energy systems under uncertainty has been presented by [7]. In this reference energy system model have been classified and different optimization methodologies, including stochastic programming, robust optimization, fuzzy programming, interval method and hybrid models have been reviewed.
Previous review papers, specifically in the area of EHs, have investigated the EH research from different points of view. Article [8] reviews different models and concepts used in the EH literature, identifying and discussing their challenges, strengths, and weaknesses. Models, concepts, and technologies of EHs in the commercial, industrial, domestic, and agricultural sectors are briefly reviewed in [9]. The authors in [10] review the topologies, design methods, and input-output relations of EHs, as well as design requirements such as cost, robustness, pollution, flexibility, and resilience. The study also briefly reviewes the mathematical models of EHs. Article [11] reviews articles about multi-carrier energy systems, published from 2007 to 2017. The reviewed papers have been summarized based on their environmental and economic considerations, applications, operation and planning, as well as their modeling approaches. Article [12] presents the EH as a solution concept to energy-related challenges and proceeds to present a literature review on the operation, planning and expansion planning problems of EHs, summarizing the types of objective functions, planning horizon length and granularity, and different forms of energy carriers considered. The performance of existing energy storage systems and multi-energy system designs have been reviewed and compared by [13]. The state-of-the-art research on EHs is linked to the concept of energy-positive neighborhoods in [14]. The authors in [15] present a bibliometric analysis and a systematic review based on a qualitative synthesis of the EH research studies.

B. MOTIVATIONS AND CONTRIBUTIONS
As it was clear from the literature review section, the conducted research in the area of energy hubs have been investigated and reviewed from various aspects. Moreover, approaches for decision making under uncertainty in power systems have been reviewed carefully. To the best of the authors' knowledge, no comprehensive review has been previously conducted on the application of mathematical optimization approaches towards solving operational and planning decision problems of EHs under uncertainty. This paper investigates the relevant approaches, including stochastic programming (SP), robust optimization (RO), information gap decision theory (IGDT), chance constrained (CC) optimization etc, and identifies the various frameworks proposed  VOLUME 11, 2023 in the literature. The main purposes of this work are 1) to present the relevant modeling frameworks, methods and uncertainties in a systematic way, discussing each method and categorizing the relevant literature, while performing a critical review on the advantages and disadvantages of each approach ; and, 2) to extract the research trend in the ares of energy hubs. In other words, it has been intended to show which mathematical optimization methods were popular to date and why, and which mathematical optimization methods are going to become more popular among researchers in the energy hub domain and why. In addition, we also discuss aspects related to the design of EHs' inputs, outputs, and devices.

7210
In summary, the contributions of this paper are highlighted as follows: • presenting a comprehensive review on the application of mathematical optimization approaches in the operation and planning of EHs under uncertainty.
• categorizing and summarizing the research conducted in the field of operation and planning of EHs from the mathematical optimization's point of view.
• revealing research trends in this field and presenting the most possible future direction from the mathematical optimization's point of view.

C. PAPER ORGANIZATION
The rest of the paper is organized as follows. The methodology of the literature review is given in Section 2. The mathematical modeling and the relevant optimization methods applied to EHs are presented in Section 3. Hybrid methods are presented in Section 4. Section 5 presents a summary of the related literature and discusses the emerging research directions. Finally, the conclusions are given in Section 6.

II. METHODOLOGY
In this section, the justification of the optimization methods is presented. The justification considers application of only the mathematical optimization techniques. After that, the flow of adequate papers selection based on research keywords (considering EH and mathematical optimization techniques), type of database, date and language of publication is presented.

A. SELECTION OF OPTIMIZATION METHODS
Optimizing an EH's operation and planning decisions can be formulated as an optimization problem, involving one or multiple objective function(s), subject to several technical and non-technical constraints. In this review, we focus on mathematical optimization approaches that provide formal guarantees regarding global system optimality and constraint satisfaction. It bears mentioning that a significant body of the EH literature adopts approaches based on evolutionary or meta-heuristic algorithms. These approaches are able to manage non-convexities in the objective functions and constraints, but do not meet our criterion of global optimality and constraint satisfaction guarantees. Hence, they are briefly reviewed in Appendix A for the interested reader, and are not further discussed in this paper. In addition, Machine Learning approaches have received great attention recently for the operation and planning of EHs. Since they are in a different category of study (it is not in the scope of mathematical optimization), they have been added to Appendix B for interested readers. The mathematical optimization techniques qualifying for detailed review, are categorized in Figure 2.

B. METHODOLOGY OF ARTICLE SELECTION
The methodology of the review investigation was based on a structured selection with the Scopus database used as the main resource, encompassing all WoS indexed journals, all IEEE and IoP conferences, and other verified venues. Our search method used the keywords ''energy hub'' or ''multicarrier energy systems'', in combination with the names of the relevant optimization techniques (SP, IGDT, RO, and CC). The considered articles were accessed before 31.12.2021 and their language was English. The summary of articles selection criteria is presented in Figure 3.

III. MATHEMATICAL MODELING AND OPTIMIZATION TECHNIQUES
A. STOCHASTIC PROGRAMMING Stochastic programming (SP) considers the optimization of the expected value of an objective optimized over N (i.e. two or more) decision stages [16]. The uncertainty over future realizations of the problem's parameters is captured using a number e of realization scenarios, with each scenario bearing a probability π e , where e e=1 π e = 1. The method uses an optimal-in-hindsight decision x e for each scenario e and constraints scenarios with common history up to stage k to have the same up-to-k decisions. Such constraints are called non-anticipativity constraints and make sure that, at any decision stage k, the algorithm only utilizes information that has been revealed up and before k.
The general form of an N-stage SP optimization problem is formulated in (1). The objective function f (·), typically relates to minimizing the total economic cost of the system or maximizing the RESs integration. Finally, the equality and inequality constraints g, h refer to the system's operational and modeling constraints.
Subject to: The main challenge of SP models is to handle the complexity arising from multiple scenarios and decision stages. Ideally, a small set of scenarios would need to be identified, that captures enough information to achieve a satisfactory level of efficiency while keeping the optimization problem manageable within realistic timeframes and computational VOLUME 11, 2023

resources. A review of related cases and approaches follows.
A two-stage stochastic programming approach is considered in [17] towards optimizing the operational decisions of a multi-carrier energy building. The authors in [18], optimize an EH's planning decisions, while [19] uses SP to make profit-maximizing decisions on behalf of an EH managing entity. Generalizing to a cluster of multiple energy hubs, [20] used two-stage SP to assess the system's economic efficiency.
The work in [21] focuses on small scale EHs, namely residential buildings, and controls the building's lighting, ventilation, cooling and heating systems as well as charge and discharge of electric vehicles. In [22] a real case study of a building-level smart EH is considered, minimizing the weighted sum of the EH's energy bill on one hand, and an emissions' penalization factor on the other.
The application of SP techniques in a real microgrid with demand response capabilities is presented in [23].
Minimizing the expected stacked operational costs of the electricity and gas networks is the objective adopted in [24]. A stochastic-interval optimization was applied in [25] towards achieving cost reduction for an EH with flexible loads. In [26], the authors present a comprehensive uncertainty modeling framework and scenario generation method towards representing the uncertainty of demands (for electricity and heating), wind speed, solar irradiance, and prices of energy carriers including electricity and natural gas. The objective adopted was the maximization of the EH's operator's profit.
The study [27] formulates the scheduling problem of an EH as a stochastic program. The EH features P2G storage, a CHP unit, WT power, boiler, electrical and thermal storage, as well as demand response (DR) capabilities. The formulated objective is to meet all demands (electrical, heat, and gas) at the least possible cost. The problem is cast as a mixed-integer linear programming problem. The coupling point between the electricity and gas carriers is the P2G storage system, since it can convert energy from one carrier to the other through hydrogen. In [28], the authors formulate a stochastic security-constrained unit commitment problem, considering the coordinated operation of price-based DR and HES system in the presence of wind generation uncertainty. The authors also consider price responsive loads and demonstrate their case study on a 6-bus and on a 24-bus system. The results show the effects of joint consideration of a HES system and a price-based DR program. In [29], the authors present a mathematical formulation for the optimal planning problem of an EH considering operational constraints and multiple objectives related to investment and operational costs, reliability, and emissions. The planning decisions are made under uncertainty of wind generation, electricity prices, and electricity demands. The EH features a transformer, a CHP, G2HE, and HESS.
In [30], system featuring WT, EES and HESS, and electrical and thermal DR programs, is optimized. The system's uncertainties include demands, market prices, and wind speed. An additional degree of freedom is enabled by considering the existence of a market for allocating heat demand. In [31], SP is used towards evaluating the impact of ESSs on a local multi-energy system. The linepack model of gas pipelines, is leveraged towards enabling the pipelines to act as a buffer, i.e. energy storage. The problem is cast as a mixed-integer linear program, solved by the CPLEX solver in the GAMS environment. In [32], a SP formulation was leveraged towards dealing with several sources of uncertainty, stemming from day-ahead market prices, real-time market prices, and WT generation. In [33], a SP formulation was deployed similarly to deal with the uncertainties of WT and PV generation. In [34], a networked version of EH operation has been taken into account with multiple uncertainties of PV generation, electricity prices as well as demands.
In conclusion of this section, we highlight that EH management systems entail various uncertain parameters related to each energy carrier, namely energy (electricity, gas, heat) demands, wholesale (electricity, gas) prices, generation output (from wind and solar panels for electricity) and more. The EH receives the uncertainty realization as input and decides upon its control actions over the available controllable resources (generally presented in Fig. 2). The resulting actions, i.e. the outputs, concern the control of energy flows across the different energy carriers. The inputs, controllable resources, and outputs discussed vary among different research studies. Table 1, summarizes the model (in terms of inputs, resources, output) adopted in each paper that uses SP as its optimization framework.

B. INFORMATION GAP DECISION THEORY
Information gap decision theory (IGDT), which was introduced by Ben-Haim [57], [58], is a non-probabilistic decision-making method to maximize the system's robustness under severely uncertain circumstances. In this method, uncertainties are represented using expected values and a sensitivity analysis is performed to determine the effect of expected value perturbations on the optimization results [57]. To date, several uncertainty models, such as envelope-bound models, energy-bound models, Minkowski-norm models, slope-bound models, and Fourier-bound models [59], have been developed to model the uncertainty parameters. The most well-known uncertainty model is the envelope-bound model, which is presented as follows: where e is the uncertain parameter andẽ is its expected value. Let α denote the uncertainty horizon which determines the upper and lower bounds of the fractional deviations of the uncertain parameter comparing to its expected value by (4). The performance level of the decision can be expressed in terms of profit, cost, or other relevant functions [60]. Let f (x; e) denote the objective function that has to be maxi-7214 VOLUME 11, 2023 mized, where x is a decision variable, and e is an uncertain parameter. In conventional optimization approaches, f (x; e) is maximized. In the IGDT approach, we are after maximizing robustness, as in: where f critical is the critical (i.e. worst acceptable) performance level of the decision. In the IGDT method, the decision making variables, x, are selected such that the safe interval of the undetermined parameter, α, is maximized considering the fact that the performance level has to be more than its critical value. To do this, if the maximization objective is the profit, in this approach instead of determining decision making variables to maximize the profit, they are selected in a way that the horizon of uncertainty in which the minimum acceptable amount of profit is assured, is maximized. Works using IGDT as the mathematical optimization tool to make decision under uncertainty in EHs are summarized in Table 2. While the basic IGDT has the ability to manage one uncertain parameter, the method has been extended to manage multiple sources of uncertainty. From this point of view, the works summarized in Table 2 can be categorized into three classes depending on the number of uncertain parameters that each one considers.
In the first category, there is only one uncertain parameter. Article [70] used IGDT to model the uncertainty of plug-in hybrid electric vehicles' power consumption during trips under risk-averse and risk-seeking strategies. The electricity price uncertainty [72], the electric load uncertainty [74], and wind power generation uncertainty [68], [69] have also been modeled using IGDT. Thermal demand, electricity demand, and energy production of renewable resources are (separately -one at a time) considered as uncertain parameters in [75].
The second and third category concern handling multiple uncertain parameters. In the second category, although the uncertain parameters are multiple, they are all modeled by one uncertainty horizon. Parameters including the market's energy price, loads, and renewable energy sources' output power are assumed to be uncertain in [67]. In article [40], the output of the renewable generation units and the electricity loads are considered as the uncertain parameters. Real-time electricity market prices and wind turbine generation are the uncertain parameters that have been handled in [73]. The main shortcoming of this approach is that it cannot handle different scales of uncertainty for different parameters. In this context, consider the operator of an EH, needing to make a decision under the uncertainties of the wholesale market's price and the uncertainty of wind's power output. The forecasting error of the wholesale market's price and wind's power output are not generally equal. Consequently, using a single uncertainty horizon cannot guarantee finding the most robust solution since the variations of the uncertain parameters around their expected values are not the same.
Studies in the third category, deal with multiple sources of uncertainty. The authors in [63] consider uncertain WT and PV generation, energy demands (including electrical, thermal, and cooling), and day-ahead electricity prices. The uncertainty horizon of electricity prices was maximized using an enhanced version of the conventional ϵ-constraint method, named AUGMECON [76]. Towards this goal, the uncertainty horizons of demands and renewable generations were calculated by a lexicographic optimization technique [77]. In this approach, the management of multiple uncertain parameters increases the computation time dramatically; consequently, a trade-off between the computation time and the Pareto set density emerges. In article [66], the uncertain parameters are the electric demands, cooling demands, heat demands, battery charging station's demands, PV's power output, wind's power output and electricity prices. The authors proposed maximizing the minimum uncertainty radius of the uncertain parameters as a way of handling multiple sources of uncertainty. However, this re-introduces the same drawbacks as in the studies of the second category discussed above. In [71] a multi-objective IGDT approach has been proposed to model the uncertainties of demand, wind's power output, and PV's power output. In this approach, the uncertainty horizons of the three uncertain parameters are considered as the objective functions of the multi-objective optimization method. To do this, a modified version of the directed search domain II method [78] has been used as the multi-objective optimization method.

C. ROBUST OPTIMIZATION
In RO, the system designer uses an uncertainty set to define the possible region of uncertainty realization, A RO approach aims at determining a solution to an optimization problem that is valid for any realization of the uncertain parameter within the predefined uncertainty set. The obtained solution gives the optimal solution for the worst-case realization of the uncertainties [79]. A generic formulation of the RO approach is represented as follows: e ≤ e max : δ ≥ 0 Eq. (7) declares a general cost function f (·), with e and x being the uncertain parameters and decision variables, respectively, while Eq. (8) represents the relevant inequality constraints. The goal is to minimize the cost function (outer minimization) under the worst-case uncertainty realization (inner maximization). The maximum deviation of the uncertain parameter e is given by (9). Eq. (10) shows that the uncertain parameter can deviate up to the maximum value of e max , while ζ and δ denote the dual variables of the respective constraints. Using duality theory, the inner maximization VOLUME 11, 2023  problem of the uncertainty deviation can be written as min ζ,δ e max δ + eζ (11) subject to (8) and, where (12) and (13) represent the duals' constraints. This is a general method towards rendering the original RO problem solvable by commercial solvers. The reviewed papers using RO are divided into three categories; 1) Classic RO, 2) Adaptive robust optimization (ARO) and 3) Distributionally robust optimization (DRO). In what follows, each category is seperately discussed.

1) CLASSIC RO
This is the most widely used RO approach. The static RO normally addresses the mathematical optimization to reach the worst-case realization of uncertainties, disregarding binary and integer variables. Therefore, it has attracted the attention of researchers due to its simplicity compared to the ARO and DRO. The typical objective of this category concerns economical aspects, i.e., minimizing/maximizing the total cost/profit of the operation and planning of EHs. Several articles have only taken into account the uncertainty of the renewable energy sources, i.e., wind power and PV generation [104], [115], [124]. In addition, [80] utilized classic RO to deal with the uncertainty of PV generation in a real case study. Other studies have only focused on the uncertainty of the electricity market prices e.g. [122], [125], [130], [131]. Article [99] deployed RO to cope with the uncertainty of electricity prices in a decentralized model, where the distribution network and EH are managed separately using the alternating direction method of multipliers (ADMM). Similarly, the article [91] has only considered the market prices' uncertainty, while the main goal of the paper is to include and assess the application of hydrogen energy. In [84], the authors focus an designing a CHP in an EH, again under prices' uncertainty. Reference [4] considered the operation of a P2X unit under uncertain electricity market prices.
Uncertainty on demands is another source of attention for studies of this category. Articles [79], [106] have addressed the uncertainties of various kinds of demands, e.g., electricity, heating, and cooling. For instance, in [106] a robustness constraint is formed, and the corresponding robust optimization model is equivalently transformed into a mixed-integer programming model. Only one of the mentioned papers in TABLE 3. RO application to EH issues. VOLUME 11, 2023 this category, namely [81], has considered the uncertainty of the driving pattern of electric vehicles (EVs).
The rest of the papers in this category have studied a combination of uncertainties. Article [97] used RO to deal with the uncertainties of wind and PV generations, as well as electricity market prices. The wind power and PV have been considered by one equation, only in the constraints. Thus, the worst-case realization was only developed for the price uncertainty. The uncertainties of electricity market prices, wind, and demands have been studied by the RO also in [86] and [117], with extended affine arithmetic and RO being used respectively. The solution to the RO problem has been developed through an evolutionary algorithm in [85], where the authors consider a combination of the mentioned uncertainty sources.

2) ADAPTIVE ROBUST OPTIMIZATION
In an ARO the optimization variable vector, which may also contain both discrete and continuous variables, is made once the uncertain parameters vector is realized. In this case, this reaction after the realization of uncertain parameters is called recourse which leads to using ARO. On this occasion, the RO framework requires extension to handle the computational complexity, pointing to ARO methods. Generally, the papers in this category have used iterative methods of column and constraint generation (C&CG) algorithm or Benders decomposition. Articles [89], [101], [107], [109], [112], [113], [126] have used ARO to handle the binary variables stemming from the start-up and shut-down generator constraints in the economic dispatch and unit commitment problem. Regarding uncertainty, [109] copes with the uncertainty of wind power generation, while [101] has taken into account both renewable energy and price uncertainties. The unit commitment constraints have led the authors in [82] to use ARO in the presence of uncertainties of wind power, loads and electricity prices, where all the uncertainties have been addressed by one interval of deviation. In [118] a two-stage model was proposed for an ARO problem. The planning of an ESS was investigated in the first stage and the uncertainty of the loads and market prices were investigated in the second stage. The Benders decomposition method is utilized to handle the computational complexity. Similarly, in [105], [111], and [120] a two-stage planning and operation method was suggested, where the planning of an EH is investigated in the first stage and the operation problem is solved in the second stage. The worst-case realization of the uncertainty of real-time electricity market prices was modeled in the second stage of an ARO in [1]. Finally, [94], [98] are the papers in this second category that use C&CG to deal with the uncertainties.

3) DISTRIBUTIONALLY ROBUST OPTIMIZATION
Stochastic programming (discussed in section III-A) assumes a known probability distribution function describing the uncertainty realization. However, in many cases, it is difficult or impossible to know a prior probability distribution function. In the DRO formulation, the solution is taken with respect to the worst-case probability distribution of an uncertain parameter, within a certain set of possible distributions. DRO is based on an ambiguity set and is less conservative compared to classic RO and ARO. In addition, there are various kinds of methods to solve DRO problems. Two main of them are moment-based and Wasserstein distance methods. Article [83] has used moment-based DRO focusing on cyber-resilience alongside economical objectives taking into account the uncertainty of wind power generation. The other moment-based DROs are [102], [119], that consider the uncertainty of renewable energy generation. The work in [100] used the Wasserstein distance for the DRO formulation, the considered uncertainty again being renewable generation. In addition, [3] has used a CVar-Based DRO to deal with the uncertainties of loads and renewables. Article [108] has applied the DRO approach to a water-energy nexus management problem, which indicates an operational problem of a networked multi-EH system under renewable energy uncertainty.

4) ADAPTIVE DISTRIBUTIONALLY ROBUST OPTIMIZATION
The combination of the second and third categories above, results in the ADRO framework. Only one paper, [96], has investigated a two-stage data-driven optimization framework for an EH. The paper deals with the uncertainty of PV production using a moment-based ADRO.

D. CHANCE CONSTRAINED
Chance constrained programming was originally developed by Charnes and Cooper [132]. The formulation enforces the so-called ''chance constraints'', which limit the probability of violating a certain constraint to be below a certain threshold value. The indicated CC approach may be used based on the following general formulas for function minimization (14), deterministic constrains (15), and separate chance constraints (16), [132], [133]. In equations (14) -(15), x is a decision vector, and e is an uncertain multi-dimensional parameter vector. The CC formulation is a constrained optimization problem (as shown in 16), where this time e is s a random vector defined on some probability space and η is a tolerance probability. The CC method allows for a probability to violate the constraints in the presence of uncertainties [134]. In contrast to the other mentioned methods, CC is only applied to the constraints. In other words, a CC application provides a feasible solution to the original problem with a probability at least η.
Among the existing works, [135] utilized the CC method to cope with the uncertainty of the driving pattern of EVs. It applied the CC to the constraints of the state of charge of the EVs. Article [136] implemented the CC on power flow equations due to the presence of the fluctuation of wind power generation. Similarly, article [137] has deployed CC to transmission constraints in the presence of renewable energy uncertainty. Papers [138], [139] have developed CC programs to consider the constraints of power flows and pipe-line flows. All mentioned works have focused on the operational cost of EHs. However, [103] used CC programming to deal with the uncertainties of wind and PV generations as well as electrical loads in a planning problem. It addressed the constraints of the power balance and capacity limit of EH devices. Table 4 summarizes the studies that applied CC programs for EHs.

IV. HYBRID METHODS
Various papers in recent years have used hybrid methods to deal with uncertainties. In the majority of them, each uncertainty source is investigated by one of the mentioned methods. Table 5 summarizes the use of hybrid methods in the relevant literature.
The relevant hybrid methods can be divided into five categories:

A. COMBINATION OF SP AND RO
The authors in [35] have taken into account the uncertainties of the EV driving pattern using SP, while the uncertainty of the electricity market prices was dealt using RO. The uncertainty of wind power is managed through RO in [42], and the uncertainties regarding electricity, heat, and gas demands through SP. The uncertainties of wind speed and electricity prices are handled by SP and RO respectively, in [44], where the main objective is to minimize the total operation cost of the EH. Reference [56] considers EVs' uncertainty parameters using RO, whereas the uncertainties of demands, market prices, reserve requirements, and renewable power are handled via SP. Among other works, [36], [39], [45], [46] have utilized SP to deal with the uncertainty of wind power (or PV generation), while deploying RO to cope with the worst-case realization of electricity market prices. Uncertainty of electricity demands is again handled using SP in [36], [45], and [46] and EV driving patterns (arrival and departure times of EVs to/from parking) are handled by SP in [36], [45], and [46]. A timely observation is that uncertain energy demands tend to be modeled using some known probability distribution (allowing for SP methods to apply), whereas wholesale market prices can be more volatile and/or unpredictable, consequently handled using RO formulations.

B. COMBINATION OF SP AND IGDT
In [38], a hybrid method combined IGDT and SP to deal with the uncertainties of wind generation, energy demands, and electricity market prices. The SP has been developed to take into account the uncertainties of the demands and electricity prices, and IGDT has been used to cope with the uncertainty of wind power. This is the only hybrid IGDT-SP method that addresses the uncertainty of the wind speed by IGDT. Other hybrid IGDT-SP papers, including [41], [51], [53], [65], have deployed SP to take into account the uncertainty of wind generation. In addition, [51] uses the IGDT to reach the worst-case realization of electricity prices and demands. Study [41] has developed the IGDT to deal with the uncertainty of the demands of EVs, and SP to face the uncertainty of PV generation. Reference [53] forms a non-linear objective function by having the electricity price uncertainty handled via IGDT, which creates a bi-linear term. The other uncertainty of this article is the energy demands which were handled using SP.

C. COMBINATION OF RO AND IGDT
Only one paper - [61] -has focused on the combination of RO and IGDT in the operation and planning of EH so far. The authors deployed RO to deal with the uncertainty of electricity market prices and IGDT to take into account the uncertainty of wind power. The main purpose of the combination of the two mentioned methods is to propose a linear and tractable objective function. It is worth mentioning that, including market prices in the IGDT method results in a nonlinear framework, as discussed in the previous paragraph.

D. COMBINATION OF CC WITH OTHER METHODS
CC is quite different than the other mentioned methods, in the sense that it only deals with the constraints instead of the objective functions. A few papers have considered the combination of CC alongside other methods. Article [48] used a hybrid method combined by SP and CC for configuration problems of the EH involving the uncertainties of energy demands where CC capture the probability of load curtailment. Paper [49] proposed an operational model for an EH in the presence of uncertainties of PV and energy demands (i.e., electrical, thermal, and cooling). The proposed framework used SP, while the constraints regarding demands; satisfaction were modeled using CC. Similarly, heating, cooling, and electricity have been considered as outputs in [50] where, in addition to the operation of the EH, the authors proposed a model for thermal services. The uncertainty of renewable energy, electricity demands, and ambient temperature has been considered using SP, while the CC part handles the thermal service quality. Paper [47] also represented the demand satisfaction by a CC due to the existence of the fluctuation of PV generation and electricity demands in a multi-hub operation problem.
Finally, article [92] has investigated the energy management problem of the retailers under the uncertainties of renewable energy and electricity demands. The main balancing constraint has been modeled as a chance constraint, while the proposed framework is robust against the mentioned uncertainties.

E. OTHER HYBRID METHODS
Under this category, we refer to methods that generaly include at least one of the mentioned methods and do not fall into the previous categories. A two-stage stochastic method has been proposed in [37], where the first stage formulates the uncertainties connected to the demands and renewable power using the two-point estimate method, while the second stage VOLUME 11, 2023   deploys the IGDT to take into account the uncertainty of gas prices.

V. SUMMARY OF RESEARCH TREND
In this Section, we present a high-level bibliographic analysis of the area of mathematical optimization applied to EH decision problems. The number of published papers in the area of EHs is shown in Fig. 4. As it is obvious from the figure, the total number of papers is almost continuously increased. With the exception of the limited number of articles published between 2011 and 2014, it is clear that in the early years (2015 to 2016) the greatest focus was on applying SP methods to EH problems. The main disadvantage of these methods is that they need a lot of information and their computational burden is high. This problem forced researchers to drastically reduce the number of scenarios, which could lead to inappropriate management of uncertainties. For these reasons, the use of robust methods, such as IGDT and RO, which require less information and has a lower computational burden than statistical methods, has been steadily increasing since 2017. However, these methods are very conservative, increasing the average cost of operation and planning of EHs. In contrast, hybrid methods that borrow the advantages of both classes naturally come as the next step in the area. This observation suggests a trend that expects the application of hybrid methods to be increasing in the years to come. This trend can indeed be observed by the steady increase of papers applying hybrid methods in recent years.

VI. CONCLUSION
In this paper, the application of mathematical optimization approaches to the operation and planning of EHs has been presented and reviewed. The main part of the paper has been dedicated to investigate and categorize existing mathematical methods for dealing with the uncertainties pertaining EHs, including stochastic programming, information gap decision theory, robust optimization, and chance-constrained methods. The literature has also been categorized on the basis of the models' inputs, controllable EH resources considered, and corresponding outputs of the EH management module. Our conclusions can be summarized as follows • SP is an effective approach that has been developed to consider the uncertainties (adhering to known probability distributions) in an EH or multi-carrier energy system, e.g. demands and electricity prices, wind speed, solar irradiation.
• Studies that have used IGDT to optimally make decisions under uncertainty have been categorized into three classes depending on the number of uncertain parameters. Even though part of the literature uses IGDT towards modeling EH management problems, the method is not deemed particularly suitable due to the fact that it cannot properly manage multiple uncertain parameters.
• The investigation of the works that used robust optimization methods demonstrates that the most common approach points to the application of classic and adaptive robust optimization, mostly in order to cope with the uncertainties of electricity market prices and sometimes also renewable energy sources.
• The chance-constrained method has been used in a few papers, in most of which it is used to constraint the probability of load shedding.
• Hybrid methods have been widely deployed recently to harvest the advantage of various methods in dealing with uncertainties. Most characteristically, multiple studies model the uncertainty of demands using a known probability distribution (pointing to SP), while the uncertainty of market prices (which are more volatile) is handled via robust optimization.
• In terms of input, output, and controllable resources, the penetration of hydrogen energy and as result, the related devices such as electrolyzer within EHs, as well as water management as an additional system that interacts with the energy systems, can be seen in several recent works. Further recent considerations include the application of P2X devices, e.g., power to gas (P2G), which allow EHs to use the surplus renewable energy to generate natural gas. Considering the extracted research trend in this area, the following directions for future work will be highly interesting: • to study how the newly introduced technologies (namely, hydrogen and P2X) penetrate into the envisioned EHs in the years to come.
• to study the operation of waste heat that comes from the conversion of energy carriers in the operation of EH.
• to develop EH use cases and real business models due to the expected substantial growth in this area.
• to develop some real-time operational models to enable the operation of the real-world EHs.
• to investigate the interaction of energy markets in EHS, i.e., electricity, natural gas, heat, and carbon.
• to study the application of machine learning and reinforcement learning to the operation and planning of smart EHs.

APPENDIX B APPLICATION OF MACHINE LEARNING TO OPERATION AND PLANNING OF ENERGY HUB
Machine learning methods have been applied to both optimal operation [175], [176], [177] and optimal planning [178], [179] of EHs. They have a huge potential to improve the EHs' energy management, particularly when real-time control is needed [180], [181], [182], [183]. A distributed energy management approach based on multi-agent reinforcement learning has been applied to residential EHs to minimize the operation cost of EHs [184], [185]. In [180] a Markov decision process has been developed for real-time management of EHs by minimizing carbon emission and energy cost. In [186] deep reinforcement learning has been used to develop a datadriven model-free framework for the energy management of EHs. Reference [187] has proposed linear programming of the reinforcement learning approach to overcome the environmental complexity of the problem. Multi-agent deep deterministic policy gradient has been used in [188] to develop a reinforcement learning-based approach for optimal operation of multi-EH. In [189] a deep Q-learning has been used to maximize the long-term profit of the EH's prosumers in both the local energy market and wholesale energy market. A hybrid data-driven and model-driven framework have also been developed by using a machine learning algorithm in [190].