Scheduling of Energy Hub Resources Using Robust Chance-Constrained Optimization

This paper develops a robust chance-constrained model for handling the uncertainties of generation and consumption in multi-carrier energy hubs. The proposed model incorporates corresponding loading factors for each type of electrical, heating, and cooling loads. This is done to assess the maximum loadability of the whole system. In this respect, the chance-constrained approach is implemented for the feasibility assessment of the operation problem with uncertainties. The uncertainties which are assumed here include the forecast errors of electrical, heating, and cooling load demands, and the volatile solar power generation. The overall problem formulation is developed in the mixed-integer linear programming (MILP) framework. The standard chance-constrained approach is converted to a deterministic optimization model by utilizing the Big M method. The main objective of the proposed model is to maximize the loadability index with uncertainties while addressing the permissible risk index of the decision-maker. The studied energy hub comprises electrical, heating, and cooling loads, and the energy flow technique is adopted in this paper to model the load balance equations. The simulation results are presented for different scenarios while addressing features of the proposed model for the summer and winter seasons. Furthermore, the developed model is evaluated for different scenarios and a comparison is made with the information-gap decision theory (IGDT) method.


Indices/Sets
Energy stored in the EES system (kWh). P PV s,t Solar power generation (kW).

I. INTRODUCTION
Energy hubs (EHs) are recently emerging technologies including multiple energy carriers where they can be interconverted to other types to satisfy the electrical, cooling, and heating load demands in an economic manner while providing more flexibility for the energy system as a whole. One key challenging issue is the usual day-ahead resource scheduling problem, which is the target of this paper. Our aim is to model and solve the day-ahead operation problem of one EH subject to intrinsic uncertainties in the electrical, heating, and cooling load demands. Consequently, the resulting problem is a three-dimensional stochastic optimization problem whose solution aims to decrease the potentially excessively high operating costs of this scenario. As to be discussed throughout this paper, an effective optimization tool is used which is called ''robust chance-constrained programming'' to tackle the problem in a computationally efficient manner.

A. MOTIVATION
One of the most challenging issues in the short-term operation of EHs is handling the uncertainties from volatile renewable energy sources and energy consumption by the end-users. This paper presents a chance-constrained optimization model augmented by a loadability index to increase the robustness of serving the electrical, heating, and cooling loads of the EHs. The main advantage of this approach is that the operating points of the hub assets can be optimally determined while the risk of operation due to the load shedding can be minimized. To deal with the energy conversion and energy transition between the hub assets, the energy flow method is employed, because it is more effective than the matrix-based representation of the multi-carrier energy system, specifically when there are dynamic energy storage devices. The mathematical formulation of hub operation and the proposed robust chanceconstrained model is worked out as a mixed-integer linear programming (MILP) approach, and thus, the computational burden of the mentioned problem can be reduced.

