Bi‐level optimization of the integrated energy systems in the deregulated energy markets considering the prediction of uncertain parameters and price‐based demand response program

The proliferation of multienergy systems (MESs) in recent years has led to the higher efficiency of energy consumption in different sectors. MESs are generally studied in the context of energy hubs (EHs). Most of the previous works in the field of EH operation and scheduling are at the scale of residential households or communities. There is a research gap in proposing a framework in which the EH actively interacts with different energy networks. This article, however, is focused on the optimal operation and scheduling of the EH and two integrated energy grids such as power and heat distribution systems. The EH consists of combined heat and power (CHP), electric heat pump (EHP), thermal and electric energy storage systems, which are operated by an independent entity. The proposed optimization problem is based on the bi‐level optimization method in which the EH is the leader and two networks are the followers. Nonlinear variables are also linearized using binary expansion and Big M methods. As another research novelty, uncertainties of load demand and photovoltaic (PV) systems' output power are taken into account using a mixed machine learning (ML)‐based forecasting method. To enhance the flexibility of the system, a price‐based demand response (DR) program is proposed based on the predicted load pattern of the mixed ML method. Simulation results show the efficiency of the optimization model. Also, the DR program has a significant impact on the system's cost by 40% improvement in the profit of the system.


| INTRODUCTION
In recent years, the cooperation of different generation systems has been facilitated thanks to the development of new technologies. [1][2][3] This trend has led to the proliferation of multienergy systems (MESs) in the current energy infrastructures. 4,5 Generally speaking, MESs incorporate different energy vectors such as electricity, district heating, and natural gas into the energy infrastructures. In fact, the crucial concept of the MESs is in the operational interactions of energy (vectors) carriers utilized in the system. 6,7 In contrast to traditional systems where a single vector of energy is analyzed (take e.g., power systems operation), in the MESs, two or more types of energy carriers are usually investigated, making the system structure more complex and computationally demanding. 8,9 To tackle the uprising problems, a suitable tool as the energy hub (EH) is proposed, which provides appropriate mathematical principles and services to convert, store, trade, and/or consume various energy carriers. The providing flexibility of EHs is significant because of the benefits of different devices and facilities, namely, combined heat and power (CHP) systems, gas/ electric boilers, electrical and thermal energy storage systems, and electric heat pumps (EHPs). Besides, these EHs typically advantage from primary energies such as electricity, heating, and natural gas. 10,11 Availability and synergic effects of the resources positively support decision-makers in the energy markets with more options for making optimal decisions. [12][13][14] The concept of MESs and EHs have been comprehensively investigated in previous years. In Sadeghi et al., 15 a systematic review is performed in the area of EHs operation and expansion planning problems. Another review paper in Mohammadi et al. 16 also researched the dominant structures used in model EHs. In addition, challenges, strengths, and weaknesses of EHs in recent studies were identified thoroughly. Authors of Mansour-Saatloo et al. 17 addressed the optimal self-scheduling of hydrogen-based EH considering demand response (DR) program and fuel cell-based hydrogen storage system. In a similar study, 5 researchers developed a self-scheduling EH model to consider the day-ahead and regulation markets. The study is modeled as stochastic mixed-integer linear programming (MILP) to handle the uncertainties of regulation market prices, wind turbine (WT) production, and electric vehicles patterns. In Shahrabi et al., 18 optimal operation and planning of the EH are performed using Quantum Particle Swarm Optimization (QPSO) algorithm. Environmental aspects of the system are also considered in this study by minimizing fuel consumption and emission production. In Miao et al., 19 compromise programming (CO) is proposed to optimize multiobjective problems considering the economic and environmental issues of a residential EH. They used CO to handle the problem. In Mokaramian et al., 20 a similar study is conducted for a residential EH, which also considered the conditional value at risk (CVaR) approach to enhance the EH operation and efficiency. In Faraji et al., 12 the authors researched on proposing a residential EH to address electrical and heat loads. The authors responded to renewable uncertainties and considered the equipment reliability of EH-based N-1 contingencies.
Most of the mentioned studies are conducted on a residential scale, which has a minor impact on the decisions of the upstream grid. Nevertheless, it is of great importance to analyze and perform the studies from the distribution-level perspective. There are limited studies in this area; for instance, the authors of Jadidbonab et al. 21 investigated the interaction among MESs and EHs thermal and electricity markets. To handle the uncertainty, they implemented the information gap theory (IGTD) method. In another research, 22 power and thermal distribution networks are studied from the view of a compressed air storage system (CAES) as the EH. In Arabkoohsar and Andresen, 23 energy storage is modeled that is able to produce both electricity and thermal energy. The proposed energy storage participates in both heat and electricity network markets. However, due to the limited operation of the energy storage, it cannot be considered an EH such that introduced in the literature. The aim of the study was to minimize the power purchased from the upstream grid and heat pumps power demand in the thermal network-both systems proposed in Jadidbonab et al. and Hu et al. 21,22 are managed and operated by a central unit. This approach may not be efficient in real-world practices since both thermal and electrical networks are operated by two individual entities. In Li et al., 24 a coordinated model is proposed for heat and electricity markets, which responds to the decisions made by the EH operator. However, the most crucial research gap in Li et al., 24 is that the uncertainty of load consumption had not been considered. Najafi et al. 25 presented a stochastic bilevel problem for the decision-making of an EH manager, which maximizes its profit by selling electricity and heat. On the other hand, the clients minimize their total cost of energy procurement. In addition, CVaR is used to decrease the unfavorable impacts of the uncertainties. In addition, there are some studies that are focused on the planning and sizing of EH equipment. In Zou et al., 26 an EH planning problem is proposed in which the upper level optimized the capacity of energy equipment, and the lower level optimized the combination of each energy carrier. In Ghasemi et al., 27 a multistage planning is proposed to determine the investment candidates with the best capacity for the equipment of integrated system. To solve the problem, the authors developed an MILP model and bi-level decomposition approach. In Najafi et al., 6 authors developed a medium-term management of EH in restructured power systems. Table 1 shows the taxonomy of recent studies in the area of EH operation and scheduling.
In this study, we tried to consider the fact that thermal and electricity networks are managed by two different operators, that is, power distribution network operator and thermal distribution network operator. The proposed EH is also considered as an independent unit controlled by a third-party entity that is able to submit its offers and bids to the thermal/electricity markets. The EH earns revenues by proposing "energy packs," that is, quantities and prices of electric and heat energy submitted in the corresponding markets. The proposed model is based on bi-level optimization method in which the EH is the leader (upper level), and the networks are followers (lower level). Moreover, to address the uncertainties of load and PV outputs, this article considered a mixed ML-based forecasting method. Using the predicted load values from the mixed ML-based forecasting method, a DR program based on an RTP scheme is also proposed in this article to improve the load profile of the power distribution system. Two case studies are suggested according to the implementation of the DR program. The optimization results of this study show the efficiency of the proposed model. The list of contributions of this article are as follows: • Proposing a bi-level optimization model for the EH dispatch problem in the presence of power and heat networks. • Considering prediction of uncertain parameters using an ML-based forecasting method. • Developing a price-based DR program for active load demand.
The rest of the article is as follows. In Section 2, the EH and markets structure are introduced. In Section 3, optimization model including bi-level model of EH and two markets are provided. The forecasting methodology and DR program are also described in this section. In Section 4, prediction and optimization results are presented. Section 5 also summarizes the main findings of the study.

