Multi‐stage coordinated operation of a multi‐energy microgrid with residential demand response under diverse uncertainties

Correspondence Yan Xu, School of Electrical and Electronics Engineering, Nanyang Technological University, 50 Nanyang Ave, Singapore, 639798. Email: eeyanxu@gmail.com Abstract This paper presents a two-stage stochastic operation method for a multi-energy microgrid (MEMG). The method can optimally schedule distributed generators, electric boilers, electrical chillers, and storage devices under the system technical constraints. Various uncertainties from renewable energy generation, electricity tariff, and load demands are handled by this method. At the same time, the indoor temperature is controlled to guarantee the thermal comfort of customers. At the first stage, the unit commitment of all generators and charging/discharging power of energy storage devices are obtained. At the second stage, real-time generation outputs from the generators are optimized given the uncertainty realization. The developed model is formulated as a mixed-integer linear programming (MILP) problem with the objective of minimizing the daily operating cost. Several case studies are carried out to examine the effectiveness and performance of the proposed method. Simulation results show that the proposed method can effectively save the operational cost and, at the same time, maintain the robustness of the MEMG and customers’ thermal comfort.

Usually, the unit commitment, which decides the on/off status of all the generators, is performed a day in advance. After the unit commitment scheme, the economic dispatch should then be decided for allocating energy generations to demands while satisfying system constraints [5]. The optimal MEMG operation has been intensively studied by researchers from the world over to improve the system energy utilization efficiency [6][7][8][9]. However, most of their studies are based on the deterministic optimization method, which implies that in their studies, the forecasted generation consumption of renewable energy sources and multi-energy demands are assumed sufficiently accurate.
Recently, the integration of renewable energy sources such as wind and solar power has led to the rapid development of MEMGs. Owing to the high volatility of renewable power sources, renewable generation outputs are generally difficult to predict with precision, posing significant challenges to microgrid operation [10]. Moreover, it has been indicated by some studies and reports [11,12] that the uncertainties from electricity tariffs and power loads also pose significant challenges for MEMG operation. Accordingly, an increasing number of researchers are focusing on considering the diverse uncertainties in the system operation. In ref. [13], a microgrid operation framework based on the model predictive control method is proposed. The prediction errors are considered in the method to handle uncertainties. A scenario-based stochastic optimal operation method is proposed to accommodate the uncertainties from renewable resources in microgrids in refs. [14,15], where all the scenarios are randomly generated and subsequently reduced to lower the computational burden. In ref. [16], a two-stage stochastic approach is used to handle the uncertainties in the system operation. In addition, risk management for power generation and trading is considered. Though the above existing studies consider the uncertainties from the renewable energy sources, they are only based on the traditional singleenergy microgrid. However, in the MEMG, the uncertainties arise from not only the power system but also the thermal one. Accordingly, the above studies are not sufficient for emerging MEMGs.
Moreover, as an important indicator in system operation, customers' thermal comfort is rarely investigated in the related research. Currently, the proportion of residential energy usage is approximately 30-40% of global energy consumption [17] in which the thermal consumption (heat or cooling) comprise a substantial portion. It is reported in ref. [18] that household loads have a significant effect on the seasonal and daily peak demands. In the study, a comfort-based game-theoretic approach that minimizes the system operating cost while keeping the users' comfort level is proposed in ref. [19]. In ref. [20], a home energy management algorithm for managing household loads is proposed based on the consumers' current comfort level. From the studies in refs. [19,20], which are confined into the deterministic study as in refs. [13][14][15][16], we can conclude the following: considering the thermal comfort of all the residents could contribute to an economic and friendly system operation; the customers' thermal comfort indicator is associated with the thermal demand response. Accordingly, considering the customers' thermal comfort in the system operation is similar to ensuring full usage of the thermal demand response scheme.
It should be noted that considering the multi-energy coordination is key for the improvement of the energy utilization efficiency. Further, to guarantee reliable and economic operation, the diverse uncertainties from both single-energy and multienergy systems should be considered. Finally, it is essential to consider customers' thermal comfort (demand response) in the MEMG operation, especially one based on residential buildings.
Given these insights, a two-stage (day-ahead stage and intraday stage) stochastic optimization approach for an MEMG with an objective of minimizing the daily system operational cost is proposed in this paper. In the proposed approach, multiple energy carriers such as heat, cooling, and electricity are integrated at the system level. The outputs from renewable generations, electricity tariffs, and load demands are modeled as stochastic scenarios using the scenario generation and reduction technique. A thermal comfort model is further incorporated to guarantee the customers' energy usage comfort level.
The contribution of this study compared with the existing studies is summarized as follows: i) A two-stage coordinated stochastic operation method is presented for an MEMG. Compared with current studies on the operation of single-energy microgrids [5,6], this study simultaneously coordinates the generation and consumption of the electricity, cooling, and thermal energy in the MEMG. ii) Residential thermal demand response is employed to guarantee the customers' thermal usage comfort so that overall economic benefits and dispatch flexibility can be markedly improved. The customer's thermal comfort level is an essentially important indicator for the multi-energy research and better utilization of this indicator could lead to more economic operations; however, this is not sufficiently resolved in other studies. iii) Diverse uncertainties from renewable energy sources, electricity tariffs, and multi-energy loads (including power, cooling, and heat loads) are fully considered to guarantee a reliable system operation. iv) Sensitivity analysis for the trade-off between the number of scenarios and solution time is performed. Accordingly, the effectiveness of the proposed method can be clearly justified.
The remainder of this paper is organized as follows: Section 2 presents a model of MEMG and mathematical formulations for the proposed system operation method. Section 3 deals with the problem of uncertainty modeling, scenario generation, and reduction. The case studies, simulation results, and analysis are presented in Section 4 while Section 5 draws the conclusion and provides a scope for future studies.

