A Day-ahead and Day-in Decision Model Considering the Uncertainty of Multiple Kinds of Demand Response

The uncertainty of demand response (DR) will affect the economics of power grid dispatch due to the randomness of participating users’ intentions. According to the different working mechanisms of price-based demand response (PBDR) and incentive-based demand response (IBDR), the uncertainty models of two types of DR were established, respectively. Firstly, the fuzzy variable was used to describe the load change in PBDR, and the robust optimization theory was used to establish the uncertain set of the actual interruption of the interruptible load (IL). Secondly, according to the different acting speed of the two types of DR, they were deployed in the two-stage scheduling model with other output resources; then based on the fuzzy chance constrained programming theory and multi-stage robust optimization theory, the dispatch problem was transformed and solved by the bat algorithm (BA) and the entropy weighting method. Consequently, intraday running costs decrease with increasing confidence of day-ahead, but increase with day-in reliability, and the economics of the model were verified in the improved IEEE33 node system.


Introduction
In 2010, the six ministries and commissions, including the National Development and Reform Commission and the Electricity Regulatory Commission, jointly issued The Measures for Power Demand Side Management [1]. Based on the user electricity information collected and analyzed by grid companies [2], many provinces gradually implemented time-of-use pricing. In addition, thanks to the update and development of self-control technology and real-time information collection technology, power equipment control devices can realize self-manipulation to adjust the demand. The vigorous development of the demand side management of China's power market is conducive to suppressing the fluctuation of output caused by the large amount of renewable energy connected to the power grid [3].
Although the grid-related institutions analyze and classify the user's electricity habits in advance, the time-of-use pricing is fundamentally based on users' voluntary participation and changing their own electricity usage behavior [4], which results in unpredictable demand response (DR). There are four types of incentive-based demand response (IBDR) [5]: direct load control (DLC), interruptible load (IL), demand side bid (DSB) and emergency demand response (EDR), but in the gradually opening electricity market in China, the latter two are relatively rare. The IL is different from the DLC. The way it works is that the user has to turn off the power after receiving the notification. Due to the randomness of the user's wishes, it is possible that some users will not strictly obey the agreement. This is the root cause of the uncertainty of the IL [6].
In [9], Monte Carlo method is used to simulate the charging and discharging behavior of EV. In [10], demand response and electric vehicles are aggregated into virtual power plants to participate in market scheduling, and multi-scenario method is used to simulate the uncertainty of market electricity prices. The multi-scene method expresses various possibilities of uncertainty with a huge amount of scenes, which can be simplified by Latin Hypercube Sampling (LHS) for solving [11]. Reference [12] has established scenarios for wind, photovoltaic and load prediction errors, using LHS to reduce the scene.
Although the multi-scene method is simple, this method requires a lot of data and has strict requirements for the solving ability. Therefore, many scholars hope to deal with the uncertainty in the power system by stochastic optimization [13][14][15]. In [13], considering the effects of uncertainty such as outdoor temperature, basic power load, occupancy rate and unbalanced price, a two-stage stochastic programming method is adopted. Reference [14] obtains the probability function of the virtual plant deviation based on IBDR and price-based demand response (PBDR) from historical data. In [15], the probability function is used to describe price-based demand response, and the power flow is solved by the point estimation method.
The stochastic method reduces the solution pressure while maintaining the uncertainty complexity, but the probability distribution of uncertainty is difficult to obtain [16]. Robust optimization takes the worst case into account in the solution process and derives decision results based on this, thus there is no need to consider the probability function [17]. In [16], the scheduling is divided into pre-scheduling phase and re-scheduling phase. In each phase, the worst case of uncertainty is considered. The combination of column-and-constraint generation (C&CG) and Benders decomposition is used to solve the optimization problem. In other studies [18,19], the robust interval variable is used to describe the uncertain part of the IL, and it is transformed into a deterministic problem by using the robust peer-to-peer conversion method. In [20], wind turbine output and price-elastic demand curve work jointly, and total social welfare is maximized in the worst case considering wind turbine output and demand response.
There are two common limitations in the above literatures: (1) the uncertainty of DR is modeled by one single method, but there is not only one kind of DR. The accuracy of one single model is insufficient and (2) when considering the scheduling with DR, the scheduling plan only involves single time scale, and the speed advantage of different types of DR cannot be fully applied.
The main contributions of this paper are as follows: • According to the different mechanisms of price-based demand response and incentive-based demand response, the fuzzy function and robust interval variable are used to describe the uncertainty of time-of-use model and the IL.

•
Considering the response speed of different types of resources, and coordinating with the output resources such as wind turbine and gas turbine, the day-ahead and day-in interactive output decision model is established.

•
The fuzzy stochastic optimization method and the multi-stage robust optimization method are used to deal with the uncertainty of different demand responses respectively, and the interactive decision model is transformed into a deterministic model that can be processed by traditional algorithms.
The rest of the paper is organized as follows: in Section 2, mathematical models of different types of demand response uncertainty are established; in Section 3, a two-stage mathematical model that considers demand response, micro gas turbine, wind turbine and other output resources is established. In Section 4, the model is transformed into a deterministic model by using fuzzy chance constrained programming theory and two-stage robust optimization theory. Section 5 presents the simulation results of the example. Section 6 draws several conclusions.

Uncertainty Models of Two Types of DR
Demand response can be divided into price-based demand response and incentive-based demand response. This paper mainly considers the time-of-use model of PBDR. The common types of incentive-based demand response include direct load control (DLC) and the IL. Considering the difference in mechanism of action, this paper analyzes and models the uncertainty of the above demand response.