| EH AND MARKETS STRUCTURE
An EH can be considered as a black box consisting of input and output vectors. 28 In real practice, it is possible to build an EH by integrating energy conversion units like CHP systems, EHPs, thermal and electric storages, and boilers. EHs also include other equipment and resources, and their utilization depends on the purpose of the study. The schematic of the proposed EH is described in Figure 1. The suggested EH is able to link the power network to the heat network. The input gas of the EH is bought from a gas facility, and electricity is purchased from the power network. The distributionlevel EH can supply thermal and electrical energy to heat and power networks, respectively. Nevertheless, the residential-level EHs generally use the provided energy for their self-consumption. 24,29 As shown in Figure 1, the input electrical energy is consumed for supplying the EHP and storing it in the electric storage system. The EHP is clean technology for producing heat energy from electrical energy. 30 The CHP consumes natural gas to produce heat and electricity, which can be stored in thermal energy storage. As mentioned above, the EH is connected to both power and thermal distribution networks. This has been schematically illustrated in Figure 2. The EH submits energy packs (quantities and prices) to heat and electric markets and subsequently earns revenues according to the submitted values. The energy packs are submitted according to the anticipation of market-clearing results. The revenues of the EH are obtained by consulting a payas-bid agreement. As indicated in Figure 2, the two independent operators, that is, electricity distribution system operator (DSO) and heat market operator, are responsible for clearing the corresponding markets to minimize total generation cost. Moreover, due to the storage and energy conversion capabilities of the two markets, the bidding strategies between the operators are totally correlated. It is assumed that the EH operator is rational and avoids bidding high electricity prices, which results in a lower share of the market by the EH. This is because the electricity and heat market operators will resort to local production to supply the demands. From the market point of view, the EH can play two roles: a prosumer for the power market and a producer for the thermal market. From the standpoint of the power market, the EH can offer electricity to the market or consumers for economic purposes (for storing in an electric storage system or utilizing the EHP). The EH does not import heat from the thermal network. The reason is that arbitrage opportunities in the heat network are less competitive since the gap between peak and valley is not considerable compared with the electricity market where real-time prices, for example, can bring several arbitrage opportunities. Therefore, the understudy EH plays the role of heat supplier to the thermal distribution network. This is conducted by the EH thermal equipment such as EHP, CHP, and thermal storage system. In addition to all of these, the DSO benefits from the prediction module to forecast uncertain parameters and DR program, which bring technoeconomic merits to both systems. These two capabilities of the DSO are discussed further. Furthermore, we assumed that the EH and two other markets would reach certain agreements on the lower and upper bounds of offering prices to guarantee market fairness. Nevertheless, if congestion and other security constraints were considered, the EH might could posse market power in certain cases, which is not our aim in this article. In this structure, the focus is on the strategic bidding of the EH by leveraging the flexibility offered by electricity storage systems and thermal storage systems as well as crossarbitrage potentials among electricity, heating, and gas markets. The EH gets payment following a pay-as-bid agreement.
The power and heat networks are connected through the EH. Power flow status in the power distribution network is formulated by the linearized branch flow model. The EH submits its electricity (heat) offering/ bidding prices and quantities to the DSO (heat market operator), who subsequently clears the electricity (heating) market for minimizing total production cost. The EH then receives the energy contracts and gets payment according to the price bids. The bidding strategies in the two markets are spatially and temporally correlated owing to the energy conversion and storage capabilities.

