An Optimal Energy Hub Management Integrated EVs and RES Based on Three-Stage Model Considering Various Uncertainties

In order to make significant progress in the operation of power systems and minimizing cost and air pollution, multi-carrier energy systems based on the economic and emission problem have been proposed by implementing various solutions and using different sources. In this study, an energy hub (EH) system includes a wind turbine (WT), photovoltaic (PV), and EVs that are exchanged energy with energy and reserve market considering energy, thermal, and gas demand response program is proposed. But these units due to uncertain behavior create a big problem in balancing between the demand and generated energy. Thus, the uncertainty of WT, PV, load, and electricity market price are modeled with Mont-Carlo method. Besides that, all parameters of the EVs with uncertainty behavior are modeled using a new method that classified EVs based on the capacity of the battery and other features. This method is employed a stochastic optimization approach to simplify the uncertainty modeling for increasing the system reliability. Hence, two objective functions, namely economic cost, and environmental cost are considered. Finally, a three-step strategy is introduced to solve the multi-objective problem as a single objective function. Finally, EH management performance is investigated by implementing the proposed method and elements.


INDEX TERMS
In recent years, in the power system structure, a multi-carrier power system has been implemented instead of the single entity generating [1], [2]. On the other hand, increasing the emission pollution problem due to the use of fossil fuels is one of the serious reasons [3]. Therefore, energy hubs are supposed to play a critical role in the future for solving power system problems [4]. RESs such as PV and wind farms generate clean and unlimited energy [5], and the use of renewable energy is an effective way to promote energy transformation in the energy hub [6]. Today, the use of electric vehicles is increasing due to reducing air pollution and exploitation as consumers and producers (as a battery) [7]. EVs are also profitable for car owners, in addition to the economic and environmental benefits of the power system [8]. Therefore, integrating EVs and energy hubs is a good choice for minimizing operation costs and air pollution. Renewable resources have uncertain behaviors due to climate change, which has a significant impact on the performance and reliability of the EH. On the other hand, EVs include uncertainty behaviors [9]. The most suitable operation of the EH will be obtained using a suitable optimization method and integrating the energy market with the reserve market [10] considering uncertainties in the presence of EVs.