Uncertain Model of Price-Based Demand Response
The demand side management with the adjustment of electricity price as the main means, according to the principle of demand elasticity in economics, changes the electricity price to affect the user's consumption behavior, that is, the electricity usage. The electricity price elasticity coefficient is generally used to predict the user's response at a certain price level, but this method is interfered by many factors, resulting in the uncertainty of the time-of-use model.
The uncertainty of the price-based demand response depends on two factors mainly, the uncertainty of the baseline load and the uncertainty of the electricity price elasticity coefficient.
Baseline load forecasts often take into account many factors, including but not limited to population, industrial development status, and policy factors. These factors are often difficult to quantify and have unpredictable effects, so the prediction of the load is inevitably deviated from the actual load.
The core of time-of-use is to use the user's sensitivity to price to change the power consumption behavior. The essence is that users voluntarily participate and adjust their own electricity usage. But the behavior of users is difficult to predict and quantify, which is the root cause of the uncertainty of demand response.
The introduction of the electricity price elasticity matrix quantifies the impact of the baseline load and the electricity price elasticity coefficient on the demand change after PBDR. Compared with other PBDR models, it can predict demand changes more accurately. In general, the demand after PBDR is as follows: where E TOU = E TOU.f , E TOU.p , E TOU.g T represents the column vector of demand after PBDR; E TOU.f , E TOU.p , and E TOU.g represent the demand of the peak, flat and valley periods after PBDR; E 0 = E 0.f , E 0.p , E 0.g T is the column vector of the demand in each period before PBDR; E 0.f , E 0.p and E 0.g represent the initial demand of the peak, flat and valley periods; c f and ∆c f represent the original electricity price and electricity price change at the flat period; c p and ∆c p , c g and ∆c g are similar. M is the electricity price elasticity matrix, and its expression is as follows: where x,y = p,f,g, m xy is the electricity price elasticity coefficient. It shows the ratio of the rate of change of electricity price in y period to the rate of change of demand in x period. ∆E p , ∆E f , ∆E g represent the change of demand at the peak, flat and valley periods; when considering baseline demand uncertainty and user response uncertainty, the following relationship exists: where, E 0.x and m xy represent the baseline demand and elastic coefficient considering uncertainty. The membership function is an important means to deal with fuzzy parameters. In various membership functions, the trapezoidal membership function is better than the triangular membership function [21], which is closer to the actual situation. Therefore, this paper uses the trapezoidal membership function to deal with the uncertainty of price-based demand response. The membership function of the baseline demand and elastic coefficient are as follows: where m xy1 , m xy2, m xy3 , m xy4 , E x1 , E x2 , E x3 , E x4 , are membership function parameters. Its function graph is shown in Appendix B. When m xy varies in [m xy2 , m xy3 ], the function value is 1, and when m xy deviates from the expected interval, its function value decreases as the deviation increases. It can be seen that the function value indicates the extent to which the evaluated target deviates from the expected range.
The price elasticity coefficient reflects the users' sensitivity to changes of electricity prices. When the price change is small, the users are not sensitive to price changes, with low participating motivation, so the response is highly uncertain. When the electricity price is different from the usual electricity price greatly, the electricity price has a greater impact on the user's electricity consumption behavior. The participating motivation is greatly increased, so the response uncertainty is small. Therefore, k x and k y are introduced to indicate the relationship between the price change and the membership function parameters. Compared to constant parameters, membership function parameters that vary with electricity price change can represent different uncertainties. The membership parameters have the following relationship with the electricity price: where, k x and k y represent the growth rate of uncertainty. When the electricity price changes greatly, m xy1 tends to m xy2 , the uncertainty is small, the electricity price changes little, m xy1 is far away from m xy2 , and the uncertainty becomes bigger.

Uncertain Model of Incentive-Based Demand Response
Most of the current research considers the incentive-based demand response as a special output unit. As the output unit, the incentive-based demand response has no limitation on the climbing rate, no starting cost, and the response is fast. It is an ideal resource for dispatching.
The incentive-based demand response is divided into the IL and DLC. Generally, the interrupt capacity, interruption time, compensation price and maximum response capacity and duration in the whole scheduling period are specified in the form of contracts. DLC is controlled by the power grid company. It can be directly interrupted by the dispatching center in a very short time when it is called. The uncertainty is mainly caused by communication lag and equipment failure. It is not discussed in this paper.
Although the IL has a contractual agreement, the essence is that the users spontaneously participate in the response after receiving the notification from the dispatching center. Therefore, there is a response uncertainty caused by the randomness of the user's willingness to participate. This problem is the scope of this paper.
After the IL users receive the response notice, due to the high cost of breaking the contract, even if the users have the willingness to break the contract, the users will generally comply with the agreement of the interruption duration and partially default on the interrupted capacity. The actual interruptible power that user j can respond to at period t is expressed as follows: where,P t IL. j indicates the actual interruptible power of the user j at time t, P t IL. j indicates the planned interrupted power of the user at that time, and P t IL. j indicates the indeterminate portion of the interrupted power of the user. p t IL. j is the upper limit of the absolute value of the deviation between actual interrupt amount and planned interrupt amount in the historical data. τ j is the uncertainty coefficient, indicating the degree of uncertainty of the response of user j.
In addition, in the case that residential consumers are not eligible to directly participate in DR market [22], load aggregators (LA) can aggregate small and medium-sized loads such as residential load and commercial load, and participate in market bidding competition on their behalf [23]. The power grid dispatch center signs contracts with LA, which can be equivalent to users with large IL resources, from the perspective of the grid.
In order to analyze the influence of the uncertainty of the IL on the scheduling in a system, we define Γ as the robustness factor of the system. The interrupted power of the systemP t IL. j can be expressed as: where, M represents the number of users participating in IL; e t j represents the state of user j during period t which is 1 when the user is called. Changing the robustness factor Γ can change the uncertainty of the IL in the system. The robustness factor determines the uncertainty of IL, which can be set by the dispatch centre and reflects the dispatch centre's tendency toward risk. When the robustness factor is 0, each user strictly abides by the contract; when the robustness factor increases, the actual interrupted amount may seriously deviate from the planned value.