| Objective function (OF)
The OF of the optimization problem is indicated in (1), which maximizes the EH profit. Part_A of the OF includes cleared quantity of heat and gas demand, and Part_B of the OF consists of the electricity generation and electricity demand. The EH submits quantities and prices of thermal and electric energies in the respective market and gets paid in accordance with the offering prices. Different energy prices in peak and valley hours precipitate arbitrage opportunities.

| EH constraints
Different from residential EHs, the considered EH here is able to sell electricity and heat to the power and heat networks, respectively, at its output side. In the EH, electricity can be used to charge the energy storage, or consumed by a heat pump that produces heat. Natural gas is burnt by a CHP unit to produce electricity and heat. Heat can be stored in a thermal storage system if necessary. The operational constraints of the understudy EH in the day-ahead energy markets in as (2)- (12): The architecture of energy hub (EH) and energy markets integration Lower and upper boundaries of other variables.
In the above formulations, Equations (2) and (3) show the EH electric and thermal energy balances, respectively. Equation (4) indicates the amount of electric demand that needs to be supplied. Equations (5) and (6) represent the state of charge (SOC) of electric and thermal energy storage systems, respectively. It is assumed that the initial and final SOC of energy storage systems is identical (SOC (0) = SOC (24)). Equations (7)-(8) denote binary variables of the electric and thermal energy storage systems. This is because simultaneous charging and discharging of energy storage systems are impossible. Equations (9)-(11) show lower and upper bounds of the main variables. For other variables like charging and discharging power of energy storages, the lower and upper boundaries are collected in (12).