B. LITERATURE REVIEW
Different studies of assessing and analyzing EHs abilities due to many advantages such as economy minimization, emission pollution mitigation, high flexibility, market infrastructure, and energy management in various sections of power systems have been done [2], [11], [12]. EH is made from two parts as an input and output, which contain multi carriers in the input section and multi demands in the output section. Employing renewable resources, suitable converters, energy storages, and redundant connection between EH input and output, supply of demand, reliability, and energy management are enhanced, the operation cost and emission pollution are minimized [13]. Implementing a good framework of EH and an appropriate scheduling program is essential to obtain the best results for EH operation. Alessandra Parisio et al. proposed a robust optimization method of EH operations [14]. Also, multi carriers' input, stored and distributed in the EH construction is implemented to supply energy demands. However, renewable sources are a good option for input in the EH framework.
In [15], EH as a combined energy system containing renewable energy resources was proposed considering demand response programs to supply electrical, cooling, and heating demands. But, the stochastic nature in the production of the RESRs has not been considered. In [16], a new mathematical method for analyzing the integrated EH with a deregulated market was reported. This method is based on a mathematic program with an equilibrium constraints model to investigate the strategic performances of a profit-driven EH in both the heating and electricity markets. Since the economic issue is important in EH performance, at the same rate, the pollution emission is important, too. So, Lu et al. in [17], an optimal load program model for a community EH represented. Both EDRP and TDRP to minimize the total cost of the energy hub, including system economic and environmental costs were considered. One of the most effective methods to improve the performance of EH and optimize the operation cost and emission pollution is using an appropriate energy storage system. In [18], and EH management mechanism was introduced by employing two new energy storages are called ice storage conditioner (ISC) and solar-powered compressed air energy storage to minimize operation and emission pollution costs considering combined CHP, photovoltaic, and wind turbine. The other important factor, which has the most effect on EH performance is using a robust scheduling method. Shams et al. [19] suggested a min-max min robust method for the short-term scheduling of microgrids integrated with the natural gas network to model the uncertainty of WT generation, electrical, and thermal loads introduced. The proposed model was solved and linearized using the column-and-constraint generation process that divides the framework into a master problem and a subproblem. It should be better if the improved EH structure behind a suitable scheduling method was analyzed simultaneously. In [20], novel energy storage integrated with an energy hub in an industrial park was proposed and a bi-level optimization scheduling method was introduced to analyze the optimal bidding and programing in the daily market. The main aim of this method is to minimize the day-ahead operation cost.
Since renewable sources are subject to weather conditions, their output power is not constant and has uncertainty. Therefore, it is essential to study the causes of uncertainty with appropriate methods. In [21], an MINLP method for the optimal scheduling of the energy hub was modeled in the presence of uncertainty factors. In addition, the effect of implementing DRPs, determine the size of the energy hub, and implementing a complete energy flow of loads are introduced.
Using the risk-averse method is a good option for modeling the uncertainty of the energy hub to increase reliability. A novel look-ahead risk-constrained based on the look-ahead scheduling technique in the daily market of the energy hub was proposed to reduce the total cost of the energy hub in the presence of DRP in [22]. Mokaramian et al. [23], analyzed the performance of the energy hub using a multi-objective function as environmental, economic, load deviation, and risk level objective functions conditional value at risk and secondorder stochastic dominance. Also, a new storage model as PHESS has been utilized in EH to optimal energy management considering DRP for water, gas, electricity, and thermal.
In [24], a spatiotemporal association of wind speed, solar irradiation, and price responsive loads in the presence of cooperative and non-cooperative modes of energy hubs using a game-theoretic method is proposed, which has participation by the energy market. In addition, mathematical scheduling with equilibrium restraints was applied to solve the highly competitive behavior of the energy hubs in this market in the presence of system losses. Stochastic operation of energy hub in the uncertain environment and using downside risk restrictions was suggested to minimize the risk-in-cost versus risk control parameter in the heating market with heat DRP in [25]. In [26], a multi-objective method was represented to optimize the operation cost and risk-averse in the energy hub. The conditional value at risk method was implemented to control the harmful effects of the uncertainties.
In [27], a multi-agent system is introduced to model the operation of an EH with electric vehicles and modeled various penetration rates and charging patterns of an electric vehicle. A complete random dispatch algorithm was integrated into a smart charging strategy to obtain the maximum capacity and potential of the vehicle to the grid. XinhuiLu et al. [28] presented a novel structure of optimal load dispatch model for a community EH in the presence of thermal and electrical demand response programs. In this model, the gas boiler, CHP unit, PV array, heat storage unit, WT, the random admission of large-scale electric vehicles was considered, and using the Monte-Carlo approach, the uncertainty of EVs is simulated. Also, a robust optimization method was proposed to simulate the price uncertainty. To improve the environmental and economic aspect of the energy hub, a multi-objective optimization method in the presence of compressed air energy storage was integrated with a battery energy storage system in [29]. Also, in the structure of residential EH, a plug-in electric vehicle-based demand response program, renewable energy generation units, thermal energy storage, solar heat collectors, and hot water storage was used. Finally, the load demand uncertainty of EVs is modeled using their corresponding coming-exit time, miles of daily moved, and type of vehicle. The optimal Pareto front of the multi-objective issue was achieved by the ε-constrained approach for different robustness alterations.
Regarding the reviewed literature in EH, some of the drawbacks have seemed that should be improved as follow: • In some studies, renewable sources, fuel cells, and EVs as energy storage have not been used to improve environmental and economic issues.
• Different DRP for water, gas, thermal, and electrical loads have not been considered.
• The impact of all parameters of EVs uncertainties on the operation of the energy hub has not been esteemed.
• Uncertainty of renewable sources, load uncertainty, market price, and EVs have not been considered, simultaneously.