Two-stage Interactive Decision Model Considering Uncertainty of Demand Response
In the model of this paper, the operating cost is considered as the objective function in the day-ahead and day-in stage. In addition, in order to avoid the situation that the electricity price is too high or too low, the user transfer coefficient is also taken as one of the objective functions of the day-ahead stage. The constraint part includes the upper and lower limits and the climbing constraints of the conventional unit. The cost in day-ahead stage include the operating costs of wind turbine and photovoltaic unit, the operating costs of gas turbines, the cost of purchasing electricity from the main network, and the negative value of revenue from the sales to the microgrid. The expression is as follows: (c t WT P t WT + c t pv P t pv + c t mt P t mt + c t buy P t buy − c t pri P t MG )∆T (10) where, T represents the number of operating cycles; c t WT and c t pv represent the unit cost of wind turbines and photovoltaic units in the period of t; P t WT and P t pv represent the output of wind turbine generation and photovoltaic units in the period. c t mt represents the unit cost of the gas turbine in the period, and P t mt represents the output of the gas turbine in the period. c t buy indicates the unit price of electricity purchased from the main network in the period, and P t buy indicates the power purchased from the main network in the period. c t pri represents the unit price of selling electricity to the microgrid in the period, and P t MG represents the power purchased by the microgrid in the period, ∆T is the duration of the scheduling period.

User transfer coefficient
The main purpose of the time-of-use mechanism is to reduce peak-to-valley difference of load curve, alleviating the pressure of units. However, excessively high or too low electricity prices in some periods will bring inconvenience to the normal life of ordinary users, and reduce the enthusiasm of users to participate in DR. Therefore, the user transfer factor is included in the objective functions in the day-ahead stage: where, S is the user transfer coefficient, Pload t 0 is the load in the period of t before the implementation of the time-of-use price, Plaod t is the load in the period after the implementation of the time-of-use.

Restrictions
In the period of t, there is the following relationship: where, P t L is the fuzzy parameter of the load after the implementation of the time-of-use pricing. According to the fuzzy chance constrained programming theory, the decision made should satisfy the constraint at a pre-set confidence level.
In order to avoid the system's power shortage, Equation (12) can be converted into the following formula: where, α indicates the confidence level that the constraint needs to meet. The upper and lower limits of the output of gas turbine and the climbing constraint are as follows: ∆T · γ down < P t mt − P t−1 mt < ∆T · γ up (15) where, P mt.max and P mt.min represent the upper and lower limits of the unit's output, γ up and γ down represent the upper limit of the unit's output ramp rate, and ∆T should T be italic like before? represents the duration of a scheduling period.

Objective Function
Taking into account the IL's advantages of speed, it is included in the intraday scheduling resources. The objective function is the operating cost during the intraday scheduling, with constraints of the power balancing, the IL users and so on.
To minimize the running cost during the intraday scheduling period, the objective function is as follows: where P t IL represents the planned power of the IL in the period of t, c t IL represents the unit price of compensation for the IL in the period; P t R represents the output of peak load regulation (PLR) resources called in the period, c t R represents the unit price at which the PLR resources is called, and ∆P t buy is the difference between the actual and the planned power of purchased electricity from main network, and should be a positive value, c t def is the unit price of compensation to the main network.

Restrictions
The system's the power balancing constraints is as follows: whereP t IL represents the actual interrupted power of the system in the period t. Because the load is a fuzzy parameter, the equation constraint that satisfies the confidence β is: User j who participates in an incentive-based demand response has the following constraints: where, P min.j and P max.j respectively represent the lower and upper limits of the response; D j is the upper limit of the total response time. The upper and lower limits of the output of the conventional unit and the climbing constraints are the same as Equations (14) and (15).