| Power market clearing constraints
In this study, the power distribution network has a radial topology that could be modeled using a linearized branch flow model. 31 The complete formulation of the linearized branch flow method can be found in reference, 31 and here, we avoid repetition. The linearized branch flow method has advantages over DC branch flow, which is frequently used in transmission networks. For instance, reactive power consumption is considered in this model, while the DC branch flow model neglects this factor. The market-clearing model for power distribution system according to the linearized branch flow method is formulated as below 24 : In the above quadratic market clearing problem, the DSO clears the market minimizing the total cost of the system. In the first term of (13), the generation cost of the system is minimized. The second term is the paid cost to the EH, which needs to be minimized from the power market point of view. Equations (14) and (15) indicate inequality constraints and lower/upper bounds of the linearized branch flow model variables, respectively.

| Heat market clearing constraints
The heat network is conceptually different from the power network. It is a network of supply and return pipelines with identical topologies. Operationally, hot water, which is prepared by thermal resources, is injected into a specific node of the network and circulated from the supply pipelines to the return pipeline. Clearly, the returned water has a lower temperature, requiring energy to get higher degrees. The mathematical modeling for heat networks is proposed in Aslani et al., 32 consisting of hydraulic and thermal parts. Here, we focus on the heat market clearing problem, which is quadratically stated as below 24 : ( 18) In (16), the total cost of the heat network is minimized. The cost of heat work includes heat production of conventional boilers and the price paid to the EH. Equations (17) and (18) indicate the heat network thermal flow status.

| Karush-Kuhn-Tucker (KKT) conditions for market clearing problems
From the Equations (1)- (18), it can be observed that a bilevel optimization problem is formulated in which the EH is the leader and heat and electricity markets are the followers of the problem. Moreover, from Equations (13)- (15) and (16)- (18), one can conceive the convexity of the quadratic optimization problems. Therefore, by replacing KKT optimally conditions of the aforementioned problems, the problem stated in (1) transforms into a single-level optimization problem. In the below, KKT conditions are applied to market-clearing problems.  (13)-(15), the following formulation can be considered: In (19)- (20), the primal variables' feasibility constraints are denoted; in (21)- (22), the dual variables' feasibility constraints are indicated; and in (23), the complementarity condition is shown.

| KKT conditions for heat market clearing problem
Similarly, for the heat market clearing problem, the following KKT conditions can be implemented: In (24)-(25), the primal variables' feasibility constraints are denoted; in (26)- (27), the dual variables' feasibility constraints are indicated; and in (28), the complementarity condition is shown.

| Linearization method
The overview of the presented optimization problem shows mixed-integer nonlinear programming (MINLP). More specifically, the OF introduced in Equation (1), as well as the market-clearing problem in Section 3.4, are nonlinear and evidently hard to be solved and hence, need to be linearized in an appropriate way. This article applies a linear disjunctive method 33 to linearize slackness and complementary constraints. This method is generally used for linearizing the market-clearing problem. Furthermore, for bi-linear product terms like the terms shown in the OF α β ϑ ( , , ), the binary expansion method introduced in Pereira et al. 34 is implemented. The linear disjunctive and binary expansion methods are formulated in (29) and (30), respectively.
In Equation (29), M should have a large value which would make the transformation accurate. In (30), y and x are two continuous variables, and 2 K discrete points for approximating possible values of y in the interval y y [ , ]. low up Hence, the right-hand side of (30) is the linear expression of xy. More information about the method is provided in. 33,34 The linearized OF of the problem is denoted in the following formula:

| Mixed ML-based forecasting method
Accurate forecasting of uncertain parameters such as load demand data brings several economic advantages to the DSO, who is responsible for the operation of the distribution system. In contrast, inaccurate forecasting could waste energy resources and reduce the competitiveness of the markets. Therefore, precise prediction of uncertain parameters is an essential task for DSOs in the day-ahead operation of the electricity market. In this article, instead of using a single ML algorithm, a mixed approach is utilized. A combination of ML algorithms, as introduced by Faraji et al., 35 is used for forecasting uncertain parameters. In this approach, uncertain parameters of the system are collected at the initialization phase. Then, the collected parameters are forecasted by different ML algorithms such as artificial neural networks (ANNs) and adaptive neuro fuzzy inference system (ANFIS) in the computational phase. Day-ahead forecasted data by the prediction algorithms are extracted and evaluated by the actual data to achieve the bestforecasted pattern. For each type of uncertain parameter and forecasting algorithm, the pattern with the minimum Root mean square error (RMSE) is selected for day-ahead energy dispatch. The best pattern is calculated by comparing actual and predicted data as (37). This criterion shows that the most efficient ML algorithm is taken under consideration for each uncertain parameter so that the accuracy of forecasting improves greatly. In contrast to, 35 we considered RMSE in substitute of MSE, which is a more appropriate term than MSE for comparing prediction accuracy. Due to the hybridization nature of this approach, this model could be used to predict several uncertain parameters following learnable patterns. The structure of the introduced approach is illustrated in Figure 3.
This article utilizes three efficient ML algorithms such as multilayer perceptron (MLP)-ANN, radial basis function (RBF)-ANN, and fuzzy C-means clustering (FCM)-ANFIS to obtain the best-predicted patterns for the main goal, that is, optimizing the energy dispatch problem. The mentioned ML algorithms are adequately discussed and well-presented in previous studies, and therefore, we avoid repetition to save space. For more information, readers can refer to Carolin Mabel and colleagues. 36,37 Furthermore, the predicted solar irradiance data are used for calculating the output power of the PV system, formulated as below 32,38 :

| Price-based DR program
Generally, two popular DR programs are applied in the electricity markets: incentive-based DR programs and pricebased DR programs. 39 The purpose of both programs is to shave consumers' load profiles to have more flat load profiles, which is ideal for power system operators in general. In incentive-based DR programs, users are offered predefined incentives for their reduction in power consumption. One of the popular incentive-based DR programs is the time of use (TOU) program, in which the distribution network operator (DSO) defines two or more price values based on the time of day. 40 According to this pricing regime, peak periods have the highest prices and off-peak hours have the lowest prices. While in this study, a pricebased DR program is proposed which considers real-time pricing (RTP) for shaving consumers' load profile. Therefore, retail electricity prices in bus 2 of the network are considered for evaluation. The proposed formulation of the F I G U R E 3 The structure of the prediction method MOGHADAM ET AL.

| 2779
proposed DR program is provided in Equations (37)-(41). Figure 4 shows the steps required for the utilization of the DR program. It is worth mentioning that the DR is implemented on the active load profile. The proposed RTP DR is based on the forecasted values of load demand, so we can use the best-predicted load pattern for this purpose, as indicated in (37). Accordingly, real-time prices are obtained using the total and average load values as in (38) and (39). The DR program is then implemented in (40)