C. CONTRIBUTION
In this study, two objectives are considered to optimize the total cost of EH performance as emission pollution and operation cost. The proposed model of EH includes a CHP, boiler, WT, PV, Fuel cell, ESS, HSS, and EV considering the uncertainty of WT, PV, load, price market, and EVs as shown in Fig. 1. In addition, a demand response program is implemented for gas, heat, and electricity demand. All uncertainty factors of EVs as starting and stopping time of charging/discharging period, the overall battery capacity, initial energy, and their number, are not yet addressed for an EH or a collector that uses these resources integrated with renewables to contribute to different electricity markets. So, in this study, we introduce a novel method to implement the storage capacity of EVs as energy storage in EH structure in the day-ahead market. Here, three factors of EVs uncertainties are considered contains EVs capacity uncertainty, coming and leaving time, the amount of energy left in the battery when it becomes obtainable, and their total number.
Finally, the fundamental contributions of this study are mentioned as follow: • Participating EH with energy and reserve markets • Integrating the EVs with day-ahead and reserve markets • Proposing a new method to model all uncertainty parameters of EVs.
• Considering demand response programs for gas, electricity, and thermal loads • Solving the multi-objective problem using a novel threestage solution algorithm.

II. ENERGY HUB MODELING
This study presents an EH including different energy carriers with the ability to convert, storing and transferring. The energy hub as shown in Fig. 1 is including utility grid, natural gas network, fuel cell, and renewable sources (wind turbine, PV) in the input section that integrated with ESS and TSS as a storage system, CHP, and boiler in converter section to supply load demand of gas, electricity, heat, and EVs that operate as a mobile storage. The main aim of this study is to minimize operating costs and emission pollution problems. Also, the uncertainty of RES, load market, price of the energy market, and three factors of EVs uncertainties i.e.: EVs capacity uncertainty, coming and leaving time, and their total number are modeled using the Mont-Carlo method. Noted that in this study, EH is assumed as a small city. Electricity and natural gas as input carriers are considered to supply the load demand. Wind turbines and PV can produce electrical energy to use by the consumer. The EVs can operate as bilateral equipment and in the charging mode, is assumed as a consumer, and in discharging mode, the extra power is saved and flow to the EH to use via consumer, in this mode the EVs work as energy storage. The CHP unit in the EH framework is fed from natural gas and produces heat and electrical energy simultaneously. In addition, the auxiliary boiler feeds the natural gas and generates heat to the heating line.
The CHP unit has an important role in EH performance because this unit involved three carriers of EH as natural gas, heat, and electricity. In this study, the CHP unit works in quadratic mode and contains three operational zones, which are shown in Fig. 2.
The relation between the heat and power in CHP operation as shown in Fig. 2, is a non-convex operation zone. When the thermal is produced by the CHP unit, the amount of power generation is limited and conversely. In order to simplify the optimal operation of the CHP unit, the non-convex zone is divided into triple zones which are shown by Zone 1, 2, and 3. Convexness is made by counting three binary variables for the CHP. These zones are convex, and any point on any line in each zone is related to that zone [30].

A. OBJECTIVE FUNCTIONS
The main aims of this study are divided into two parts. The first aim is to minimize total operation cost and the second aim is to optimize the emission pollution cost of the EH in the day-ahead market.