MATHEMATICAL MODEL OF MEMG
In this study, power sources in an MEMG can be divided into an uncontrollable generation (wind turbines [WTs] and photovoltaics [PVs]) and a controllable generation (energy storage devices, fuel cells, CCHP plants, and electric boilers). The schematic diagram of a typical MEMG is presented in Figure 1.
In this paper, we assume the studied MEMG is operated in gridconnected mode; thus, it can exchange power with the upperlevel grid while considering electricity tariffs, load requirements, and available generation capacities. In this section, the generation outputs of diverse components and their economic considerations, such as fuel cost, operation, and maintenance costs are modeled in detail. Practically, a CCHP plant comprises three units: a gas-fired micro-turbine, an absorption chiller, and a heat recovery unit [21]. CCHP plants can attain a high efficiency because the waste heat released from the electricity generation is absorbed and then utilized by the recovery unit/absorption chiller for heating or cooling. Owing to its capability of simultaneously generating electricity, heating, and cooling, the CCHP plant is a key component in an MEMG. The relationship between power and heat/cooling power of a CCHP plant at time t can be expressed as follows: where i and t denote the indexes of bus and time periods, respectively; H i,t CCHP ∕C i,t CCHP represents the heat/cooling power of the CCHP plant; r h ∕r c and COP h ∕COP c denote the heat recovery rate and the coefficient of performance of the heat recovery unit or absorption chiller, respectively; P i,t CCHP and CCHP denote the power output and electrical efficiency of the CCHP plant; L represents the heat loss ratio during the absorption process.