| RESULTS AND DISCUSSION
In this section, the results of the proposed optimization problem are presented in two different case studies. In the first case study, it is assumed that the electricity pricing is based on the TOU prices, while RTP and DR prices are taken into account in Case 2. Before presenting the results of the case studies, primary details of the study are presented in this section. As shown in Figure 5, a 32-node heat network and an IEEE 33-bus system are considered for distribution systems. The EH is linked to the power network through bus 2 and is linked to the heat network through Node 31. Moreover, two gas turbines connected to buses 7 and 13 supply the electricity needed by the DSO, and gas boilers in Nodes 1 and 32 provide the required heat of the network. Two PV plants (500 kW) are installed in buses 2 and 26 of the electricity network to support clean energy to the system. In addition, the EH independently purchases gas from the gas company to compensate for the reactive power, and to maintain the voltage profile of the network, two static var generators (SVGs) are equipped into buses 3 and 12. Essential parameters of the system are provided in Tables 2 and 3. Heat load profiles are presented in Figure 6.
Other important outputs of the optimization model are feeder capacity (3 MW) and the maximum gas inflow rate delivered to the EH (1.5 MW). The upper and lower bounds, as well as average heat offering prices ϑ ( ) t , are respectively 30, 12, and 25 $/MWh. In this study, 128 discrete points (K = 7) are considered in the binary expansion scheme. The offering electricity price β ( ) b must not be greater than the highest day-ahead electricity price multiplying (2.5/2), and the price of import electricity bid φ ( ) b must not be less than the retail price at bus 2 in the period t Δ . Furthermore, charging and discharging efficiencies of electricity and thermal energy storage are considered as 98%. The heating and electricity conversion rates of CHP and HP efficiency assumed 0.65, 0.35, and 3, respectively.
Moreover, the optimization model is solved using GAMS software and CPLEX solver while forecasting is performed in MATLAB software, and the predicted data are linked to GAMS. The specification of the laptop used is Corei5-CPU and 8GB RAM.

| Prediction results
One essential step in this article is to obtain the predicted values of uncertain parameters, that is, electricity load demand and PV output power, as the inputs of the dispatch model. In this subsection, prediction results are provided based on the proposed method. First, we tried to present the results in a more general way and then compare the day-ahead prediction results and select the best patterns of each type according to the presented method. Based on the literature, time-series forecasting is an appropriate option for having an accurate prediction of uncertain parameters. In the time-series method, a time delay (depending on the data) is typically applied, and prediction is conducted. We used 1-year hourly load demand and solar irradiance data (8760 h) for performing an accurate prediction. Hence, a time delay of 24 h is considered to perform the prediction. In addition, 70% of data is used for training, 15% for testing, and 15% for validating the data. First, we analyze the general prediction results over the designated horizon (Figures 7 and 8). Figure 7 shows the error bar graph for different ML methods. According to this figure, RBF-ANN and ANFIS algorithms demonstrated satisfying results, while it can be seen that the ANFIS model outperforms RBF-ANN. This can also be observed in Figure 8, where output and target values are depicted appropriately. More specifically, Figure 8 shows actual and predicted data, including training, testing, and validating data. Similarly, we can observe the prediction results of the solar irradiance, which can directly be used to compute the output power of PV systems. Figure 9 shows the error bar of the solar irradiance data and represents the superiority of the MLP-ANN and ANFIS models. In fact, the utilized algorithms demonstrated better performance with solar irradiance data. One reason is that the variability of solar irradiance data is much less than networks load demand data. As a result, it is much easier for ML algorithms to learn the patterns and cope with the nonlinearities of the data. Nevertheless, ANFIS is the most accurate model for predicting the solar irradiance data. The results of Figure 10 also confirm the accuracy of the ANFIS model. Previous claims can be verified by the information provided in Table 4. For instance, RMSE, standard deviation (STD), and R values for load demand and solar irradiance data show the superiority of the ANFIS model in comparison with other ML algorithms. While the prediction results show the general advantage of the ANFIS model, it is crucial to evaluate the best pattern for the understudy day.
As mentioned above, the need for finding the bestpredicted pattern for the understudy day drives us to consideration of an accurate statical criterion to distinguish between patterns. To do so, daily patterns are collected based on the predicted data, and the RMSE of each pattern is calculated in the presence of actual load and solar data. Figure 11 is provided for this regard and shows predicted and actual daily patterns to be used for finding the most accurate patterns of load demand and solar irradiance. As can be observed from Figure 11A, the major deviations occurred between 11:00-14:00 and 19:00-21:00. According to the suggested approach, RMSE of different ML algorithms such as MLP-ANN, RBF-ANN, and FCM-ANFIS is calculated as 0.0922, 0.0691, and 0.0651, respectively. Based on the obtained values by RMSE calculations, FCM-ANFIS will be selected as the best load demand pattern yielding the lowest RMSE compared with the other models. Moreover, the same procedure will be conducted on the selection of the best solar irradiance pattern among the predicted ones ( Figure 11B). RMSE values for this parameter achieved 0.0291, 0.0247, and 0.0242 for MLP-ANN, RBF-ANN, and FCM-ANFIS, yet showing the advantage of FCM-ANFIS over the other models. So, in this article, predicted patterns both by the ANFIS model are selected as the inputs of the dispatch problem.