B. ECONOMIC OBJECTIVE FUNCTION
The economic objective function includes two terms to optimize the total cost of the EH as follow [23]: The first term denotes the operation cost of the EH, and the second term represents the reserve market cost. The energy hub system contains RESs, boiler, CHP unit, thermal and electrical energy storage, and EVs. Demand response program is considered for all types of load demand including gas, thermal and electric demand with a specified cost. So, the operation cost of EH is presented as follow: The energy hub system purchase energy from the upstream grid of electric and gas station, which is calculated as follow: In order to manage energy usage in peak load duration by shifting the energy consumption to low periods, DRP is implemented in the energy hub to cost-minimizing. The main effect of the demand response program is based on the balance between energy consumption and production in each period. Electrical, thermal and gas demand program are expressed as follow: Also, the operation cost of the storage system is made from two-part including the electric storage system and thermal storage system as follows.
The operational cost of the renewable sources is equal to the set of operation costs of the WT, PV, and an FC that is a function of output power.
The operation cost of the CHP unit is a quadratic function of the generated power and heat and given by.
The operation cost of the boiler is a function of output power and can be defined as follow: In order to provide a suitable economic program for electric vehicles, the objective function is defined based on the amount of charging/discharging energy of electric vehicles. The required energy for charging EVs is provided by EH, P dem j , with the specific cost, and at the same time when EVs have been used as an energy storage system at time t with specific degradation cost.
which is calculated as follow: The second term of Eq. (1) is the cost of purchasing energy from the reserve market and is formulated by Eq. (10).
The environmental issue of the energy hub is critical. So, the emission pollution cost of the three most effective gases (CO2, NO2, and SO2) which are generated through EH equipment is considered in the objective function. The emission cost of EH is produced by the upstream grid, boiler, CHP unit, and consumed residential gas [25].
The constraint of up and low load levels of the EH is formulated by Eq. (13). The increasing and decreasing load level cannot be done at the same time which is modeled in Eq. (14). The EH must save the overall electricity consumption constant in the forecast cycle, which is restricted by Eq. (15).

2) LIMITATION OF THERMAL DR PROGRAM
The equal and unequal limitations of thermal DRP are as follow: Similar to the EDRP, Eq. (17) formulates the limitation of increased/decreased load of the energy hub. The restraint to cope with the increasing/decreasing of the thermal load level at the same time is modeled by Eq. (18). The amount of increasing/decreasing thermal load should be equal which is formulated by Eq. (19).

3) LIMITATION OF GAS DR PROGRAM
The equal and unequal constraints of gas DRP are as follow: The exact definitions of the previous sections apply to GDRP.

4) LIMITATION OF ESS
Another important factor in the EH framework is to save the extra produced power by units. So, ESS is used to deal with energy losses. The mathematical formulation and limitation of ESS are stated as [20]. The electrical energy balance restraint of the ESS unit is given by Eq. (24). The saved electrical energy by the ESS unit should be restricted in a confident range, which is given by Eq. (25). The maximum charging/discharging powers of the ESS unit are restricted by Eqs. (26) and (27)

5) LIMITATION OF TSS
Thermal demand in the energy hub is at a high level, so it is necessary to save the extra produced power by units. Thus, TSS is used to deal with energy losses. The equal and unequal constraints of TSS are [20].
The exact definitions of the previous sections are implemented to TSS.

6) ENERGY BALANCE CONSTRAINT
The produced or bought energy from the upstream grid in alltime should be equal to the demand for energy. So, electrical, thermal, and gas balance constraint are expressed as follow, respectively.
A nominal rate is definitely for purchasing power from the upstream grid that should supply transformers nominal capacity constraint and total bought gas limitation from the gas grid.

8) MODEL OF BOILER AND CHP
The boiler is fed by natural gas and it is necessary to meet the upper and lower level of natural gas input constraints [30].
As stated in the previous sections, the operation of CHP is a nonlinear behavior with the triple-zone of electricity and thermal production. The generated heat through the CHP unit is frequently proportional to its electric generation. In addition, heat operation constraints will be guaranteed if the electrical performance constraint is confirmed [31]. Limitations associated with the CHP unit with three operational areas as shown VOLUME 10, 2022  in Fig. 4 are formulated as follows: Using the big-M method, the upper and lower range are defined viaM and M , that represent inactive at the predicted time [30].