Transformation of Uncertainty Model of PBDR
In the day-ahead stage, only the uncertainty of the load after PBDR needs to be considered, which is composed of the fuzzy variable of baseline load and the fuzzy variable of electricity price elastic coefficients. In the day-ahead model, the most important thing is the conversion of the opportunity constraint.
Taking the load of peak period as an example. After obtaining the new price, the membership function parameters of the elastic coefficients are obtained by Equation (6). The following relationship can be obtained from the expanded formula (1): where, E TOU.p represents demand at the peak period considering the uncertainty. In this formula, three fuzzy parameters of the electricity price elasticity coefficients and the initial demand are included, and the four fuzzy parameters have their respective membership function parameters, which avoids the approximation process of setting the initial demand as the historical average and ensures the universality and accuracy of the model. In addition, the four fuzzy parameters are linear with each other, which contributes to further processing. After obtaining the demand of peak period E TOU.p , the load consumption of each hour is obtained according to the ratio l t which can be obtained by the following formula: The opportunity constraint of the fuzzy constraint problem can be transformed into a crisp equivalence, and then the traditional solution method can be used to solve the transformed problem [24].
For the fuzzy variable δ and decision variable χ, its membership function is µ. If the function g(χ,δ) makes g(χ,δ) = h(χ)-δ true and when the fuzzy chance constraint is based on the credibility "Cr" as in Equation (22), it has a relationship as Equation (23) [24].
When α ≥ 1/2, the equality constraint (12) can be transformed into the crisp equivalence as follow: where, r 1 and r 2 are the membership parameters of the fuzzy parameter P L . However, P t L contains four fuzzy variables with a linear relationship. Therefore (24) can be converted into the following form: where,E p1 , E p2 , m xy1 , m xy2 (x = p, y = f,p,g) are membership function parameters. The process of transforming the original constraint into a crisp equivalence form deals with the uncertainty in the equality constraint, which greatly reduces the difficulty of the solution. So far, the day-ahead scheduling optimization model with the chance constraint has been transformed and can be solved by the traditional optimization method. In intraday optimization, robust interval variables are the focus of processing. Similar to fuzzy stochastic optimization, the general form of robust optimization is as follows: where, f (χ,ξ) is a objective function, g(χ,ξ) ≤ 0 is a constraint, andξ is a robust interval variable. Robust optimization methods can be roughly divided into static robust optimization methods and multi-stage robust optimization methods. Compared with the multi-stage optimization method, the results of the static robust optimization method tend to be conservative. Therefore, multi-stage optimization method is mainstream.
The robust interval variableξ can be expressed as follows: where, ξ represents the deterministic part ofξ and ξ represents the indeterminate part ofξ. Since the essence of robust optimization is to solve the optimal solution under the worst possibility, if g(χ,ξ) = h(χ) +ξ, then the Equation (27) is brought into Equation (26), which can be transformed into the following form [18]: For the same reason, the same processing can be performed on the intraday optimization problem with robust variableP t IL . In order to show the process more clearly, the fuzzy variable P t L will not be converted to a crisp equivalent. The intraday constraint is as follows: In order to minimize the risk of power shortage in the system, the following formula should be satisfied: Converting (29) to (30) is an expression of the robust idea that constraints should be still satisfied under harsh conditions. Equation (29) is equivalent to the following formula: Bringing (9) into the above formula: F is defined and should satisfy the following relationship: then (32) is converted into (34): Energies 2019, 12, 1711 10 of 26 G t is defined and satisfies the following relationship: and (28) can be converted into the following form: (14) and (15), Equation (19)   (36) It can be seen that Equation (36) should be a min-max optimization that minimizes the pessimistic value of the objective function. Since there are different decision variables in different optimizations of the min-max problem involved in this paper, it can be simplified to a two-stage optimization: after setting the robustness factor, the maximum value of the robust variables is obtained first, then, the solution to minimization problem of the intraday scheduling can be obtained, with the constraints being satisfied even in the worst case. A pseudo-code of PBDR and a simplified diagram of the IBDR are attached in Appendix C.

Model Solving
After the crisp equivalent conversion and transformation of robust optimization, the uncertainty optimization problem can be solved by deterministic optimization.
For the global optimization problem in the day-ahead stage, commonly used algorithms include particle swarm optimization (PSO), genetic algorithms (GA), etc. The bat algorithm (BA) is a heuristic intelligent optimization algorithm, which essentially uses the tuning technique to control the dynamic behavior of the bat group, and adjusts the relevant parameters to obtain the optimality [25]. Compared with other intelligent optimization algorithms, BA not only has the characteristics of less parameters and simple model, but also has advantages in solving [26]. It is often used for power system optimization scheduling [27][28][29].
The iterative formula and update formula of the algorithm are as follows [30]: where, f i is the pulse frequency of bat i, f min and f max are the upper and lower limits of the pulse frequency; s t i is the speed of bat i at t;x t i is the position of bat i at t;ρ is a uniformly distributed random variable on [0, 1]; x * is the current optimal bat position.
In order to process multi-objective functions in the day-ahead stage, entropy weights are introduced to transform multi-objective optimization into single-objective optimization [31]. The entropy of n evaluation indicators of m evaluation objects is as follows: where, v kX is the value of the indicator X(1,2, . . . ,m) of the object k(1,2, . . . ,n).
The entropy weight matrix of the evaluation index W and the weight w k are calculated as follows: There is no doubt, a weighted sum does not guarantee obtaining unsupported non-dominated solutions. However, although the Pareto optimal set can provide effective support for decision-making, the final decision still depends to a large extent on the designer's subjective judgment [32]. The weighted sum method combined with the entropy weight method is not only simple and efficient, but also has higher objectivity than the Pareto method [33,34]. The intraday phase is modeled by YALMIP and the CPLEX solution is called.