B. LITERATURE REVIEW
The optimal operation of multi-carrier energy systems is one of the prominent research topics in recent years, with an extensive focus on minimization of the total operating cost of the EHs while serving the electrical, heating, and, cooling loads. Some papers have been published on the modeling of EHs with deterministic load demand [1], [2], [3]. However, there are still challenging issues regarding the mathematical modeling and problem formulation of hub assets and their functionality, such as the computational efficiency of the model and the capability to address uncertainties. The effects of uncertainties in the generation and consumption on the optimal dispatch of the hub assets have remained a controversial topic in the specialized literature. VOLUME 10, 2022 Specifically, there has been relatively little work done about uncertainty handling in the short-term operation of multicarrier energy systems. Several studies, carried out thus far on the operation of EHs while addressing the stochastic behavior of the problem can be found in [4], [5], and [6]. A scenario-based optimization framework considering the stochastic behaviors of natural gas (NG) and electricity tariffs in line with the uncertainties of electrical load forecasting has been presented in [4]. The minimization of the overall operating cost and environmental emissions in a weighted sum framework with the concept of conditional value at risk (CVaR) method was suggested. A two-stage stochastic optimization model was developed in [5] to determine the optimal reserve capacity and energy scheduling for an EH. The expected cost of energy serving has been introduced as the main objective function. The uncertain parameters in the mentioned study were solar and wind power generation forecasts as well as the electrical and thermal loads. To promote energy efficiency, the interdependency among the natural gas and electrical energy systems was studied, and the reliability and security assessment for the given EH has been carried out. A bi-level stochastic optimization approach has been presented in [6] to show the effectiveness of multiple EHs with a simplified functionality of the control strategy for energy transactions. To reduce the unwanted impacts of uncertain parameters on the optimal operating points of hub assets, the CVaR approach has been adopted in the mentioned research. A stochastic optimization model taking into account multi-energy systems operation and power exchange with the electricity market has been investigated in [7]. The uncertain day-ahead and real-time clearing prices have been treated as scenarios for the proposed model. The expected operating cost and the potential risk using CVaR approach have been considered in the aforementioned model to optimize the dayahead EH scheduling within the energy market. The mathematical optimization problem has been established as a robust optimization to address the secure operation in the worstcase scenario. The main focus of the research is to handle the conservativeness and computational load of the problem in a stochastic-robust coordinated manner. The optimal scheduling of an EH in a stochastic programming paradigm has been studied in [8] addressing the heating market impacts on the EH operation and heat demand response in line with the electricity demand response. The uncertainties of price and wind power generation have been studied in the problem formulation. Recent developments in the field of robust optimization techniques have led to an increasing interest in renewing the mathematical representation of EH operation with considerable uncertainties in both input and output parameters, i.e. generation and consumption sides. The wellknown information gap decision theory (IGDT) approach has extensively been applied to the problem of optimal operation scheduling of EHs [9], [10], [11], [12], [13]. A new mathematical formulation has been developed in [13] to address the uncertainties of the EH operation problem while the combined heat and power (CHP) unit has been modeled with a convex feasible operating region (FOR); however, the cost function of the CHP has been represented as a non-linear function. A risk-averse technique has been introduced in [11] dealing with the optimal operation of a multi-carrier energy system considering the impacts of plug-in electric vehicles in the model. An IGDT-based mathematical representation of the robust optimization problem has been investigated in [10] taking into account the uncertain load demand of plug-in electric vehicles. An IGDT-based model for investigating the impacts of the interdependency of heat and electrical power production in CHP units has been addressed in [9]. The electrical and thermal energy storage devices have been considered in the model and the uncertainties due to the forecast errors of the electricity price, renewable power generations, as well as load demand have been characterized in the simulations. Besides, the load flow constraints have been accommodated in the mathematical problem formulation. A new robust optimization model based on the extended affine arithmetic approach has been introduced in [14] to tackle the optimal operation problem of EHs in the presence of heterogeneous multiple uncertainty sources. A robust optimization approach was developed in [15] using the control approach to cope with the bounded uncertainties on EH parameters. The obtained solutions reported were feasible for all values, for a given subset, of uncertain parameters in the given simulation case study. Another robust optimization technique has been presented in [16] for robust scheduling of multi-carrier EHs with techno-economic/environmental limitations, affected by the market price uncertainty and demand response mechanisms. The proposed model was evaluated in different timebased demand response programs. Price fluctuations were addressed in the mentioned study and the suggested model guaranteed the minimum global optimal operating cost. A multi-objective robust optimization framework has been proposed in [17] handling the risks of combined demand and supply uncertainties for smart residential end-users. A tradeoff between the model robustness and the solution robustness has been made in this study with a substantial cost saving achieved by applying the robust optimization approach. Another robust optimization method has been developed in [18] to deal with the robust scheduling of residential buildings in the presence of price uncertainty. A precise margin has been adopted for market prices instead of electricity price forecasting to characterize the uncertainties of hourly electricity prices. A comprehensive simulation has been carried out to address the functionality of the proposed model in a case study with 10 smart buildings. A novel hybrid robuststochastic optimization technique has been investigated in [19] for the bidding strategy of large-scale prosumers. The robust and stochastic optimization approaches have been incorporated respectively to model the uncertainties arising from load and market prices. Jamalzade, et al. [20], presented an optimal operational strategy for EHs by employing a hybrid stochastic-interval optimization method to capture the demand uncertainties. In order to apply deterministic programming to the EH management problem, the authors of Ref. [21] used the Cornish-Fisher expansion approach to convert probabilistic constraints into deterministic ones. The created model was solved using the interior point approach due to the robust performance of the presented model. A twostage risk-oriented stochastic p-robust optimization-based scheduling strategy was developed in Ref. [22] for EHs. In the first phase, the EH may trade energy on the day-ahead electricity and thermal markets. In the second phase, the EH controls the surplus/deficit of renewable energy output on the real-time electricity and thermal markets. Reference [23] proposed a distributionally-robust optimization framework for the day-ahead scheduling of EHs, aimed at maximizing social welfare. A two-stage optimization model has been suggested incorporating interval optimization, and the wellknown weighted-sum and fuzzy satisfying methods have been used for solving the mathematical optimization problem. Co-optimization of the energy and reserve markets has been carried out in [24] and the operation problem of the multi-carrier energy systems has been dealt with as a twostage robust optimization model. Demand response programs have been applied to increase the flexibility of the operation problem. The reported results confirmed that the proposed model can effectively increase clean energy production as a holistic goal of demand response incorporation. So far, chance-constrained programming has also been used for the operation of EHs. In this regard, to improve the flexibility of the multi-EH system including CHP and CCHP systems, and achieve a reliable operation for the system, a novel optimization framework was provided in Ref. [25] using chanceconstrained programming and multi-objective optimization. A chance-constrained optimization approach was used in Ref. [26] to tackle an MILP optimization problem. The optimization was developed, aimed at supplying the electricity and cooling loads of a data center and the hydrogen demand of a neighboring hydrogen fuel station. In this regard, renewable power curtailment, along with other assets is used in the EH. Using a chance-constrained optimization method, the uncertainties associated with the use of renewable energy curtailment have been accounted for. Ref. [27] proposed a chance-constrained optimization technique for the uncertain operational planning of EHs. By employing a convexification technique, the nonlinear formulations of energy and gas flows were handled and relaxed. A comparative overview of the recently published papers is given in Table 1. Besides, a comprehensive review of the optimization methods used for the EH management problem has been provided in Ref. [28].