9) SOURCES MODEL
• Wind turbine The produced power by the wind turbine is proportional to wind speed and is given by [23].
• Photovoltaic The imported power by the PV system is proportional to the solar irradiance as follow [32]: • Fuel Cell The power made through the FC in each period t∈ {1, 2, · · · , NT}, H FC 2t , which is proportional with the hydrogen feeding in peck period is expressed as [32]: Several constraints related to fuel cell performance are:

10) EVs CONSTRAINTS
The constraints of EVs participation in EH operation and their real-life restrictions are defined as: Eq. (52) represents the final battery power of the EV when the plug-in time is ended, which is equal to the total of the demanded energy through j_th EV and the initial power of the EV battery [33]. The plug-in time of EV is in the interval of [t b,j , t e,j ], and the demanded energy of the EV can be defined by the difference between needed energy for charging the battery and used energy by EVs in the day ahead and reserve markets.
For charging/discharging an EV battery, a certain bound is considered equal to the battery conditions and existing equipment minus the charging and discharging facilities. Also, the charging and discharging operation cannot be done at the same time.
The amount of energy in EVs batteries at the arrival and departure time is formulated as follow: The power of the j_th EV at time t is related to three factors including the initial state of the battery, charging/discharging power of the EV at the same time. The starting and stopping time of the charging and discharging of the battery is also formulated as follow:

E. STOCHASTIC SIMULATION FOR MODELING THE UNCERTAINTY
Some of the energy carriers' parameters in the EH have random and vague behavior. Therefore, it is hard to predict this parameter's behavior. Finding an appropriate method to forecast these parameters' function with high reliability is significant.
The output power of WT and PV is related to climate changes, and their forecasted power includes uncertainty. So, using a stochastic model can try to achieve the most exact model. Moreover, the price of the energy market and load demand are affected by consumer demand and climate behavior. Like WT and PV, it is essential to implement an effective stochastic method to model the uncertainty of these parameters.
Finally, for the first time in the EH performance, all of the EV parameters that have uncertain behavior consist of coming and leaving time, the energy remaining in the battery when it becomes available, and the number of EH is considered.
A scenario-based stochastic model is implemented to model the uncertainty behavior of all parameters.

F. SCENARIO GENERATION
In this model, a PDF of each parameter is used and shown in Fig. 3. Scenarios of the PDF curve are divided into seven levels. Each is of PDF contains a probability rate. It should be noted that the center of the fourth scenario is the same predicted value.

1) SOLAR IRRADIANCE UNCERTAINTY
Since the output power of the solar is a function of the solar irradiance and weather variation. Beta PDF is used to model the uncertainty of solar irradiance. Beta PDF is a function of two parameters as ω and ψ [35]: where the ψ and ω variables are computed as follow: 2) WIND SPEED UNCERTAINTY Similar to the solar operation, wind speed is a function of weather changes. So, the Weibull distribution function is implemented to model the wind speed uncertainty [36].

3) UNCERTAINTY OF MARKET PRICE AND LOAD
Since the exact value of consumer demands and energy market price is not definite. A normal PDF function is used to model the load and market price uncertainty.

4) EVs UNCERTAINTY
In this study, EVs operate in a dual state as both consumer and energy storage. The EV has three factors of uncertainty behavior. The saved energy in the battery at the start, their availability, and needed energy is unclear when the battery is in the charge or discharge state. Similar to another uncertainty parameter, EH can detect the charging state and behavior of EVs in the parking lots.
In order to simplify the uncertainty modeling of the EVs, like the [37], EVs are categorized into different types. Each type of electric vehicle has its specifications. These futures include the obtainability time, initial state of the battery, and VOLUME 10, 2022 capacity of the battery. In this method, 1000 different EVs are classified into 100 types. So, in addition to the small number of EVs and uncertainty parameters are studied, the EVs' parameters are linked together, characterizing them causes the chance to analyze the uncertainty in a more straightforward and improved way.
Each type of EV in the EH has grouped base on the {Eu,0, [tb,u, te,u], Emax,u} by assuming that u ∈U.
Implementing this method is caused to reduce the sources of the uncertainty parameter to the number of EVs. Here, the number of EVs from each group is assumed as an uncertain parameter.
The uncertainty of the EVs with w scenarios and type u are formulated as follow: To model the uncertainty of EVs numbers, the normal distribution method with zero mean and σ n as the distribution error is implemented.