Study Data
In order to verify the effectiveness of the proposed model and algorithm, this paper simulates the example on the improved IEEE33 node system [35]. Among them, the B25-B32 node constitutes the microgrid part. The node B0 is the contact node of the distribution network and the superior distribution network. The load and grid structure parameters are shown in Tables A1-A6 in Appendix A; the installation cost, operation and maintenance cost, and output upper and lower limits of wind turbines, gas turbines, etc. are shown in Table A7 in Appendix B; the bat algorithm parameters are shown in Table A8 in Appendix B. The elasticity coefficient of electricity price in each period is shown in Table A9 in Appendix B; the membership function parameters of the baseline demand and the price elastic coefficient are shown in Table A10 in Appendix B; most of these data come from the literature [35]. As shown in Figure 1, the B3 node, the B7 node, the B12 node, the B14 node, and the B22 node each have their own LA, aggregating the load of each node to participate in the DR market. The maximum interrupt capacity is 15% of the load of each node and the interruption time of each node is less than one hour per day. The negotiated price of electricity purchased from the superior distribution network is 0.52 yuan/kWh; the unit price of IL compensation is 1 yuan/kWh; the unit price of the PLR is 1.2 yuan/kWh. This example is solved on the MATLAB (2015b, The MathWorks, Natick, MA, America) platform.  Table B1 in Appendix B; the bat algorithm parameters are shown in Table B2 in Appendix B. The elasticity coefficient of electricity price in each period is shown in Table B3 in Appendix B; the membership function parameters of the baseline demand and the price elastic coefficient are shown in Table B4 in Appendix B; most of these data come from the literature [35]. As shown in Figure 1, the B3 node, the B7 node, the B12 node, the B14 node, and the B22 node each have their own LA, aggregating the load of each node to participate in the DR market. The maximum interrupt capacity is 15% of the load of each node and the interruption time of each node is less than one hour per day. The negotiated price of electricity purchased from the superior distribution network is 0.52 yuan/kWh; the unit price of IL compensation is 1 yuan/kWh; the unit price of the PLR is 1. In the day-ahead scheduling, one day is divided into 24 periods on average; in the day-in scheduling stage, one day is divided into 96 periods on average, that is, the load and renewable energy output of the next period are predicted every 15 min and based on this, the operation is optimized. Table 1 shows the division of peak, flat and valley periods in Hebei Province, China. The load, wind turbine, photovoltaic output and purchased power by microgrid are shown as Figure 2.  In the day-ahead scheduling, one day is divided into 24 periods on average; in the day-in scheduling stage, one day is divided into 96 periods on average, that is, the load and renewable energy output of the next period are predicted every 15 min and based on this, the operation is optimized. Table 1 shows the division of peak, flat and valley periods in Hebei Province, China. The load, wind turbine, photovoltaic output and purchased power by microgrid are shown as Figure 2.

Influence of Uncertainty Model on Optimization Results
In order to verify the validity of the uncertainty model of PBDR and IBDR in the two-stage interactive decision model, this paper sets four scenarios to compare the operating cost in two stages: 1) Uncertainty is not considered in the day-ahead and intraday scheduling 2) Only consider the uncertainty of intraday scheduling 3) Only consider the uncertainty of the day-ahead scheduling 4) Uncertainty is considered in both day-ahead and intraday scheduling Since the time margin of intraday scheduling is small, the day-in reliability level should be higher [36]; on the contrary, the deviation of the day-ahead scheduling can be adjusted by the intraday scheduling, so that a smaller confidence can be taken in day-ahead stage. In addition, since the maximum robustness factor set in this paper is 5 and the minimum is 0, the robustness factor takes the intermediate value. Set the day-ahead reliability 0.6 α = , the day-in reliability = 0.9  [36], and the robustness coefficient Γ =3 ; the operating costs of two stages in the four scenarios are shown in the following Table 2. It can be seen that in the scenarios where uncertainty is considered in the previous stage, the scheduling cost is higher than other scenarios. This is because the uncertainty of price-based demand response is considered, and the estimation of the load change tends to be conservative, leading to the predicted load curve with poor effect of reducing the peak-to-valley difference. So the running cost is high.
However, considering uncertainty in the intraday scheduling phase reduces the operating costs of the intraday phase. This is because the uncertainty of the intraday phase is mainly due to the uncertainty of the IL, which can be considered as the uncertainty of outputs. And because the interrupted power will not be greater than the scheduling plan, so the dispatch center will increase

Influence of Uncertainty Model on Optimization Results
In order to verify the validity of the uncertainty model of PBDR and IBDR in the two-stage interactive decision model, this paper sets four scenarios to compare the operating cost in two stages: (1) Uncertainty is not considered in the day-ahead and intraday scheduling (2) Only consider the uncertainty of intraday scheduling (3) Only consider the uncertainty of the day-ahead scheduling (4) Uncertainty is considered in both day-ahead and intraday scheduling Since the time margin of intraday scheduling is small, the day-in reliability level should be higher [36]; on the contrary, the deviation of the day-ahead scheduling can be adjusted by the intraday scheduling, so that a smaller confidence can be taken in day-ahead stage. In addition, since the maximum robustness factor set in this paper is 5 and the minimum is 0, the robustness factor takes the intermediate value. Set the day-ahead reliability α = 0.6, the day-in reliability β = 0.9 [36], and the robustness coefficient Γ = 3; the operating costs of two stages in the four scenarios are shown in the following Table 2. It can be seen that in the scenarios where uncertainty is considered in the previous stage, the scheduling cost is higher than other scenarios. This is because the uncertainty of price-based demand response is considered, and the estimation of the load change tends to be conservative, leading to the predicted load curve with poor effect of reducing the peak-to-valley difference. So the running cost is high.
However, considering uncertainty in the intraday scheduling phase reduces the operating costs of the intraday phase. This is because the uncertainty of the intraday phase is mainly due to the uncertainty of the IL, which can be considered as the uncertainty of outputs. And because the interrupted power will not be greater than the scheduling plan, so the dispatch center will increase the planned amount of the IL. Although this has caused the increase of IL compensation cost, the output gap caused by the uncertainty of IL is reduced, leading to the reduction of called PLR. Therefore the total operating cost of the day-in stage is reduced.
In order to analyze the intraday cost in different scenarios, the outputs of IL, PLR and additional electricity purchased in the four scenarios during 11:00-17:00 are shown in Figure 3: as the Figure 3, during 11:00-12:00, the IL resources in the scenarios almost reached the maximum interrupt amount. This result from the uncertainty increases caused by the increasing demand in the peak period. Similarly, the additional electricity purchased is also larger in this period of the four scenarios. However, the occurrence of the peaks of the called amount of PLR are not specific to the peak period. This is due to the fact that PLR is called on account of the uncertainty of load and the IL. output gap caused by the uncertainty of IL is reduced, leading to the reduction of called PLR. Therefore the total operating cost of the day-in stage is reduced. In order to analyze the intraday cost in different scenarios, the outputs of IL, PLR and additional electricity purchased in the four scenarios during 11:00-17:00 are shown in Figure 3: as the Figure 3, during 11:00-12:00, the IL resources in the scenarios almost reached the maximum interrupt amount. This result from the uncertainty increases caused by the increasing demand in the peak period. Similarly, the additional electricity purchased is also larger in this period of the four scenarios. However, the occurrence of the peaks of the called amount of PLR are not specific to the peak period. This is due to the fact that PLR is called on account of the uncertainty of load and the IL.