| Implementation of the DR program
In this subsection, the proposed DR program is implemented on the network's load demand. This requires predicted load values from the mixed MLbased method to be imported to the program. As early discussed, the prediction results of the ANFIS model for load demand data were found more satisfying and hence, selected for implementing the DR program. TOU retail prices of the distribution network are presented in Figure 12A, consisting of three periods according to the time of day. Generally, the off-peak period includes 23:00-8:00, the mid-peak period includes 9:00-13:00, and the peak period includes 14:00-22:00. Figure 12B shows the two electric load profiles. The new electricity prices and load profile are achieved using the TOU and forecasted load profile (presented in Figure 12B). As shown in Figure 12A, electricity prices have real-time changes and correspond to the consumers' load pattern. The aim is to shave the load profile over the day, which is responded to in Figure 12B. It is essential to mention that the DR program is formulated for the active load demand of the electricity network. This is mainly because the flexibility of active load demand is considerably higher than reactive load demand in the residential sector.

| Optimization results
This section provides corresponding results of the presented optimization model, including EH dispatch and coordinated operation of distribution and heat networks. To understand the impact of uncertainty forecasting and DR programs on the outcome of the optimization model, two scenarios are considered in this study. The results of both scenarios are presented and discussed in two subsections. The suggested scenarios are as follows: (1) Scenario 1: Optimal dispatch of integrated energy grids considering uncertainty prediction. (2) Scenario 1: Optimal dispatch of integrated energy grids considering uncertainty prediction and DR program.

| Scenario 1
In this scenario, the proposed optimization model is solved using predicted load and PV values. Furthermore, the TOU electricity pricing system is utilized, as  Figure 12. According to the results, the profit of the operation of the EH is obtained as 566.5749 $. This profit is obtained from the revenues of selling electricity to the electricity distribution network (615.5576 $) and heat to the thermal distribution network minus the gas purchasing cost from the gas company (445.0667 $). Moreover, the operation costs of the power distribution network and heat distribution network are achieved at 4116.1 $ and 770.8931 $, respectively. Overall, the operation cost of the power distribution system is considerably higher than the thermal distribution network. To better understand the costs of each network, two pie graphs are provided in Figure 13, indicating the cost share of each system in the networks. As can be observed from Figure 13A, purchasing power from the upstream grid is the main cost of the grid (72% of total cost). The EH has the second-largest operation cost of the network (15%), and gas turbines 1 and 2 have lower share operation costs. The costs share related to the heat provision by EH and heat boilers are shown in Figure 13B. In contrast to the electricity distribution network, the highest operation cost in the thermal distribution network is dedicated to the EH, and gas boilers 1 (31%) and 2 (18%) are the other sources that impose additional costs to the system. Therefore, the role of the EH is significant in both networks and can have a substantial impact on the system costs.
In Figure 14, the charging and discharging of electric energy storage and thermal energy storage systems are indicated. From Figure 14A can be seen that the battery storage charges during off-peak periods (according to the TOU pricing) and discharges mainly on-peak periods. In this way, the battery is F I G U R E 8 Actual and predicted profiles of load demand data for (A) MLP; (B) RBF; (C) ANFIS optimized to minimize the operation cost of the EH. In addition, the thermal energy storage is discharged at 15:00 and charged at 24:00 to appropriately respond to the thermal load demand. The gas boilers in the structure of the system may not be able to supply all the load demand, and hence, the EH resources, including the thermal energy storage system, can play an important role.
In Figures 15 and 16, offering prices and quantities of electrical and heat energy are shown. Since HP benefits from a higher efficiency than the CHP system, the heat offering prices in intervals 1-5 are lower than the remaining time of the day during which the heat is provided by the CHP system. Heat quantity offers H ( ) b and electricity quantity offers and bids P P ( , ) db gb are equal to the cleared values. Through the cross-arbitrage among the gas market, heat market, and electricity market, the EH gains a profit of $566.0865. Therefore, the role of the EH is considerable in balancing and economic operation of the networks.