G. K-MEANS CLUSTERING METHOD
In this paper, five uncertainty parameters are modeled. The PDF is divided into nine levels, and each level represents various scenarios and errors. In the first step, 1000 scenarios for each parameter are generated with the mentioned method and using the k-means clustering approach the scenario numbers are reduced to 10 scenarios in the MATLAB environment.

H. MULTI-STAGE SOLUTION METHOD
In this study, two functions as economic and emission pollution objectives are considered. A multi-stage solution is proposed that consists of three stages [38].
First stage: payoff table establishment Second stage: fuzzy linearization Third stage: Objective weight calculation.

1) PAYOFF TABLE ESTABLISHMENT
The payoff table is implemented to investigate the relationship between the various objective functions. This table is used to calculate the weight coefficients. The stages of the establishment are as follow: The objective function is chosen as the optimization objective obji (i = 1, 2 . . . , I) where the i_th objective function value is objii and another objective functions values are objik(i, k = 1, 2, . . . , I).
The pay of table is structured as Table.1 based on the obtained value of the objective functions. Then the lower and upper value of the objectives is used for objective weight competition and fuzzy linearization.
The fuzzy linearization and objective value competition are described in the next section.

2) FUZZY LINEARIZATION
In this study, the descending half-line membership function is used to process the objective functions of the minimum oper- ation cost and the minimum carbon emissions. The specific method is as follows: Eq. (63) represents the half-life reducing half gradient membership function and is utilized to procedure the objective functions of the minimum operation cost, which is shown in Fig. 4.

3) OBJECTIVE WEIGHT CALCULATION
The objective authorization approach is applied to compute the weight coefficients. The pre-processed objective function decision matrix [ob ji k]I×I is calculated by Table 1. Then, to compute the weight of the objective function, the entropy weight approach is used by implementing the different degrees of the objective function. So, the weights of any objective could be calculated by the entropy weight approach [39] as: where ϕ = 1/ln(n) is the constant represented with the sample number, where Ei ∈ [0, 1]. r ij meets 0 < r ij < 1, and r ij = 1 when rij = 0, rijln(rij) = 0.
• Compute the weights of the objective function obji as follows: Based on the weight coefficients of various objective functions, the total single objective function can be formulated as follows: where, OBJ represents the total single objective function, which can achieve the best acceptable solution for matching all objective functions.

III. NUMERICAL RESULTS AND SIMULATION
In this section, the optimum scheduling of the proposed EH as shown in Fig. 1 are studied in the energy and reserve markets. It includes renewable sources as WT, PV, and fuel cell, storages as ESS and TSS, and EVs in both states of storage and generator.