C. CONTRIBUTIONS
The main contributions of this paper with respect to the literature are as follows: Investigating the EHs by using the energy flow model; the energy flow model can be used to address the dynamic features of the energy storage devices and the input-output functionality of each asset can be extensively addressed. • Proposing a MILP model to tackle the EH operation problem; the optimal scheduling of EH's assets in the presence of diverse producers, converters, storage systems, and consumers is a complicated optimization problem. Therefore, the scenario-based stochastic optimization problem is formulated as the standard MILP model in this paper to reduce the computational complexity.
• Handling the uncertainties using the robust chanceconstrained approach; The stochastic optimization model is tackled as a robust chance-constrained optimization problem addressing the forecasting errors of both generation and consumption profiles. Furthermore, the robustness of the solutions is evaluated by using the loadability index.

D. PAPER ORGANIZATION
The remainder of the paper is organized as follows where the fundamentals of the chance-constrained programming are described in Section II and, Section III comprises the VOLUME 10, 2022 mathematical modeling of the EH day-ahead operation problem. Simulation results are given and discussed in Section IV and lastly, concluding remarks are included in Section V.

II. CHANCE-CONSTRAINED OPTIMIZATION APPROACH
One of the biggest current challenges relates to the tackling of large-scale problems. The chance-constrained method is one of the most efficient approaches used in stochastic optimization problems with a high uncertainty level. It constrains a number of more unlikely scenarios so that the decisionmaker can choose the level of reliability and risk intended and considered adequate. This type of solution is robust, but in real-dimension problems, it may become difficult to solve. The typical formulation for this type of problem is as follows: It is noteworthy that the decision vector, the uncertainty vector, as well as the equality and inequality constraints, are denoted by x, ξ , g, and h, respectively. In this regard, the equality constraints would be rewritten as follows by using the chance-constrained method [44]: where ε indicates the risk level to be specified by the decisionmaker. The chance-constrained programming is solved in several ways. One of the solutions may involve reformulation and solving the problem by the equivalent deterministic problem. This problem can be derived through the Big M method or through the bilinear reformulation. By employing the Big M method, a binary variable is used indicating whether the associated scenario should be considered or may be violated, and thus the problem is converted into MILP. Equation 2 can be converted by utilizing the Big M method as follows: where z ω is the binary variable representing whether or not the scenario is active; M represents the Big M parameter which should be sufficiently large. Moreover, ε represents the predefined risk level, defined by the decision-maker.