Influence of Uncertainty Parameters
In the uncertainty optimization model, the change of confidence has a great influence on the optimization result, and generally the operating cost increases with the increase of confidence [21]. Because the increase of confidence means the reliability is improved. Important ways to improve reliability include measures to improve the reliability of the element or increase the redundancy of the system, so the economic cost will inevitably increase. In the two-stage joint scheduling model, not only the confidence, but also the robustness factor control uncertainty. How the two parameters affect the scheduling results jointly requires further exploration. After changing the robustness factor and the confidence level of intraday stage, the two-stage optimal scheduling is performed. The planned interrupted amount of the IL and the PLR output at every hour are obtained as in

Influence of Uncertainty Parameters
In the uncertainty optimization model, the change of confidence has a great influence on the optimization result, and generally the operating cost increases with the increase of confidence [21]. Because the increase of confidence means the reliability is improved. Important ways to improve reliability include measures to improve the reliability of the element or increase the redundancy of the system, so the economic cost will inevitably increase. In the two-stage joint scheduling model, not only the confidence, but also the robustness factor control uncertainty. How the two parameters affect the scheduling results jointly requires further exploration. After changing the robustness factor and the confidence level of intraday stage, the two-stage optimal scheduling is performed. The planned interrupted amount of the IL and the PLR output at every hour are obtained as in Figure 4. reliability include measures to improve the reliability of the element or increase the redundancy of the system, so the economic cost will inevitably increase. In the two-stage joint scheduling model, not only the confidence, but also the robustness factor control uncertainty. How the two parameters affect the scheduling results jointly requires further exploration. After changing the robustness factor and the confidence level of intraday stage, the two-stage optimal scheduling is performed. The planned interrupted amount of the IL and the PLR output at every hour are obtained as in Figure 4.
As shown in the Figure 4, when the intraday confidence is changed with determined robustness factor, the running cost increases as the day-in reliability increases. At the same time, in the case of day-in confidence being determined, increasing the robustness factor also leads to an increase in operating costs. This is because the essence of increasing the robustness factor is to increase the uncertainty, and the estimated power that cannot be interrupted by contract is increased, leading to an increase in PLR resources which affect the operating costs. In addition, by comparing Figures 4ac, it can be seen that compared with the robustness factor, the change of confidence has a greater impact on the result. This is mainly because the change of the robustness factor only affects the actual interruption of the IL. The selection of reliability will change the value of the load in the crisp equivalence constraint, thus the allocation of resources has to be changed.
When considering the day-ahead scheduling stage, the robustness factor and the confidence level of the day-ahead scheduling are changed, and the influence of the two parameters on the result is explored. When different robustness factors are adopted in the intraday scheduling phase, the dayin reliability α is 0.9, and the confidence of the day-ahead scheduling is changed respectively. The output of the PLR and the IL at every hour in the intra-day scheduling are as follows (see Figure 5):  As shown in the Figure 4, when the intraday confidence is changed with determined robustness factor, the running cost increases as the day-in reliability increases. At the same time, in the case of day-in confidence being determined, increasing the robustness factor also leads to an increase in operating costs. This is because the essence of increasing the robustness factor is to increase the uncertainty, and the estimated power that cannot be interrupted by contract is increased, leading to an increase in PLR resources which affect the operating costs. In addition, by comparing Figure 4a-c, it can be seen that compared with the robustness factor, the change of confidence has a greater impact on the result. This is mainly because the change of the robustness factor only affects the actual interruption of the IL. The selection of reliability will change the value of the load in the crisp equivalence constraint, thus the allocation of resources has to be changed.
When considering the day-ahead scheduling stage, the robustness factor and the confidence level of the day-ahead scheduling are changed, and the influence of the two parameters on the result is explored. When different robustness factors are adopted in the intraday scheduling phase, the day-in reliability α is 0.9, and the confidence of the day-ahead scheduling is changed respectively. The output of the PLR and the IL at every hour in the intra-day scheduling are as follows (see Figure 5): As can be seen from Figure 5, with the increase of confidence in the day-ahead stage, the output of IL and PLR are reduced significantly. This is because the estimation of load change after PBDR will be more conservative with the increase of confidence in day-ahead stage, which reduce the pressure on the output of the day-in stage. With the increase of the intra-day robustness factor, the output of the IL and PLR also increase slightly, which fully demonstrates the effectiveness of the two-stage interactive decision model. interruption of the IL. The selection of reliability will change the value of the load in the crisp equivalence constraint, thus the allocation of resources has to be changed.
When considering the day-ahead scheduling stage, the robustness factor and the confidence level of the day-ahead scheduling are changed, and the influence of the two parameters on the result is explored. When different robustness factors are adopted in the intraday scheduling phase, the dayin reliability α is 0.9, and the confidence of the day-ahead scheduling is changed respectively. The output of the PLR and the IL at every hour in the intra-day scheduling are as follows (see Figure 5): As can be seen from Figure 5, with the increase of confidence in the day-ahead stage, the output of IL and PLR are reduced significantly. This is because the estimation of load change after PBDR will be more conservative with the increase of confidence in day-ahead stage, which reduce the pressure