A. INPUT DATA
In this part, all the essential data for implementing the shortterm optimized scheduling of the proposed EH is presented. The main aim of the proposed structure is to minimize the economic and emission pollution costs using participating in the day-ahead energy and reserve markets. All of the needed data are categorized in two sections a real-time data and forecasted data. The needed forecasted data as WT speed, solar irradiance, load demand, market price, and number of EVs based on the real-time data are made using the Mont Carlo approach in 600 scenarios, and the by k-mean cluster-   ing method these data are decreased to ten scenarios. In Fig. 6, wind turbine speed scenarios with the red line, and 10 reduced scenarios with different color is plotted. Similar to the wind turbine, other mentioned uncertainty parameters are shown in Figs.7-9.
The last uncertain parameter is the EVs number of each type of electric vehicle based on the battery capacity in the interval of [5 90] kwh and also charging/discharging enabling which have enhanced all time significantly. Some of the parameters used to classify electric vehicles are analyzed in [40]. The obtained necessary data are integrated with the proposed method of scenario generation. 11 types of popular  EVs, based on the availability, capacity, and another critical behavior parameter, in the German market with sales statices are introduced in Table 2 [40]. Moreover, their percentage is used as a range in [0 1]. In every iteration, a random number among the 0 and 1 states which type of EVs is chosen. Based on the total capacity of parking lots, the number of each type of EV is produced. The number of different types of EVs for the day-ahead in each scenario, the overall number of EVs in each scenario, and the reduced 10 number scenarios are depicted in Fig. 10. Fig. 11 shows the energy, gas, and heat demand of the EH (kW). Table 3 lists the essential parameters of TSS, ESS, and WT. The real-time energy demand and reserve market price are exposed in Fig. 12. The data associated with the various units in the energy hub (CHP, boiler, gas, grid, PV, etc.) is given in Table 4. The emission factor of NO2, SO2, and CO2 in each unit of EH in kg/kWh is presented in Table 5. In Fig. 13, the required parameters of the CHP unit are defined. The operation cost of numerous grids and units are given in Table 6.
Finally, the proposed EH model with stochastic constraints has a nonlinear behavior that is modeled using the equivalent MINLP method and is solved with CPLEX in GAMS.
In order to assess the proposed method, three cases are classified as follow: Case 1: Simple EH without EVs and reserve market Case 2: Considering the EVs and uncertainty factors Case3: Considering EH with EVs and reserve market.

Case 1:
In this case, the simple structure of the EH is assumed and both objective functions are minimized using the multi-stage solution method. At first, to investigate the   difference between novel multi-stage method and other methods presented in [31], the results are shown in Table 7: So, applying a new method, in addition to responding quickly and easily, reduces the operation and emission pollution costs. The total operation cost of EH without and with EDRP, TDRP, and GDRP is $2581.700, $2487.808, respectively. Also, the emission pollution of EH without and with EDRP, TDRP, and GDRP are 10190.837kg, 10089.064kg,    respectively. It is obvious that implementing DRP is caused to reduce both objective functions. The purchasing power from the upstream and gas network are shown in Fig. 14. Also, the output power of renewable sources is shown in Fig. 15.
Case 2: In this case, the EVs are integrated with EH, and the uncertainty factor of solar irradiance, wind speed, market price, load, and EVs are applied. The total number of EVs is equal to 21. The operation cost and emission pollution values are $ 2109.250 and 8139.657 kg, respectively which is decreased about 15.2% and 19.32% rather than case 1, respectively. It is clear that implementing the EVs is caused to reduce both objective functions, as EVs don't produce   emission pollution, this cost is significantly reduced but the most effective of EVs is on the emission pollution issue due these devices do not produce any pollution. In Fig. 16, produced power by different units of EH is shown. Also, the purchased gas from CHP, boiler, and gas network are plotted in Fig. 17. Finally, the charged and discharge power of EVs,  ESS, and TSS is represented in Fig. 18. Most of the required energy for EH is provided by EVs because the emission and cost of EVs are appropriate for the ESS and TSS.
Charge and discharge power of total EVs, ESS, and TSS.  It is clear from Fig.18 that part of the rechargeable energy of electric vehicles is used to supply the demanded energy.
Case 3: Finally, in this case, EH is participating with the reserve market for energy management. In this case, when the load demand is high, EH tries to buy energy from the reserve market at a low price. So, the amount of buying energy from the upstream network with high costs is reduced. In this case, the operation cost and emission pollution value are $1638.577 and 7642.973 kg, respectively. The operation cost of case 3 relative to the case 2 and 1 is decreased about 22.31% and 34.1%. Also, the emission pollution value of case 3 related to case 2 and case 1 is reduced by about 6.1% and 24.2%. Therefore, it can be concluded that the reserve market has a greater impact on operating costs. In Fig. 19, day-ahead forecasting cost of energy and reserve markets bought energy via upstream grid and reserve market is shown. It is clear that applying the reserve market is caused to flat the profile of load.
Finally, the operation costs and emission pollution costs of the 3 cases are listed in Table. 8. As expected, implementing EVs is caused to the reduction of both objective functions relative to case 1, but the second objective is significantly reduced. Finally, participating in the reserve market is a good choice for energy management in peak time and also operation cost minimization.