WTs and PVs
WT and PV are appliances that directly convert clean energy (wind's kinetic energy and solar energy) to produce electricity. Their generation outputs are related to the ambient wind speed and solar radiation, which can be formulated as follows: where P i,t WT ∕P i,t PV denotes the power outputs of the WT and PV panel, respectively; 1 , 2 , and 3 represent coefficients of the WT; v co ∕v ci denotes the cut-out/cut-in wind speed of the WT while v t represents the ambient WT at time t; P r,WT and V rate denote the rated power output and the rated wind speed of the WT; f pv denotes the derating factor of the PV panel; G T,t and G T ,STC represent the incident solar radiations at time t and under standard test conditions. P r,PV is the rated power of the PV panel.

Energy storage
Energy storage devices have been extensively installed in microgrids. This is because they can decouple the energy generation from demands, enhancing both the flexibility and technical performance of the entire system. In this paper, we consider three categories of energy storage: electrical energy storage (EES), heat storage tank (HST), and cooling storage (CS). It is assumed that these three types of energy storage devices share similar mathematical models [22]. This is because for energy storage, the relationship between the stored energy and charging/discharging power focuses on the microgrid operation problem.
Their models can be expressed as follows: where ES refers to energy storage (including EES, HST, and CS); E i,t ES denotes the stored energy; P i,t ES,ch ∕P i,t ES,dis are the charging/discharging powers of energy storage, respectively; is the decay rate; ch ∕ dis denote the charging/discharging efficiencies of energy storage while Δt represents the unit scheduling interval.

Electric boiler (EB) and electric chiller (EC)
EB and EC are two types of ancillary units in an MEMG that aim to convert surplus electrical energy to thermal or cooling energy, so that the dispatching flexibility can be substantially enhanced. According to ref. [21], the mathematical models of EB and EC are similar, which can be formulated as follows: where H i,t EB ∕C i,t EC denotes the heat/cooling outputs; represents the electric power consumed; EB ∕ EB denotes the transferring efficiency of /by the EB and EC, respectively.

Buildings (heat/cooling load)
Recently, the advancement of information and communication technologies has provided an excellent opportunity to develop a residential demand response for future smart grids. The residential demand response is a type of demand response program which covers customers' responsiveness and their behavioral issues. Thermal energy (heating/cooling) for buildings usually satisfies up to 90% of the total thermal demand during winter and summer [23]. According to previous studies, thermal loads are usually fixed [21]; however, given thermal inertia, residents can be comfortable with constant indoor room temperature. Consequently, thermal loads that are highly dependent on the temperature can be relaxed to enhance the operational flexibility. Accordingly, the specific modeling of the relationship between the indoor temperature of buildings and thermal loads is modeled to capture the thermal dynamics. Therefore, thermal comfort can be utilized as thermal demand response.
In winter, it is assumed that the indoor temperature in a building at time t is related to the temperature in the previous time period and the heat energy obtained from an exchanger (as illustrated in Figure 2). In addition, the heat exchanger connected to bus b is responsible for all the buildings near bus b. According to ref. [24], the indoor temperature of building m can be expressed where air denotes the specific heat capacity of the air; T m,t in denotes the indoor temperature of the building; T t am represents the ambient temperature; H m,t in denotes the heat energy obtained from the heat network; R T represents the heat resistance of the building shells (including windows, doors, and roofs); H m,t load denotes the design heat load for building m at time t; T t in,min represents the minimum indoor temperature that is required to keep customers feeling comfortable during winter; H i,t load represents the design heat load at bus i; H i,t ex denotes the heat power of the heat exchanger at bus i; M denotes the number of buildings near bus i.
Similarly, for cooling energy in summer, the indoor temperature of building m can be expressed as follows: where C m,t in represents the cooling energy obtained from the MEMG; C m,t load denotes the design cooling load for building m at time t; T t in,max represents the minimum indoor temperature that is required to keep customers feeling comfortable during summer; C i,t load denotes the design heat load at bus i; C i,t ex represents the heat power of the heat exchanger at bus i.

2.2
Optimal coordinated scheduling

Objective function
In this problem, the objective is to minimize the daily operational cost, which includes the fuel costs of the CCHP plants, operation & maintenance (O&M) cost, turn-on/turn-off cost, as well as the cost of power exchange with the main grid [22].
where F gas and L HVNG denote the unit price of natural gas and the lower heating value of natural gas, respectively; g refers to all the controllable components in MEMG (namely the CCHP plant, WT, PV panel, energy storage, EB, and EC); F om g denotes the unit O&M cost; F SU g ∕F SD g denotes the unit turn-on/turnoff cost of the component, respectively; O i,t g represents a binary variable representing the on/off status of component (O i,t g = 1 implies the on status while O i,t g = 0 implies the off status); P t buy ∕P t sell denotes the purchasing/selling power from/to the main grid; F t tr denotes the electricity tariff.

Constraints
In this paper, the constraints are as follows: the power outputs limits and ramp up/down limit of components (20) and (21); maximum charging/discharging power limits and energy limits, starting/ending state of charge limit and state limit of energy storage (22)-(26); power exchange limit with the upper grid (27)-(29); linearized distribution network power flow constraints (30)-(35); energy balance limits (36)-(38) and indoor temperature constraints (39). Note that (39) implies the indoor temperature can be flexible in a range, which will not affect the customers' comfort. The following can be the interpretation of the residential demand response scheme: buy ∕O t sell are binary variables with O t buy = 1 indicating that the MEMG is the purchasing power from the upper gird; b denotes the index of the branches and B i ∕B i,i+1 denotes the sets of branches that are connected to bus i or between buses i and i + 1; P i,t ∕Q i,t represents the active/reactive power flowing on branch b; P b,l,t ∕Q b,l,t denotes the active/reactive power flowing on the lateral branch of branch b; P i,t load ∕Q i,t load represents the reactive/active power load at bus i; V i,t denotes the voltage at bus i ; r b ∕x b represents the resistance/reactance of branch b; ΔV max represents the maximum voltage deviation; P max ∕Q max denotes the maximum active/reactive powers flowing on a branch; P i,t EES,ch ∕P i,t HST,ch ∕P i,t IST,ch denotes charging powers of EES/HST/CS, respectively, while P i,t EES,dis ∕P i,t HST,dis ∕P i,t IST,dis represents discharging powers of EES/HST/CS, respectively; H i,t load ∕C i,t load are heat/cooling loads at bus i.

SOLUTION METHOD
The renewable energy generation outputs, electricity tariffs, and power/heat/cooling demands are the uncertainties to be considered in this paper [25]. In previous studies, stochastic and robust optimization methods are utilized [26] to handle uncertainties. The robust optimization method handles optimization problems in which a certain measure of robustness is sought against any uncertainty realization. Thus, the goal is to find a solution that is feasible for the worst uncertainty case. Though reliable enough, the worst-case-oriented robust optimization tends to be too conservative. This is because the possibility of the worst-case occurrence would be extremely low or never happen in the real world [27]. Using this low-possibility worst case for the system operation would not be prudent. Conversely, the stochastic approach sets the expectation of all scenarios as the objective, which would be more reasonable for providing operation references compared with using the worst case. Furthermore, with the development of more advanced forecasting techniques, the probability distribution for diverse uncertainties can be well described [28]. Accordingly, the stochastic optimization method is applied in this study.

Uncertainty modeling
The uncertainties of forecasted wind and solar power generations can be expressed as follows: where f P i,t WT and f P i,t PV represent the probability density function of beta distribution for forecasted power P i,t wT and P i,t PV , respectively; denotes the normalized factor; and represent two parameters that form the beta distribution.
As for electricity tariffs and power demands that are highly related to customers' behaviors, their randomness and forecasting errors can be modeled using the normal distribution [29]. For example, for a certain predicted P i,t load , its probability density function f P i,t load can be formulated as follows: where and are the mean and standard deviation for normal distribution, respectively. Note that for heat and cooling demands, it is assumed that the variation of the outdoor temperature follows the normal distribution, and load demands can be calculated using (9) and (12).
It should be noted the parameters for beta and normal distributions should be carefully selected according to the requirement for different cases.

Scenario generation
In this paper, we use the Latin hypercube sampling method [30] to generate a high number of random scenarios for the stochastic variations of electricity tariffs, renewable energy source generations, and load demands. Assume S scenarios are randomly generated to represent uncertainties, the cumulative distribution function is then divided into S intervals, each with an equal probability of 1∕S . A value is then randomly selected from each specific interval to obtain S scenarios.

Scenario reduction
In stochastic programming, the number of scenarios is directly proportional to the computational burden; hence, an efficient scenario reduction technique is implemented to reduce computational time. In this paper, we utilize the simultaneous backward reduction technique [29,31] to choose a certain number (R) of scenarios, which comprises the following steps: Step 1: Let S denote the set of all the initial scenarios, and S d represent the sets of scenarios that need to be deleted (initially null). s (s = 1, 2, … , S ) is the index of the initial In this first step, we compute the distances between all the scenario pairs D s,s ′ = D( s , s ′ ), (s, s ′ = 1, … , S ); s denotes the average/mean value of all the scenarios in the full scenario set S.
Step 2: For a specific scenario n, D n (r ) = minD n,s ′ , n, s ′ ∈ S and n ≠ s ′ . r denotes the nearest scenario with scenario n.
Step 3: Compute the probability D n (r ) = s * D n (r ), n ∈ S . Select d, which satisfies D d = min D n , n ∈ S .
Repeat the above steps until the required selected number R is satisfied.

Two-stage stochastic optimization method
In this paper, we propose a two-stage stochastic optimization method for the optimal scheduling problem of MEMG while  considering diverse types of uncertainties. Figure 3 shows the framework of the proposed method. In the first stage, a "here-and-now" decision x (unit commitment decisions of the controllable generators) is made before the uncertainties are realized, whereas a "wait-and-see" decision y (power contributions of the controllable generators, energy storage; power exchange with the grid and customer-side indoor temperature) is optimized in the second stage intra-day when the uncertainties are realized. Term y acts as a compensation for the difference between the predicted value in the first stage and the real value in the second stage. The first stage objective is to minimize the operating cost one day in advance while considering future uncertainties. The second stage aims to minimize the operational cost including the fuel cost of the CCHP, the O&M cost of the controllable energy units, and the cost of power exchange with the main grid when the uncertainties are realized. The two stages are coordinated through constraints and parameters. The "here-and-now" decision vector x is viewed as a state variable in the second stage.
In the first stage, the uncertainty models of electricity tariffs, load demands, and renewable energy source generations [(40)-(42)] are used to generate 2000 forecasting scenarios. Subsequently, 10 representative scenarios are selected. Using the 10 representative scenarios, the first stage problem can be formulated as a mixed-integer multi-objective linear optimization model and expressed as follows: F (x) denotes the cost of turning on/off the controllable units while ∑ n∈S n L(y n ) represents the expected operating cost in the second stage, considering uncertainties. The firststage decision variable x, namely the on/off states of controllable units (CCHP plants, EBs, and ECs) for each hour of the  The above objective function (44) can be rewritten in detail as follows: Afterward, given the first-stage commitment decisions, generation outputs of CCHP plants, EBs/ECs, and charging/ discharging powers of energy storage (including EES, HST, and CS) are optimized. The simplified objective function in the second stage is shown as follows: Here, the O&M costs of WTs and PV panels outputs are already fixed because their uncertainties are realized in the first stage. Thus, in this stage, the objective is to minimize the summation of the O&M costs of CCHP plants, EBs/ECs, energy storage, fuel costs of CCHP plants, and the power exchange cost with the upper grid.
Moreover, to employ a linear programming approach, the non-linear cost functions (17) and (18) are linearized for higher computational efficiency.

Test system
Based on the IEEE 33-bus distribution radial system, we constructed an MEMG for the case study. The schematic diagram of the tested MEMG is illustrated in Figure 4. In this section, we chose a typical day during summer and winter, respectively, in Northern China for the case study. Three types of load requirements were considered at each bus: electrical, heat, and cooling loads. Considering long-distance heat/cooling transferring results in substantial energy waste, each heat or cooling load was assumed to be served by a nearby heat/cooling source [22]. Additionally, it was assumed that on the summer day, customers required only the cooling load while on the winter day, only the heat load needs to be satisfied. The system parameters are listed in Table 1. Table 2 presents the parameters of generators and auxiliary units. The size and cost of items of energy storage are presented in Table 3 [21]. Note that if components installed at different buses are of the same type, then they are of the same size. Other detailed specifications of WTs and PV panels can be found in refs. [22,32].

First-stage simulation results
In the first stage, 2000 uncertain scenarios are initially generated based on historical values. Next, 10 scenarios are selected from the initial 2000 scenarios using the simultaneous backward reduction technique. The power output/input of the CCHP plant and EB/EC under 10 representative scenarios are shown in Figures 5 and 6. The on/off statuses of CCHP plants and EB/EC are listed in Table 4. Figures 7 and 8 depict the charging/discharging power of energy storage (the power above zero implies energy storage is charging, whereas power below zero implies energy storage is discharging). The total expected daily costs for winter and summer cases are ¥6457.081 and ¥5684.293, respectively. By the CPLEX solver, the first-stage optimization is solved within 4.2 min, which is efficient enough for the first-stage (dayahead) operation.
It can be noted from Figure 7b that the heat energy storage tank receives surplus heat energy generated by the CCHP plant during the periods 01:00-3:00, 16:00-17:00, and 23:00-24:00; releases the stored heat energy during the periods 15:00-16:00 and 19:00-22:00 when the CCHP plant and EB are unable to generate sufficient heat energy to satisfy the increasing heat load demand. It should be noted that the above energy storing/releasing process is repeated a few times during the day to decouple the electrical and heat generation of the CCHP plant.
The peak cooling load demands usually appear at approximately noon (10:00-14:00). During this period, CS dis-charges cooling power to satisfy the customers' requirements ( Figure 8b).

Second-stage simulation results
The second-stage operation is tested when the online renewable energy source outputs, customers' load demands, and electricity tariffs become available every hour. The real-time data inputs in this stage are plotted in Figure 9.
With the above input data, the simulation results for 24 h of power scheduling in the second stage are illustrated in Figure 10. The indoor temperature is demonstrated in Figure 11. Additionally, the daily costs in the second stage for winter and summer are ¥6419.249 and ¥5714.788, respectively. The computation time for the second stage is 2.7 s for the winter case and 2.4 s for the summer case. The solution time indicates that the proposed stochastic optimization is efficient and effective for the second-stage operation.
In the case study, the heat/cooling load demands are all sufficiently satisfied by the CCHP plant, EB/EC, and HST/CS, whereas electrical loads are served by the CCHP plant, WT, PV, EES, and power exchange with the main grid. It can be noted from Figure 10 that multiple energy sources/carriers are coordinated in this approach. Owing to insufficient solar irradiance during the period 19:00-21:00 (Figure 10a), PV outputs are equal to zero and WT outputs are reduced because of low wind speed. To compensate for the energy deficiency, the CCHP plant increases its electrical power output. Consequently, the heat power generation of the CCHP also increases. Further, the heat load demand concurrently decreases slightly as illustrated in Figure 9a while the surplus heat energy generated by the CCHP plant is stored in the HST (Figure 7b).
In winter (Figure 10b), the peak cooling demand appears during the period 12:00-14:00; therefore, the CCHP plant generates an increased amount of cooling power and the CS starts to  In summary, all the units are optimally scheduled under diverse uncertainties while the indoor temperature is kept within the customers' comfort range ( Figure 11).

Comparison and sensitivity analysis
To fully illustrate the benefits of our proposed two-stage stochastic operation approach, two other typical operation approaches are compared as follows: Approach 1: Single-stage deterministic operation. In this approach, no uncertainties are considered (renewable generation outputs, electricity tariffs, and load demands are deemed as accurate). The problem is optimized at a single stage. Approach 2: Single-stage stochastic operation. The uncertainties are modeled and considered, but the components are scheduled in a single stage with one representative scenario selected. Approach 3: Two-stage stochastic operation. Our proposed approach.
The above three approaches are all tested on the same 33-bus MEMG (Figure 4), and thereafter simulated using 2000 scenarios generated using the Latin hypercube sampling method. The average total cost and voltage satisfaction rate under 2000 scenarios of the above-mentioned approaches are listed in Table 5.
Simulation results show that the two-stage stochastic optimization method is efficient in solving MEMG's operation problem while considering the diverse randomness from RES, loads, and electricity price. It can be shown that although the average cost is similar among the three approaches, our approach can significantly increase the robustness of the entire system. The proposed approach can achieve the highest voltage satisfaction rate. The improvement of voltage quality is mainly because of the detailed consideration of diverse uncertainties from renewable energy generation, electricity tariffs, and customers' load demands in the first-stage day-ahead scheduling in the proposed approach. The uncertainty realization is further compensated by the second-stage online dispatch. Furthermore, a sensitivity analysis is conducted to select an appropriate number of scenarios R to ensure optimal operational performance and computational efficiency. Different R values are analyzed to show the sensitivity of the operational cost and computational time over R, which is illustrated in Figure 12. Note that the operational cost refers to the daily cost in the second stage, and the computational time includes the runtime in both stages.
According to Figure 12, as the scenario number R decreases, the computational time decreases accordingly while the operational cost increases. It can be shown that a significant portion of R results are in the lower operation cost mainly because of not only a better representation of uncertainties but also increased computational burden. Therefore, a trade-off between the computational efficiency and operational performance has to be made when it comes to the selection of scenario number R.

CONCLUSION
A two-stage stochastic operation approach for MEMG that considers diverse uncertainties is proposed in this paper, to optimally coordinate diverse types of energy sources/demands. In the first stage, the commitment decisions of the CCHP plant, EB/EC, are optimized to minimize the operation cost.
In the second stage, the real-time power output of the CCHP plant, EB/EC, the charging/discharging power of energy storage devices and customer-end thermal demand response, power exchange with the main grid, and energy storage are optimized to compensate for the first-stage decisions when diverse uncertainties are realized. The proposed method is compared with two other methods. The simulation results validate that the proposed method can achieve a comparable and reasonable economic performance and significantly higher robustness in terms of voltage satisfaction rate while satisfying customers' thermal comfort.
For future studies, first, the dynamic system operation and modeling can be studied as this paper covers only the steadystate operation. Moreover, the detailed modeling of customer's multi-energy satisfaction degrees can be further involved to formulate a multi-objective optimization problem. Finally, more advanced uncertainty consideration methods that could combine merits from both stochastic and robust optimization methods can be further utilized.