Influence of Weight of Objective Functions in Day-Ahead Phase
The day-ahead scheduling plan has a large impact on the operating costs. In this paper, in order to avoid damage to the user's interests caused by the extreme price, the objective functions of the day-ahead stage include user transfer coefficient in addition to the operating cost. The entropy weight method is used to determine the weight of the objective functions as [0.30, 0.70], and they are transformed into a single objective function. Changing the weight is equivalent to changing the objective function of the day-ahead scheduling phase. In order to explore the influence of the weight change on the optimization result, five different sets of objective function weights are taken to optimize.
It can be seen from Figure 6 that as the weight of the user transfer coefficient increases, the difference of price between the peak and flat period gradually decreases. And there is almost no difference when it comes to 0.9. This is because the pursuit of lower user transfer coefficient makes the optimal solution tend to the original price scheme, in which the difference between the flat price and the peak price is small. However, the load curve under this scheme has great peak-to-valley difference, which leads to an increase in operating costs. on the output of the day-in stage. With the increase of the intra-day robustness factor, the output of the IL and PLR also increase slightly, which fully demonstrates the effectiveness of the two-stage interactive decision model.

Influence of Weight of Objective Functions in Day-Ahead Phase
The day-ahead scheduling plan has a large impact on the operating costs. In this paper, in order to avoid damage to the user's interests caused by the extreme price, the objective functions of the day-ahead stage include user transfer coefficient in addition to the operating cost. The entropy weight method is used to determine the weight of the objective functions as [0.30, 0.70], and they are transformed into a single objective function. Changing the weight is equivalent to changing the objective function of the day-ahead scheduling phase. In order to explore the influence of the weight change on the optimization result, five different sets of objective function weights are taken to optimize.
It can be seen from Figure 6 that as the weight of the user transfer coefficient increases, the difference of price between the peak and flat period gradually decreases. And there is almost no difference when it comes to 0.9. This is because the pursuit of lower user transfer coefficient makes the optimal solution tend to the original price scheme, in which the difference between the flat price and the peak price is small. However, the load curve under this scheme has great peak-to-valley difference, which leads to an increase in operating costs. It can be seen from Figure 7 that as the weight of the operating cost increases and the weight of the transfer coefficient decreases, the operating cost decreases and the transfer coefficient increases. However, after the [0.5, 0.5], the running cost reduction rate slowed down, and the transfer coefficient increased rapidly. Excessive operating cost weights or transfer coefficient weights can cause another objective function to approach a worse situation. Among the weights of each group, [0.3, 0.7] can It can be seen from Figure 7 that as the weight of the operating cost increases and the weight of the transfer coefficient decreases, the operating cost decreases and the transfer coefficient increases. However, after the [0.5, 0.5], the running cost reduction rate slowed down, and the transfer coefficient increased rapidly. Excessive operating cost weights or transfer coefficient weights can cause another objective function to approach a worse situation. Among the weights of each group, [0.3, 0.7] can avoid the rapid increase of the transfer coefficient due to the weight change, and also ensure the lower operating cost. The result indicates the validity of the entropy weight method.

Comparison of Algorithms in Day-Ahead Phase
In order to verify the performance of BA, this paper compares the BA with the particle swarm optimization (PSO). The parameters of the BA are the same as Table B2 in Appendix; the acceleration factor of the PSO is set to c1 = c2 = 2, and the inertia weight w is set to 0.5. At the same time, the two entropy weights are the same as Table 3. Since the values of the two algorithms in the previous stage are quite different, they are normalized, and the formula is as follows.  are the maximum and minimum values in the initial solution. Keep the initial solution the same, running 10 times each. Table 4 is a comparison of the algorithms, and the convergence curves for one time are shown in Figure 8.

Comparison of Algorithms in Day-Ahead Phase
In order to verify the performance of BA, this paper compares the BA with the particle swarm optimization (PSO). The parameters of the BA are the same as Table A8 in Appendix A; the acceleration factor of the PSO is set to c 1 = c 2 = 2, and the inertia weight w is set to 0.5. At the same time, the two entropy weights are the same as Table 3. Since the values of the two algorithms in the previous stage are quite different, they are normalized, and the formula is as follows. F k.max and F k.min are the maximum and minimum values in the initial solution. Keep the initial solution the same, running 10 times each. Table 4 is a comparison of the algorithms, and the convergence curves for one time are shown in Figure 8.  As shown in Table 4 and Figure 8, among the two algorithms, BA has a stronger ability to solve, its optimal solution is better than the result of another method, and the standard deviation of the results of ten independent runs is lower; the time-consuming of PSO is shorter, but the standard deviation of the ten runs is large and that means lack of stability. Although the time spent by BA in optimization is higher, the stage in which BA is used does not have strict calculation speed requirements, so the application of BA with higher solution accuracy is reasonable.