III. PROBLEM FORMULATION A. ENERGY FLOW MODEL
The mathematical modeling of the optimal operation problem of the EH is represented in this section and Fig. 1 depicts a typical EH, including a microturbine (MT), a battery energy storage (BES), a TES, a CHP unit, a solar PV system, a boiler, and a wind turbine (WT), serving electrical loads (ELs) and Thermal loads (TLs). This system has bidirectional electrical power transactions with the electrical grid and the required natural gas is supplied through the natural gas network. • Objective Function The objective function represented in (5) comprises three items, where the first one indicates the cost due to transacting electrical power with the upstream network. This item is derived by multiplying the amount of transacted power at each hour by the corresponding market price. The second part of the cost function relates to the NG-powered generation units, i.e. the boiler and the CHP unit. The third item indicates the cost due to load shedding, either electrical, heating, or cooling loads. In this respect, the penalty factor is denoted by and determined such that the load shedding cost is minimized. Eq. (6) and Eq. (7) show the generation cost functions of the CHP unit and the boiler while taking into account the electrical and heating efficiencies, respectively, as in (5)-(7), shown at the bottom of the next page.
• CHP The FOR of the studied CHP unit is convex and the constraints of the electrical and heating power generation have been shown in (8)- (10). Constraints (11) and (12) • Boiler The constraint of heat generation by the boiler is expressed in (13) while the heat flow equation has been shown in (14). This relationship states that the heat output of the boiler can be delivered to the heating load demand,, or even the AC, to provide the cooling power.  (14) • EH The EH's feasible operating interval has been indicated in (15). The heat power equation of the EH is shown in Eq. (16) in which the electrical power-to-heat conversion equation has been stated. The efficiency of EHs is relatively high and denoted by η EH . The electrical power required by the EH would be supplied through the solar photovoltaic (PV) panel, P PV →EH . This limitation has been modeled in the relationship (17).
• EHP EHPs are devices used to provide heating or cooling power by using electrical power. These devices are not able to concurrently work in the mentioned modes. Thus, a binary variable would be defined to specify each mode. Constraint (20) is used to apply this operational limitation. The heating power generation relationship and cooling power generation relationship are expressed in constraints (21) and (22), respectively. Some research studies have used the concept of coefficient of operation instead of efficiency [45]. The input power balance equation is indicated in (23) and the output power balance equations are stated in relationships (24) and (25) H EHP C EHP • AC Relationships (26)-(29) illustrate the mathematical model of the AC. The cooling power generation constraint is stated in (26) while the heat-to-cooling power conversion equation is shown in (27). As previously mentioned, the heating load demand of the AC can be met by the CHP unit, H CHP→AC s,t and the boiler, H Boiler→AC s,t , indicated in (28). The cooling power output of the AC is directly delivered to the cooling load demand, C AC→CL s,t , as indicated in (29).
• EES The EES system is regarded as a key asset of the EH and it is modeled as (30)- (37). As constraint (30) indicates, the amount of energy available in the battery at every time slot of the scheduling period should fall within the feasible operating interval specified by the manufacturer and the battery operator. Furthermore, as stated in the relationship (31), the amount of energy available in the battery at each time slot is defined as the function of the energy available in the system in the previous slot plus the charging power and minus the discharging power, taking into account the efficiencies of these two operating modes. The hourly charging and discharging power constraints of the battery are modeled through relationships VOLUME 10, 2022 (32) and (33), respectively. It is worth mentioning that the battery can work in one of the charging, discharging, or idle modes at a time as emphasized in constraint (34). To meet the operational requirements for the subsequent scheduling period, the amount of energy available in the battery once the scheduling period is over, should meet its initial value.
• PV The hourly solar power generation by the PV panel is a function of the solar irradiance, ambient temperature, and also manufacturing characteristics of the panel. The electrical power generation equation of the PV panel is stated in (38) while constraint (39) applies the power generation limitation [46]. Besides, as expression (40) (40).
• Grid The power exchange between the hub and the main grid is characterized by using relationships (41)- (45 , as stated in (42). The amount of power transaction at each time slot is limited to the transformer, connecting the hub to the grid as shown in (43) and (44). It is worth noting that the variable showing the hourly power transaction with the grid is a positive variable. It should be noted that it would not be permitted to concurrently import/export power from/to the grid as emphasized in (45).
• Load Balance The most significant constraints of the problem of any multicarrier energy system operation problem are electrical, heating power, and cooling power balance equations stated in relationships (46) . In case the electrical power generation by the assets and power imported from the main grid do not meet the load demand, load shedding would occur. It is noteworthy that all variables are positive variables. Likewise, the heating and cooling power balance equations are expressed in (47)

B. ROBUST CHANCE-CONSTRAINED APPROACH FOR SOLVING EH OPERATION
The robust chance-constrained optimization has been deployed in this paper to ensure the robustness of the solution against the uncertainties due to solar power generation and the three load types. In this regard, the objective of the system operator is to maximize the loadability of the EH with the load demand supply, associated with the probability 1-ε, where ε is the risk index of the decision maker. The higher values of ε would ensure the higher loadability of the EH against the uncertainties. Accordingly, the balance equations of electrical, heating, and cooling power would be rewritten as (49)- (51), shown at the bottom of the next page, respectively. The Big M technique is used to transform the stochastic problem into a deterministic one. As a result, constraints (49)-(51) would be rewritten as constraints (52)-(54), and the corresponding constraints would be (55)-(57) and (58) where z is an auxiliary binary variable to ensure the balance between generation and consumption. It is noteworthy that the total loadability of the system is studied showing the worst scenario for the load demand increase.
Constraints (55)-(57) indicate that the load shedding is a positive variable where if it is true, its associated binary variable would also be true. Accordingly, the primary optimization problem is converted into the optimization of the loadability of the EH. Thus, the problem will be iteratively solved and intended to maximize α. Then, the maximum loadability would be specified, subject to minimizing the operating cost of the hub. The conceptual structure of the proposed robust chance-constrained optimization framework is demonstrated in Fig. 2.

C. ROBUST IGDT APPROACH FOR SOLVING EH OPERATION
This part describes the implementation of the IGDT approach into the deterministic framework that was provided in the previous section to accommodate the extreme uncertainty that was caused by the power demand and the intermittent nature of wind generation. The overall operating cost may be expressed as a function of the uncertain source, where X indicates the vector of choice factors [47], [48].
It is noteworthy that the IGDT model deployed in this paper uses the envelope bound model as [49] and [50], expressed in (62) and (63)   s,t fromH HL s,t less than αH HL s,t . A significant advantage of the IGDT is that it makes the decision maker, the system operator in this case, to prevent the risk of achieving the lowest expected values taking the parameter uncertainty into account. The robustness function is an effective risk assessment tool. In the IGDT approach, the robustness function is defined by the highest values that α may achieve at a cost less than the maximum expected cost shown by R c [51], [52].

RF(k, P EL s,t , H HL s,t )
= Max α (α) : Maximum cost which is not higher than a given biggest cost where RF(k, P EL s,t , H HL s,t ) indicates the input/output architecture of the system model. That is to say, it shows the operator's award for the selected values of decision variable k taking into account the uncertain parameters P EL s,t and H HL s,t . At the highest degree of uncertainty, the robustness function provides the best performance, meaning that the operating cost is less than the predetermined cost R c . Thus, the robustness function represents the performance of risk-hedging. The greater the value of this robustness function, the more solid, risk-hedging, and impervious to existing uncertainties the decision. A risk-hedging operator must adhere to a schedule that limits exposure to losses or excessive cost levels. Consequently, the robust performance may be described as in (65), shown at the bottom of the next page. R c expresses the critical cost, while R 0 represents the minimal expected cost based on the predicted input factors. σ is the cost aberration factor used to determine the greatest expected cost. The purpose of the IGDT in the robust EH scheduling for risk-hedging strategy is to maximize the uncertainty parameter α so that the desired performance is achieved.
According to the IGDT technique, the proposed optimization framework would be expressed as follows, (66)-(70), as shown at the bottom of the next page.

IV. SIMULATION RESULTS
The simulation results and the case study including an EH equipped with a solar PV panel as well as electrical, heating, and cooling loads, are discussed in this section. Table 2 represents the equipment of the EH together with the associated data [45]. The analysis has been carried out for two seasons, i.e. winter and summer. The uncertain parameters of the problem are electrical, heating, and cooling load demands besides the volatile solar power generation. It is noteworthy that the uncertainties have been characterized by generating scenarios and an efficient scenario reduction approach has been used to alleviate the number of scenarios, i.e. 10 scenarios for each uncertain parameter. Fig. 3 depicts the market price where the selling price and market price are the same in summer because of the relatively high load demand of the main grid. The selling price refers to the price of energy sold to the main grid by the EH. In the winter, the amount of load demand of the main grid is relatively lower and accordingly, the selling price would be 80% of the market price. The scenarios used in this study for the uncertain parameters in the summer and winter have been depicted in Fig. 4.

A. CASE A: STOCHASTIC OPTIMIZATION RESULTS
The problem of optimal operation of the EH is tackled in this case as a stochastic optimization problem, intended to optimize the total operating cost. In other words, the problem is investigated while skipping the risk measure, i.e. ε =0. In this relation, the problem is studied for the base case without considering the loadability index. The simulation results indicate that the expected value of the total operating costs for a typical day in the summer and winter would be $1582.766 and $1530.01, respectively. In the summer, the NG price is considerably low at 0.006 $/kWh and as a result, the CHP is scheduled to generate more electricity. Hence, the CHP unit is employed at its maximum capacity to generate electricity, and heat generation would be set at the permitted value. A fraction of the heat output of the CHP unit would be used as the input heat of the AC, while the remaining would be deployed to satisfy the heating load demand over the day. The boiler would also be employed to serve the heating load demand and for heating-to-cooling power conversion in the AC. The hourly heat output of the boiler is illustrated in Fig. 5. The AC and EHP are deployed to supply the cooling power demand in the winter. It is worth mentioning that the cost of the NG used by the AC is more than the cost of electricity consumed by the EHP to supply the heating load demand. Thus, the operator tends to employ the AC to supply the cooling load demand rather than generate heat. However, the AC alone would not be sufficient to thoroughly meet the cooling load demand and accordingly, the EHP would also be necessarily utilized to this end. The amounts of cooling power delivered by the AC and EHP in the summer are depicted in Figs. 6 and 7, respectively. The amount of electrical energy available in the battery over the day in the summer is demonstrated in Fig. 8. This device is charged during the initial time slots of the scheduling period and delivers power to the system from time slot 11, associated with high energy prices. The battery is again charged during the final time slots to meet the pre-scheduled energy value of 200 kWh at the end of the day. The amount of power transacted between the EH and the utility grid in the summer is illustrated in Fig. 9. The hub exports power to the utility grid over the initial time slots of the scheduling period as its electrical load demand is significantly low. It is noteworthy that the price of electrical energy provided by the EH taking into account the heating power generation would be approximately 12 $/MWh. Meanwhile, the minimum market price is 16 $/MWh. Consequently, it is economical to deploy the CHP unit over the day for electricity generation. Therefore, the EH is capable of selling the surplus energy to the main grid and benefits from the power transaction. The EH imports power from the main grid during hours 8-21, and sells electrical energy to the grid during hours 1-7 and hours 22-24. The hourly amount of NG purchased from the NG network for different scenarios in the summer is depicted in Fig. 10. The results obtained from simulating the problem for the winter are different due to the different load profile in this season. The amount of heating load demand is substantially higher in the winter, while the NG price is 0.0085 $/kWh. Furthermore, the electricity selling price in the winter is lower than in summer, i.e. 80% of the market price. This issue implies the fact that any transaction with the main grid to sell the surplus electrical power would be possible in case the electrical energy generation is economical with respect to the market price.
The efficiency of the CHP unit is around 40% and according to the NG price and market price, selling electrical energy generated by the CHP unit to the grid would be possible if the market price is at least 26.56 $/MWh, while the price over on-peak hours is 24 $/MWh. So it is not economically justified to sell electricity produced by the CHP system to the main grid. However, it should be noted that if all assets operate together, the electricity generation cost for the EH would be around 21.25 $/MWh, which is lower than 24 $/MWh. Thus, the hub can sell the surplus power to the main grid at time slots with a market price higher than this value. The hourly electrical power and heat outputs of the CHP unit for different scenarios in the winter are demonstrated in Figs. 11 and 12, respectively. As expected, the CHP unit generates power during the time slots at which the market price is 24 $/MWh, i.e. time slots 8-22, and it is not economical to operate at other slots. On the other hand, the capacity of the transformer linking the hub to the main grid is 300 kW and the amount of peak electrical load is lower than this amount.   So, it is not needed to use the CHP unit during the initial and final time slots of the day. The boiler is also capable of satisfying the peak load demand. The heat generation cost of the boiler is 14.16 $/MWh (8.5/0.6=14.16 $/MWh), which is lower than the electricity price over the entire day. As a result, it would not be economically justified to use the EH and EHP to supply the heating load. The efficiencies of the EH and EHP are 85% and 90%, and accordingly, the heat generation cost of these two assets would be 17.64 $/MWh (15/0.85=17.64 $/MWh), and 16.67 $/MWh       it is yet not economical to use these devices. In this regard, the heating load demand of the hub would be supplied by the boiler and CHP unit. The hourly heat generation of the boiler is shown in Fig. 13. This asset operates at its maximum capacity over the initial and final hours of the day with considerable heating load demand. The boiler along with the CHP unit provides the heating load demand of the AC over the day. It is noteworthy that the AC can supply the entire cooling load demand, and there would be no need to use the EHP. The hourly energy available in the battery can be observed in Fig. 14. As this figure depicts, the battery is charged over the initial hours of the day with low electricity prices, and it delivers power to the system during the hours with relatively high market prices. Finally, the battery is charged over the final hours of the day to meet the constraint of the final available energy, i.e. 200 kWh. The simulation results in the base case for the studied scenarios show that the system operator would be able to reliably supply the electrical, heating, and cooling load demand, and the operating points of the assets are in the permitted operating ranges. The next section evaluates the problem in the probabilistic and robust state.

B. PROBABILISTIC OPTIMIZATION RESULTS
This case study evaluates the loadability of the EH considering the problem's uncertainties. In this regard, the worst scenario faced by the system operator is intended, i.e. the concurrent increase in the three types of load demands. To this end, the load demand is continuously increased until the solution becomes infeasible. This is done for different values of ε specified by the system operator. The obtained  simulation results show that the loadability of the EH is relatively higher in the winter compared to that of the summer. This is due to the fact that a substantial fraction of the load relates to the heating load demand which can be supplied by utilizing the boiler, EH, EHP, and CHP unit. It is noteworthy that first, the heating power-related constraints cause the solution infeasibility. The electrical and cooling load demands are considerable in the summer. The cooling load demand can be supplied by using the AC and EHP, while the AC operates at its maximum capacity during onpeak hours. On the other hand, the EHP consumes electricity which in turn causes the electrical load demand to increase. Thus, an increase in the cooling load demand would directly impact the electrical load demand. The electrical load demand supply would encounter severe difficulty since the capacity of the transformer is limited to 300 kW. Table 3 includes the simulation results, obtained from studying the EH loadability in different seasons. The results verify that the loadability of the EH increases by increasing ε, i.e. a more risk-taking operational strategy. In this regard, 5% and 10% increases in ε result in 22% and 32% increases in the loadability of the hub compared to the case without any load shedding. For the 5% and 10% increases in ε during the summer, 24% and 31.5% increases in the loadability of the hub compared to the base case have been observed. Furthermore, Fig. 15 depicts the simulation results for different values of α while ε =0 in the summer and winter. The obtained results show that individually increasing any of the three load types would result in different loadability indexes in different seasons. As it is expected, the highest value of the loadability index in the winter pertains to the cooling load demand by α = 4.785. The loadability indexes for the heating and electrical loads are 1.495 and 1.486, respectively. The loadability for the heating load demand in the summer is 7.573 while the loadability indexes for the electrical and cooling loads are 0.378 and 1.270, respectively. These results mean that the total loadability index which is equal to 0.308 in the summer is highly dependent upon the loadability of the electrical loads.

C. IGDT RESULTS
The results obtained from the chance-constrained programming are compared to those obtained from the IGDT technique to validate its performance. In this respect, the robust optimization problem is simulated for two days in winter and summer. The main difference between the IGDT and chanceconstrained programming models relates to the power balance constraint that must be satisfied in all scenarios without any curtailment. The system loadability can be assessed and compared to the base case for the increase in the operating cost. Hence, the loadability can be simulated only for ε = 0. As Table 3 shows, the maximum loadability in winter and summer disregarding any increase in the operating cost would be 1.019 and 0.308, respectively. In other words, the maximum loadability has been derived without violating any constraint. Using the IGDT technique, the system operator is looking for the maximum loadability constrained to the maximum limit set for the operating cost. It is obvious that by increasing the operating cost beyond the maximum loadability, no change would occur in the system loadability. Fig. 16 depicts the numerical comparison made between the chance-constrained programming and IGDT techniques. As expected, once the cost deviation is zero, the variations in the loadablity would also be zero. On the other hand, when the system operator accepts to tolerate higher operating costs, the total system loadability would also increase in line with the increase in the total operating cost. As Fig. 16(a) demonstrates, the maximum system loadability values by using the IGDT and chance-constrained programming techniques have the same trend. It is noted that the maximum loadability in summer is 0.308 corresponding to the 45.7% increase in the base case cost. On the other hand, the maximum loadability in winter is 1.019 for the chance-constrained programming and IGDT techniques. As can be observed in Fig. 16(b), the maximum loadability by using the IGDT technique by increasing the operating cost is higher than the chance-constrained programming. However, both methods have finally led to the same maximum loadability. It is also noteworthy that the maximum loadability in this case requires 119.26% increase in the base case cost.

V. CONCLUSION
This paper investigated the problem of optimal operation of the EH by using a robust chance-constrained approach. The optimal day-ahead operation problem was formulated by utilizing a MILP model. The Big M technique was used to transform the primary probabilistic optimization problem into a deterministic one that can be solved by the available commercial MILP solvers. Moreover, the robustness of the operation model against the uncertainties was assessed by employing a loadability index within the chance-constrained framework. Besides, a comparison was also made between the chance-constrained programming and the IGDT technique. It was indicated that elevating the chance-constrained model to the robust chance-constrained one would enable the system operator to implement the optimal operational strategy with respect to the risk index. The simulation results showed that the loadability index for the worst scenario in the winter was significantly higher than that of the summer. This is due to the fact that there were various assets available in the winter to supply the heating load demand, i.e. CHP unit, EH, and EHP. On the other hand, the options to serve the electrical and cooling loads, forming the major part of the total load demand, were limited in the summer. Besides, the simulation results for different risk indexes revealed that in the winter, first, the constraints relating to the electrical and then, heating power-related ones were violated. The results showed that there was no difficulty in serving the cooling load. It is noteworthy that as the loadability index in all cases was greater than 1, the uncertainty level should be more than 100% to provide the opportunity to serve the load demand without any shedding. The loadability index in the summer is limited due to the increased electrical load demand and the first violated constraint relates to the electrical load demand. It is noteworthy that the IGDT is indeed equivalent to the chance-constrained programming in case ε =0. In other words, the chance-constrained programming is a more complete model compared to the IGDT technique. In the Nordic region with extremely cold winters, the heating load demand is considerably high, and as a result, the total loadability index highly depends on the loadability for the heating load. Accordingly, the total loadability index is mainly dependent upon the electrical loadability in the summer. Therefore, the capacity expansion in the solar PV panel, battery, transformer, and CHP unit would help to reinforce the EH.