| Scenario 2
In this scenario, the proposed DR program is implemented on the optimization model, which means that the RTP pricing scheme is used instead of the traditional TOU pricing scheme. The profit of the system is significantly increased using the RTP pricing scheme. The profit of the system is achieved 967.8131 $, which shows a 40% increase in the profit of the system compared with the previous scenario. respectively. Similar to the previous case study, the operation cost shares of the power distribution system and thermal distribution system are illustrated in Figure 17. The share of each resource of the power distribution system is different in this scenario. For instance, the share of the upstream grid is decreased from 72% to 65%. This shows that due to the utilization of the DR program and RTP pricing, the load profile is improved, and hence less electricity is purchased from the grid. On the other hand, the EH share is increased in this scenario (8% increase). Similar to previous scenarios, two additional resources, such as gas turbines 1 and 2, have the lowest share of operation costs. The thermal distribution system, however, is faced with minimum changes compared with the previous scenario. The share of gas boiler 2 is increased from 17% to 18%, while the other resources have not changed. This is due to the fact that the DR program is implemented on the electric demand, so the thermal load demand did not undergo any changes. As a result, it is not expected to have considerable differences in the operation cost of the thermal distribution network. The results of electrical and thermal energy storage systems are shown in Figure 18. As shown in Figure 18A, the electricity storage system is charged in off-peak periods, then is discharged when the electricity price and consumption are peaked in the morning (8:00). The rest of the energy is charged in nigh peak hours, started at 15:00, and gradually discharges until the peak hours. Then the battery is charged at night's peak hour.
Considering the thermal energy storage ( Figure 18B), we observe that more charging has occurred in off-peak hours where the HP can charge the thermal storage and discharges in peak demand values.
Offering prices and quantities of heat and electricity are indicated in Figures 19 and 20. Similarly, the heat cleared prices and quantities follow the previous explanation in the previous scenario. In Figure 20, the electricity price and quantities of the network are illustrated, which infers more variations than the previous scenario. In fact, due to the real-time changes in prices, the cleared prices and quantities also involve changes. This is quite visible in the generation cleared prices.

| CONCLUSION
In this article, an energy system is proposed based on power and heat distribution networks which are linked through the utilization of the EH. Therefore, a bi-level optimization problem is formulated where the EH is the leader and electricity and heat networks are the followers. Using the KKT and suggested linearization method, the problem was converted to a single-level MILP and solved by the CPLEX algorithm. Another important consideration in this study was the application of a mixed ML (MLP-ANN, RBF-ANN, and ANFIS) approach in the forecasting of uncertain load demand and solar irradiance data. The method consisted of two stages. In the first stage, general forecasting is performed using different machine learning algorithms, and in the second stage, best-predicted patterns for the understudy day are extracted. The prediction results showed the superiority of the ANFIS model according to the defined evaluation criterion (lowest RMSE for load demand [0.0651] and solar irradiance [0.0242]). In addition, the predicted patterns are used for implementing a pricebased DR program. The DR program is then used in the optimization model. The scenarios are then suggested based on neglecting the DR program and considering the DR program. Optimization results show the great impact of the efficiency of the prediction method and FR program. The profit of the system by considering the DR program increased by 40% from 566.5749 $ to 967.8131 $, which yields the importance of utilizing the DR program. DR program led to a reduction in the electricity purchases from the upstream grid (from 72% to 65%). However, the heat network had no significant changes due to the changes in the power network. For future work, the authors will consider stochastic programming for considering uncertainties.

D ͠ Best
Best forecasted load pattern D ͠ Act Actual load data D ͠ fore Forecasted load data G t ( ) Solar irradiance | 2791