IV. CONCLUSION
In this study, the short-term optimal scheduling of EH as a small city using a stochastic optimization method was applied with two objective functions in the energy and reserve market. In EH, three storages consisting of TSS, ESS, and EVs in discharging mode, were used. Renewable and nonrenewable resources such as WT, PV, FC, CHP, upstream grid, and natural gas station were considered. The economic operation cost and air pollution problem were considered as a bi-level objective considering EDRP, TDRP, and GDRP, to obtain the EH optimal scheduling in an uncertain situation. The wind speed, solar irradiance, loads, market price, and the number of EVs were assumed as uncertain parameters. The multi-stage solution method was applied to speed up system performance and simplify procedures of the bi-objective optimization problem. One of the main aids of this study is integrating the EVs with the EH considering the three uncertainties of EVs. This three-parameter is converted to one uncertainty parameter which is the number of EVs. Also, due to the uncertain behavior of some elements, implementing appropriate storage is a good option.
Three test cases are defined to analyze the effectiveness of the proposed method with considering reserve market and EVs. Case 1 was considered as a base case of the EH without EVs and reserve market, however, in this case, the effect of applying different DRPs was investigated. In case 2, EVs were integrated with the EH and both modes of consumer and producer are applied. Two objective functions (economic and environmental cost) were decreased by about 5.34% and 3.54%, compared to Case 1, respectively. Finally, the energy hub has traded with the energy and reserve markets for bidding and sailing energy. In this case, the operation and pollution emission costs were decreased compared to Case 2 by about 6.27%, and 8.63%, respectively. It can be concluded that the operation cost and pollution emission of the EH is highly minimized by implementing the proposed three-stage method and reserve market participation in an uncertain environment. AMIN SAFARI received the B.Sc. degree in electrical engineering from the University of Tabriz, Tabriz, Iran, in 2007, the M.Sc. degree from the University of Zanjan, in 2009, and the Ph.D. degree from the Iran University of Science and Technology, Iran, in 2013. He is currently an Assistant Professor with the Department of Electrical Engineering, Azarbaijan Shahid Madani University, Tabriz. He has published more than 140 articles in international journals and conference proceedings. His research interests include the application of artificial intelligence and heuristic optimization algorithms to power system design, FACTS device, power system analysis and control, renewable energy, control and management of microgrids, and smart grid. He was selected as a Distinguished Researcher of the Azarbaijan Shahid Madani University, in 2016.
HASSAN HAES ALHELOU (Senior Member, IEEE) is currently a Faculty Member with Tisheen University, Lattakia, Syria. He is included in the 2018 and 2019 Publons list of the top 1% best reviewer and researchers in the field of engineering. He has published more than 100 research articles in the high-quality peer-reviewed journals and international conferences. He has participated in more than 15 industrial projects. His current research interests include power systems, power system dynamics, power system operation and control, dynamic state estimation, frequency control, smart grids, micro grids, demand response, load shedding, and power system protection. He was a recipient of the Outstanding Reviewer Award from Energy Conversion and Management journal, in 2016, ISA Transactions journal, in 2018, Applied Energy journal, in 2019, and many other awards. He was also a recipient of the Best Young Researcher in the Arab Student Forum Creative among 61 researchers from 16 countries at Alexandria University, Egypt, in 2011. He has also performed reviews for high prestigious journals, including the IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS, the IEEE TRANSACTIONS ON INDUSTRIAL ELECTRONICS, Energy Conversion and Management, Applied Energy, and the International Journal of Electrical Power and Energy Systems. VOLUME 10, 2022