Conclusions
This paper establishes a model of uncertainty for different types of DR, and establishes a twostage interactive decision model considering the uncertainty of DR and other output resources based on different speeds of response. The fuzzy chance constrained programming theory and multi-stage robust optimization are used to transform the uncertainty problem, and the BA and CPLEX are utilized to solve the problem. Finally, the following conclusions are verified in the IEEE33 node system: 1) Joint scheduling with different DR uncertainties at different time scales can effectively reduce the operating cost. But when the uncertainty of day-ahead stage is considered exclusively, the operating cost increases. 2) Considering the uncertainty of DR in the day-ahead stage will increase the cost of day-ahead, but reduce the cost of intraday scheduling, and the reduction will increase with the increase of the confidence of the day-ahead. 3) In the intraday phase, the operating cost increases with the increase of the robustness factor and the day-in reliability level. Among them, the day-in reliability has a greater impact on operating costs. 4) In the day-ahead stage, the excessive weight coefficient gap will affect the search for the optimal solution of the multi-objective solution. The entropy weight method can provide a more reasonable weight to avoid the deviation of an objective function from the ideal solution.
There are still many limitations in this article. Firstly, the uncertainty of renewable energy is not considered; secondly, the time scale of scheduling does not extend to the real-time stage. The next phase of work will be focused on an interactive model that considers the uncertainty on both sides of the source and the load.  As shown in Table 4 and Figure 8, among the two algorithms, BA has a stronger ability to solve, its optimal solution is better than the result of another method, and the standard deviation of the results of ten independent runs is lower; the time-consuming of PSO is shorter, but the standard deviation of the ten runs is large and that means lack of stability. Although the time spent by BA in optimization is higher, the stage in which BA is used does not have strict calculation speed requirements, so the application of BA with higher solution accuracy is reasonable.

Conclusions
This paper establishes a model of uncertainty for different types of DR, and establishes a two-stage interactive decision model considering the uncertainty of DR and other output resources based on different speeds of response. The fuzzy chance constrained programming theory and multi-stage robust optimization are used to transform the uncertainty problem, and the BA and CPLEX are utilized to solve the problem. Finally, the following conclusions are verified in the IEEE33 node system: (1) Joint scheduling with different DR uncertainties at different time scales can effectively reduce the operating cost. But when the uncertainty of day-ahead stage is considered exclusively, the operating cost increases. (2) Considering the uncertainty of DR in the day-ahead stage will increase the cost of day-ahead, but reduce the cost of intraday scheduling, and the reduction will increase with the increase of the confidence of the day-ahead. (3) In the intraday phase, the operating cost increases with the increase of the robustness factor and the day-in reliability level. Among them, the day-in reliability has a greater impact on operating costs. (4) In the day-ahead stage, the excessive weight coefficient gap will affect the search for the optimal solution of the multi-objective solution. The entropy weight method can provide a more reasonable weight to avoid the deviation of an objective function from the ideal solution.
There are still many limitations in this article. Firstly, the uncertainty of renewable energy is not considered; secondly, the time scale of scheduling does not extend to the real-time stage. The next phase of work will be focused on an interactive model that considers the uncertainty on both sides of the source and the load.

Conflicts of Interest:
The authors declare no conflict of interest.

Indexes t
The index of time j The index of IL user k The index of object x The index of the period of demand in m xy y The index of the period of price in m xy X The index of indicator i The index of bat Parameters T The number of operating cycles m xy The electricity price elasticity coefficient The initial demand of the peak, flat and valley periods c p , c f , c g The original electricity price at the peak, flat and valley periods m xy1 , m xy2 , m xy3 , m xy4 Membership function parameters of elastic coefficient Membership function parameters of baseline load k x , k y The growth rate of uncertainty M The number of IL users P mt.max /P mt.min The upper/lower limits of the unit's output γ up /γ down The upper/lower limit of the unit's output ramp rate P min.j /P max.j The lower and upper limits of the response D j The upper limit of the total response time ∆T The duration of a period f min /f max The upper/lower limits of the pulse frequency m The number of indicators n The number of objects F k.max /F k.min Maximum/minimum of the object k in the initial solution Set E TOU The electricity consumption column vector after ToU E 0 The column vector of the electricity consumption before PBDR M The electricity price elasticity matrix W The entropy weight matrix of the evaluation Variables ∆c p , ∆c f , ∆c g Electricity price change at the peak, flat and valley periods ∆E p , ∆E f , ∆E g Demand change at the peak, flat and valley periods E 0.x The baseline demand considering uncertainty m xy The elastic coefficient considering uncertainty E TOU.p , E TOU.f , E TOU.g The demand of the peak, flat and valley periods after PBDR µ E0.x Membership function of baseline load µ mxy Membership function of elastic coefficient P t

IL. j
The actual power of interruption of the j user at time t The current optimal bat position.

C k
The entropy of object k v kX The magnitude of the indicator X of the object k. w k Entropy weight of object k c 1 , c 2 The acceleration factor of the PSO w The    Figure B1. Elastic coefficient membership function.
The above is the membership function of the elastic coefficient of the x and y periods, and the membership function of the baseline load is the same. Step 7. Restore the load for each time period to each hour when the confidence is α for t = 1 to T if t = (8 to 12) or (16 to 20)