Next Article in Journal
Evaluation of the Use of Energy in the Production of Sweet Sorghum (Sorghum Bicolor (L.) Moench) under Different Production Systems
Previous Article in Journal
Flexible Modern Power System: Real-Time Power Balancing through Load and Wind Power
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

School of Electrical and Electronic Engineering, North China Electric Power University, Baoding 071000, China
*
Author to whom correspondence should be addressed.
Energies 2019, 12(9), 1711; https://doi.org/10.3390/en12091711
Submission received: 16 March 2019 / Revised: 29 April 2019 / Accepted: 30 April 2019 / Published: 6 May 2019

Abstract

:
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.

1. 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].
At present, there are many studies on optimal scheduling considering uncertainty [7,8,9,10,11,12]. The multi-scene method has been welcomed by many scholars because of its simplicity and convenience [7,8,9]. Monte Carlo method is often used to simulate the various possibilities of uncertainty [8]. 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.

2. 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.

2.1. 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:
E TOU = E 0 + [ E 0 . f 0 0 0 E 0 . p 0 0 0 E 0 . g ] M [ Δ c f / c f Δ c p / c p Δ c g / c g ]
where ETOU = ⎡ETOU.f, ETOU.p, ETOU.gT represents the column vector of demand after PBDR; ETOU.f, ETOU.p, and ETOU.g represent the demand of the peak, flat and valley periods after PBDR; E0 = ⎡E0.f, E0.p, E0.gT is the column vector of the demand in each period before PBDR; E0.f, E0.p and E0.g represent the initial demand of the peak, flat and valley periods; cf and Δcf represent the original electricity price and electricity price change at the flat period; cp and Δcp, cg and Δcg are similar. M is the electricity price elasticity matrix, and its expression is as follows:
M = ( m x y ) 3 × 3 , m x y = Δ E x E 0 . x ( Δ c y c y ) 1
where x,y = p,f,g, mxy 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. ΔEp, ΔEf, ΔEg 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:
{ E 0 . x = E ˜ 0 . x M = M ˜ M ˜ = ( m ˜ x y ) 3 x 3
where, E ˜ 0 . x and m ˜ x y 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:
μ E 0 . x = { 0 , E 0 . x < E x 1   or   E x 4 < E 0 . x E 0 . x E x 3 E x 4 E x 3 E x 3 < E 0 . x < E x 4 E 0 . x E x 1 E x 2 E x 1 , E x 1 < E 0 . x < E x 2 1 , E x 2 < E 0 . x < E x 3
μ m x y = { 0 , m x y < m x y 1 or m x y 4 < m x y m x y m x y 2 m x y 4 m x y 3 m x y 3 < m x y < m x y 4 m x y m x y 1 m x y 2 m x y 1 , m x y 1 < m x y < m x y 2 1 , m x y 2 < m x y < m x y 3
where mxy1, mxy2, mxy3, mxy4, Ex1, Ex2, Ex3, Ex4, are membership function parameters. Its function graph is shown in Appendix B. When mxy varies in [mxy2, mxy3], the function value is 1, and when mxy 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, kx and ky 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:
{ m x y 2 m x y 1 = m x y 4 m x y 3 m x y 2 m x y 1 = k x c x | Δ c x | + k y c y | Δ c y |
where, kx and ky represent the growth rate of uncertainty. When the electricity price changes greatly, mxy1 tends to mxy2, the uncertainty is small, the electricity price changes little, mxy1 is far away from mxy2, and the uncertainty becomes bigger.

2.2. 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:
P ^ IL . j t = P ¯ IL . j t + P IL . j t , j [ 1 , M ]
P IL . j t = τ j p IL . j t , τ j [ 1 , 1 ]
where, P ^ IL . j t indicates the actual interruptible power of the user j at time t, P ¯ IL . j t indicates the planned interrupted power of the user at that time, and P IL . j t indicates the indeterminate portion of the interrupted power of the user. p IL . j t 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 system P ^ IL . j t can be expressed as:
P ^ IL t = j = 1 M e j t ( P ¯ IL . j t + τ j p IL . j t ) , j = 1 M | τ j | Γ
where, M represents the number of users participating in IL; e j t 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.

3. 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.

3.1. Day-Ahead Scheduling

3.1.1. Objective Function

Operating Cost in Day-Ahead

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 day = t = 1 T ( c WT t P WT t + c pv t P pv t + c mt t P mt t + c buy t P buy t c pri t P MG t ) Δ T
where, T represents the number of operating cycles; c WT t and c pv t represent the unit cost of wind turbines and photovoltaic units in the period of t; P WT t and P pv t represent the output of wind turbine generation and photovoltaic units in the period. c mt t represents the unit cost of the gas turbine in the period, and P mt t represents the output of the gas turbine in the period. c buy t indicates the unit price of electricity purchased from the main network in the period, and P buy t indicates the power purchased from the main network in the period. c pri t represents the unit price of selling electricity to the microgrid in the period, and P MG t 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:
S = t = 1 T | P l o a d t P l o a d 0 t | t = 1 T P l o a d 0 t
where, S is the user transfer coefficient, P l o a d 0 t is the load in the period of t before the implementation of the time-of-use price, Plaodt is the load in the period after the implementation of the time-of-use.

3.1.2. Restrictions

In the period of t, there is the following relationship:
P ˜ L t + P MG t = P WT t + P pv t + P buy t + P mt t
where, P ˜ L t 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:
Cr { P ˜ L t + P MG t = P WT t + P pv t + P buy t + P mt t } > α
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:
P mt . min P mt t P mt . max
Δ T γ down < P mt t P mt t 1 < Δ T γ up
where, Pmt.max and Pmt.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.

3.2. Day-In Scheduling

3.2.1. 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:
C day in = t = 1 T ( c IL t P ¯ IL t + c def t Δ P buy t + c R t P R t ) Δ T
where P ¯ IL t represents the planned power of the IL in the period of t, c IL t represents the unit price of compensation for the IL in the period; P R t represents the output of peak load regulation (PLR) resources called in the period,   c R t represents the unit price at which the PLR resources is called, and Δ P buy t is the difference between the actual and the planned power of purchased electricity from main network, and should be a positive value, c def t is the unit price of compensation to the main network.

3.2.2. Restrictions

The system’s the power balancing constraints is as follows:
P ˜ L t + P MG t = P WT t + P pv t + P buy t + Δ P buy t + P mt t + P ^ IL t + P R t
where P ^ IL t 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:
Cr { P ˜ L t + P MG t = P WT t + P pv t + P buy t + Δ P buy t + P mt t + P R t + P ^ IL t } > β
User j who participates in an incentive-based demand response has the following constraints:
{ P min . j < P IL . j t < P max . j t = 1 T e j t Δ T = D j
where, Pmin.j and Pmax.j respectively represent the lower and upper limits of the response; Dj 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).

4. Model Transformation and Solution

4.1. Transformation of the Model of Uncertainty

4.1.1. 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):
E ˜ TOU . p = E ˜ 0 . p + E 0 . p m ˜ pf Δ c f c f + E 0 . p m ˜ pp Δ c p c p + E 0 . p m ˜ pg Δ c g c g
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 lt which can be obtained by the following formula:
l t = { P l o a d 0 t Δ T E 0 . p , t p e a k p e r i o d P l o a d 0 t Δ T E 0 . f , t f l a t p e r i o d P l o a d 0 t Δ T E 0 . g , t v a l l e y p e r i o d
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].
Cr { g ( χ , δ ) 0 } > α
h ( χ ) = { sup { K | K = μ 1 ( 2 α ) } , α 1 / 2 inf { K | K = μ 1 ( 2 2 α ) } , α > 1 / 2
When α ≥ 1/2, the equality constraint (12) can be transformed into the crisp equivalence as follow:
( 2 2 α ) r 1 + ( 2 α 1 ) r 2 = P WT t + P pv t + P buy t + P mt t P MG t
where, r1 and r2 are the membership parameters of the fuzzy parameter P ˜ L . However, P ˜ L t contains four fuzzy variables with a linear relationship. Therefore (24) can be converted into the following form:
l t Δ T { ( 2 2 α ) [ E p 1 + E 0 . p ( m pf 1 Δ c f c f + m pp 1 Δ c p c p + m pg 1 Δ c g c g ) ] + ( 2 α 1 ) [ E p 2 + E 0 . p ( m pf 2 Δ c f c f + m pp 2 Δ c p c p + m pg 2 Δ c g c g ) ] } = P pv t + P WT t + P buy t + P mt t P MG t
where,Ep1, Ep2, mxy1, mxy2 (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.

4.1.2. Transformation of Uncertainty Model of IBDR

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:
{ min   f ( χ , ξ ^ ) s . t . g ( χ , ξ ^ ) 0
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]:
{ min   f ( χ , ξ ^ ) s . t . h ( χ ) + ξ ¯ + max ( ξ ) 0
For the same reason, the same processing can be performed on the intraday optimization problem with robust variable P ^ IL t . In order to show the process more clearly, the fuzzy variable P ˜ L t will not be converted to a crisp equivalent. The intraday constraint is as follows:
P ˜ L t = P WT t + P pv t + P buy t + Δ P buy t + P mt t + P ^ IL t + P R t P MG t
In order to minimize the risk of power shortage in the system, the following formula should be satisfied:
P ˜ L t + P MG t P ^ IL t P WT t + P pv t + P buy t + Δ P buy t + P R t + P mt t
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:
P ˜ L t + P MG t P ^ IL t P WT t P pv t P buy t Δ P buy t P mt t P R t 0
Bringing (9) into the above formula:
P ˜ L t j = 1 M e j t ( P ¯ IL . j t + τ j p IL . j t ) + P MG t P WT t P pv t P buy t Δ P buy t P mt t P R t 0
F ^ is defined and should satisfy the following relationship:
F ^ = j = 1 M e j t τ j p IL . j t
then (32) is converted into (34):
P ˜ L t j = 1 M e j t P ¯ IL . j t + max ( F ^ ) + P MG t P WT t P pv t P buy t Δ P buy t P mt t P R t 0
G t is defined and satisfies the following relationship:
G t = P MG t j = 1 M e j t P ¯ IL . j t P WT t P pv t P buy t Δ P buy t P mt t P R t
and (28) can be converted into the following form:
{ min C day in s . t . max ( F ^ ) + P ˜ L t + G t 0 Equations ( 14 )   and   ( 15 ) , Equation ( 19 )
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.

4.2. 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]:
{ f i = f min + ( f max f min ) ρ s i t = s i t 1 + ( x * x i t ) f i x i t = x i t 1 + s i t
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 i t is the speed of bat i at t ; x i t 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:
C k = X = 1 m B k X l n B k X / l n m , B k X = v k X / X = 1 m v k X
where, v k X 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:
W = ( w k ) 1 x n , w k = ( 1 C k ) / ( n k = 1 n C k ) , k = 1 n w k = 1
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.

5. Case Study and Discussion

5.1. 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 Table A1, Table A2, Table A3, Table A4, Table A5 and Table 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.
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.

5.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 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.

5.3. 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.
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.

5.4. 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 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.

5.5. 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; 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. 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.
N k = F k . max F k F k . max F k . min , F k . min < F k < F k . max
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.

6. 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.

Author Contributions

Conceptualization, S.S. and Q.G.; methodology, S.S. and Q.G.; software, Q.G.; validation, S.S. and Q.G.; formal analysis, S.S. and Q.G.; investigation, S.S. and Q.G.; resources, S.S.; data curation, S.S.; writing—original draft preparation, Q.G.; writing—review and editing, S.S.; visualization, Q.G.; supervision, S.S.; project administration, S.S.; funding acquisition, S.S.

Acknowledgments

This work was supported by State Key Laboratory of Alternate Electrical Power System with Renewable Energy Sources (NCEPU).

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Indexes
tThe index of time
jThe index of IL user
kThe index of object
xThe index of the period of demand in mxy
yThe index of the period of price in mxy
XThe index of indicator
iThe index of bat
Parameters
TThe number of operating cycles
mxyThe electricity price elasticity coefficient
E0.p, E0.f, E0.gThe initial demand of the peak, flat and valley periods
cp, cf, cgThe original electricity price at the peak, flat and valley periods
mxy1, mxy2, mxy3, mxy4Membership function parameters of elastic coefficient
Ex1, Ex2, Ex3, Ex4Membership function parameters of baseline load
kx,kyThe growth rate of uncertainty
MThe number of IL users
Pmt.max/Pmt.minThe upper/lower limits of the unit’s output
γup/γdownThe upper/lower limit of the unit’s output ramp rate
Pmin.j/Pmax.jThe lower and upper limits of the response
DjThe upper limit of the total response time
ΔTThe duration of a period
fmin/fmaxThe upper/lower limits of the pulse frequency
mThe number of indicators
nThe number of objects
Fk.max/Fk.minMaximum/minimum of the object k in the initial solution
Set
ETOUThe electricity consumption column vector after ToU
E0The column vector of the electricity consumption before PBDR
MThe electricity price elasticity matrix
WThe entropy weight matrix of the evaluation
Variables
Δcp, Δcf, ΔcgElectricity price change at the peak, flat and valley periods
ΔEp, ΔEf, ΔEgDemand change at the peak, flat and valley periods
E ˜ 0 . x The baseline demand considering uncertainty
m ˜ x y The elastic coefficient considering uncertainty
ETOU.p, ETOU.f, ETOU.gThe demand of the peak, flat and valley periods after PBDR
μE0.xMembership function of baseline load
μmxyMembership function of elastic coefficient
P ^ IL . j t The actual power of interruption of the j user at time t
P ¯ IL . j t The planned interrupted power of the j user at time t
P IL . j t The indeterminate portion of the interrupt capacity of the user at time t
p IL . j t Upper limit of the absolute value of the response uncertainty of the user j at time t
τ j The uncertainty coefficient of the j
P ^ IL t The actual interrupted power of the system at time t
e j t The state of user j during t period
Γ Robustness factor
CdayThe cost in day-ahead stage
c WT t The unit cost of wind turbine generation at t
P WT t The output of wind turbine generation at t
c pv t The unit cost of photovoltaic units at t
P pv t The output of photovoltaic units at t
c mt t The unit cost of the gas turbine at t
P mt t The output of the gas turbine at t
c buy t The unit purchasing price of electricity at t
P buy t The purchasing power at t
c pri t The unit price of electricity sold to microgrid at t
P MG t The power purchased by the microgrid at t
SThe user transfer coefficient
Pload0tThe load of the time t before ToU
PloadtThe load of the time t after ToU
P ˜ L t The fuzzy parameter of the load after ToU
αThe confidence level in day-ahead stage
Cday-inThe operating cost during the intraday
βThe confidence level in day-in stage
c IL t The unit price of the IL at t
c R t The unit price of PLR
c def t The unit price of compensation to the main network
P R t The output of PLR at t
Δ P buy t The difference from the actual and the planned amount of purchased power
P ¯ IL t The planned power of the IL at time t
E ˜ TOU . p The demand at the peak period considering the uncertainty
ltRatio of hourly demand to demand of per period
χ Generalization of decision variable
δ Generalization of fuzzy variable
μMembership function of δ
ξ ^ Generalization of robust interval variables
ξ ¯ The deterministic part of ξ ^
ξ The indeterminate part of ξ ^
fiThe pulse frequency of bat i
s i t The speed of bat i at t
x i t The position bat i at t
ρ Uniformly distributed random variable on [0, 1]
x*The current optimal bat position.
CkThe entropy of object k
vkXThe magnitude of the indicator X of the object k.
wkEntropy weight of object k
c1, c2The acceleration factor of the PSO
wThe inertia weight of the PSO
FkFunction value that should be normalized
NkNormalized Fk
Abbreviations
DRDemand response
PBDRPrice-based demand response
IBDRIncentive-based demand response
ILInterruptible load
BABat algorithm
DLCDirect load control
DSBDemand side bid
EDREmergency demand response
WTWind turbine
MTMicro turbine
PLRPeak load regulation
TOUTime of Use
C&CGColumn-and-constraint generation
LALoad aggregator

Appendix A

Table A1. The bus load at different time of nodes 3–9.
Table A1. The bus load at different time of nodes 3–9.
TimeB3B4B6 B7B8B9
P/kWQ/kVarP/kWQ/kVarP/kW Q/kVarP/kWQ/kVarP/kWQ/kVarP/kWQ/kVar
198.2545.2352.3225.6974.3936.53101.3453.6978.6940.4728.5914.23
290.5541.6848.2223.6871.1634.9496.9451.3675.2738.7127.3513.61
386.6739.946.1522.6670.4234.5895.9350.8274.4938.3127.0613.47
482.8738.1444.1321.6669.6134.1894.8350.2473.6337.8726.7513.32
587.6940.0546.6922.7568.5533.6693.3849.4772.5137.2926.3413.11
6108.9749.7758.0328.2774.8636.76101.9854.0379.1940.7328.7714.32
7126.9357.9767.5932.9395.8847.08130.6169.2101.4252.1636.8518.34
8147.6467.4478.6238.3107.6952.89146.7177.73113.9258.5941.3920.6
9162.4474.1986.542.14119.1158.49162.2685.97125.9964.845.7822.78
10170.0177.6590.5344.1128.7163.2175.3492.9136.1570.0249.4724.62
11181.4782.8996.6347.08136.5967.07186.0798.58144.4874.3152.4926.13
12168.9377.1789.9643.83160.2378.68218.28115.64169.4987.1761.5830.65
13168.9277.1789.9543.83152.3474.81207.54109.96161.1582.8858.5529.14
14165.1475.4487.9442.85128.7163.2175.3392.9136.1570.0249.4624.62
15174.3979.6792.8745.25130.0863.88177.2193.89137.670.7749.9924.88
16175.4880.1693.4545.53135.2766.43184.2897.63143.0973.5951.9925.88
17176.5580.6794.0245.82131.3364.49178.9194.79138.9271.4550.4725.12
18172.1678.6691.6844.68144.4670.94196.8104.27152.8178.655.5227.64
19176.5780.6794.0345.82152.3474.81207.53109.96161.1482.8858.5529.14
20210.4596.16112.0754.62145.7771.59198.58105.22154.279.3156.0227.89
21211.5896.67112.6754.91136.5867.07186.0698.58144.4774.3152.4926.13
22172.1978.6791.6944.69133.9565.78182.4896.69141.6972.8851.4825.63
23126.3557.7367.2832.79122.1359.98166.3888.16129.1966.4546.9423.36
24103.9247.4855.3426.9794.5546.44128.8168.25100.0251.4436.3418.09
Table A2. The bus load at different time of nodes 10–15.
Table A2. The bus load at different time of nodes 10–15.
TimeB10B11B12B13 B14B15
P/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVar
153.6726.7432.1415.4252.9727.12150.7570.7550.2520.32124.2555.35
251.3425.5830.7414.7550.6725.94140.7567.4245.2519.36113.2852.75
350.825.3130.4214.650.1425.67128.7564.0543.2518.4105.7350.11
450.2225.0230.0714.4349.5725.38115.4559.6339.2417.13102.7846.65
549.4624.6429.6214.2148.8124.99112.8956.2937.2616.1799.5644.04
654.0126.9132.3415.5253.3127.29107.8954.2134.2615.5795.5642.41
769.1734.4641.4219.8768.2734.95106.353.3433.2615.3294.6641.73
877.738.7146.5322.3276.6839.26118.2561.3541.5417.62105.3847.99
985.9342.8151.4624.6984.8143.42140.9574.4760.8821.39121.4758.26
1092.8646.2755.6126.6891.6546.92150.8581.0368.9523.27132.3363.39
1198.5449.159.0128.3197.2649.8156.8184.2771.7624.2137.6965.93
12115.657.669.2333.21114.0958.41149.7582.0369.8522.27130.5363.59
13109.9154.7665.8231.58108.4855.54153.0484.4971.3922.94133.465.5
1492.8646.2755.6126.6891.6546.92163.1490.0576.124.45142.2169.81
1593.8546.7656.226.9792.6347.43161.0288.9375.1124.15140.3668.94
1697.5948.6358.4428.0496.3249.32149.1482.2669.122.34129.1363.77
1794.7547.2156.7427.2293.5247.88148.2581.8268.6922.21128.7563.43
18104.2351.9362.4229.95102.8752.67134.8274.3562.4220.19117.6557.64
19109.9154.7665.8231.58108.4755.54135.6174.863.7920.31117.3457.99
20105.1752.462.9830.22103.853.15148.2381.8169.6822.21128.3463.42
2198.5449.159.0128.3197.2549.8164.2290.6576.6324.61142.7770.27
2296.6448.1557.8727.7795.3848.84154.9785.3972.8223.18134.2166.2
2388.1143.9152.7725.3286.9644.53145.7380.1668.5121.77126.1562.15
2468.2233.9940.8519.667.3334.47138.0575.9264.9320.61119.4758.85
Table A3. The bus load at different time of nodes 16–21.
Table A3. The bus load at different time of nodes 16–21.
Time B16B17B18B19B20B21
P/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVar
190.7540.2374.2530.78124.7550.256025.6845.2821.5992.144.83
288.7738.3472.6729.33123.7547.8959.6524.4741.7319.984.8841.32
387.2536.4270.5827.87121.7445.4955.9323.2539.9419.0481.2439.54
485.8633.9165.5825.94108.4342.3551.7821.6538.1918.2177.6837.8
581.5432.0160.624.49105.4839.9845.3220.4340.4119.1282.239.7
679.5430.8357.623.58100.8538.544.919.6850.2223.76102.1549.33
777.6430.3356.323.2199.337.8944.319.3658.527.67118.9857.46
887.9434.8867.6526.69111.243.5754.822.2768.0432.19138.466.84
9105.7542.3586.8532.4129.2352.8969.5727.0374.8635.42152.2773.54
10115.3546.0795.7535.25138.8357.5578.8229.4178.3537.07159.3776.96
11119.9647.9299.5836.66144.3859.8581.9730.5983.6339.57170.1182.16
12117.4546.3792.6535.75140.6358.5579.9529.6177.8536.84158.3576.49
13121.0347.7794.6936.82143.7260.3181.7130.577.8536.83158.3576.48
14127.9650.91101.9439.25153.2164.2787.132.5176.1136.01154.874.77
15126.2950.2899.6338.76151.2263.4885.9732.180.3738.03163.4778.97
16117.1946.5192.6635.85139.1258.7280.0929.6980.8738.26164.579.45
17116.4946.2591.8135.66138.4958.479.6229.5381.3738.5165.579.95
18104.9642.0483.832.41125.6753.0772.4526.8479.3437.55161.3877.96
19106.5842.2983.2932.6127.4253.3971.872781.3738.51165.5279.96
20115.4846.2592.0935.66138.5758.479.6129.5396.9945.9197.2895.31
21129.8251.25100.9939.51154.2864.787.1432.7297.5146.15198.3495.82
22121.7648.2896.2637.22145.5960.9582.230.8279.3637.55161.4177.98
23114.5145.3289.5434.94136.9257.2277.2728.9458.2327.56118.4457.22
24108.542.9284.833.09129.7154.1974.1727.447.8922.6697.4147.06
Table A4. The bus load at different time of nodes 22–27.
Table A4. The bus load at different time of nodes 22–27.
Time B22B23B24B25 B26B27
P/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVar
192.542.4920.39.4640.9421.77101.3552.3641.6923.5876.5839.87
285.2539.1618.718.7237.7320.0697.0550.1439.9222.5873.3338.18
381.637.4817.918.3436.1119.288.0345.4836.2120.4866.5234.63
478.0235.8317.127.9834.5318.3685.4644.1535.1519.8864.5733.62
582.5537.6218.128.3836.5419.2882.0242.3733.7419.0861.9732.27
6102.5946.7622.5110.4145.4123.9684.1343.4634.6119.5763.5733.1
7119.554.4626.2212.1352.8927.990.9446.9837.4121.1668.7235.78
813963.3530.514.161.5232.46103.3253.3842.524.0478.0740.65
9152.9369.733.5615.5267.6935.71126.4865.3452.0329.4395.5649.75
10160.0672.9535.1316.2470.8437.38140.9872.8357.9932.8106.5255.46
11170.8577.8737.4917.3475.6239.9158.1481.765.0536.79119.4962.21
12159.0472.4934.916.1470.3937.14156.8481.0364.5236.49118.5161.7
13159.0372.4934.916.1470.3937.14127.3165.7752.3729.6296.250.08
14155.4870.8734.1215.7868.8136.31134.5869.5355.3631.31101.6952.94
15164.1874.8536.0316.6672.6738.35133.3368.8854.8431.02100.7452.45
16165.2175.336.2616.7773.1238.58142.3273.5358.5433.11107.5455.99
17166.2275.7836.4816.8773.5738.83143.1673.9658.8933.31108.1756.32
18162.0873.8935.5716.4571.7437.86142.7473.7458.7233.21107.8656.15
19166.2475.7936.4816.8773.5738.83133.7269.0955.0131.11101.0452.61
20198.1490.3343.4820.1187.6946.28147.4576.1860.6534.31111.4158.01
21199.290.8243.7220.2288.1746.53165.0685.2767.938.4124.7264.93
22162.1173.9135.5816.4571.7537.87142.7373.7458.7133.21107.8556.15
23118.9554.2326.1112.0752.6527.79123.5163.8150.8128.7493.3348.59
2497.8344.621.479.9343.322.85113.1858.4746.5626.3385.5244.53
Table A5. The bus load at different time of nodes28–32.
Table A5. The bus load at different time of nodes28–32.
Time B28B29B30B31B32
P/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVarP/kW Q/kVar
144.6923.59105.2450.3773.6853.6942.4920.4795.6649.28
242.7922.59100.7848.2370.5651.4140.6919.691.647.19
338.8220.4991.4143.756446.6436.9117.7883.0942.81
437.6819.8988.7442.4762.1345.2735.8317.2680.6641.55
536.1719.0985.1740.7659.6343.4534.3916.5777.4139.88
637.119.5887.3641.8161.1644.5735.2716.9979.440.91
740.121.1794.4345.266.1148.1838.1318.3785.8444.22
845.5624.05107.2951.3575.1154.7343.3220.8797.5250.24
955.7729.44131.3362.8691.956753.0225.54119.3761.5
1062.1632.81146.3970.07102.4974.6859.128.47133.0668.55
1169.7336.81164.2178.6114.9783.7866.331.94149.2676.9
1269.1636.51162.8677.95114.0283.0965.7531.68148.0476.26
1356.1429.63132.263.2792.5567.4453.3725.71120.1661.9
1459.3431.32139.7466.8997.8471.2956.4227.18127.0265.44
1558.7931.03138.4566.2696.9370.6355.926.93125.8464.83
1662.7633.13147.7970.73103.4775.459.6728.75134.3369.2
1763.1333.32148.6671.15104.0875.8460.0228.91135.1269.61
1862.9433.22148.2270.94103.7775.6259.8428.83134.7369.41
1958.9631.13138.8666.4697.2270.8456.0627.01126.2265.02
2065.0234.32153.1173.28107.1978.1161.8229.78139.1771.7
2172.7838.42171.3982.03119.9987.4469.233.34155.7980.26
2262.9433.22148.2170.94103.7675.6159.8428.83134.7269.4
2354.4628.75128.2561.3889.7965.4351.7824.95116.5860.06
2449.9126.34117.5356.2582.2859.9647.4522.86106.8355.03
Table A6. The resistance parameters of the IEEE 33 network.
Table A6. The resistance parameters of the IEEE 33 network.
Branch Number Starting Node End Node ResistanceReactanceBranch Number Starting Node End Node ResistanceReactance
/Ω /Ω /Ω /Ω
1B0 B1 0.09220.04717B23 B24 0.7860.564
2B1 B13 0.4930.251118B5 B6 1.0590.9337
3B13 B14 0.1640.156519B6 B7 1.030.74
4B14 B15 0.45120.308320B7 B8 0.80420.7006
5B15 B16 0.3660.186421B8 B9 1.0440.74
6B16 B17 1.5041.355422B9 B10 0.50750.2585
7B17 B18 0.38110.194123B10 B11 0.19660.065
8B18 B19 0.40950.478424B11 B12 0.97440.963
9B1 B2 0.8960.701125B5 B25 0.37440.1238
10B2 B3 0.8190.70726B25 B26 0.31050.3619
11B3 B4 0.70890.937327B26 B27 1.4681.115
12B4 B5 0.2030.103428B25 B28 0.3410.5302
13B2 B20 0.18720.618829B28 B29 0.54120.7129
14B20 B21 0.28420.144730B28 B30 0.5910.526
15B21 B22 0.71140.235131B30 B31 0.74630.545
16B22 B23 0.7320.57432B31 B32 1.2891.721

Appendix B

Table A7. Unit parameters.
Table A7. Unit parameters.
Unit TypePower OutputYears of UseInstallation Cost (Wanyuan)Operation & Maintenance (yuan/kW·h)
Lower LimitUpper Limit
WT0700102.3750.0296
PT0500206.650.0096
MT0350151.6670.03
Table A8. Parameters of Bat algorithm.
Table A8. Parameters of Bat algorithm.
ParametersPopulation SizeNumber of IterationsParametersFrequencyPulse Frequency Increase Factor
data100100data[0,2]0.9
ParametersPulse LoudnessLoudness Attenuation CoefficientParametersMutation ProbabilityRandom Coefficient
data[1,2]0.9data0.90.5
Figure A1. Elastic coefficient membership function.
Figure A1. Elastic coefficient membership function.
Energies 12 01711 g0a1
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.
Table A9. Electricity price elasticity coefficient.
Table A9. Electricity price elasticity coefficient.
TimePeakFlatValley
Peak−0.01730.03830.052
Flat0.0192−0.04260.0577
Valley0.02020.0448−0.0607
Table A10. The membership function parameters of the baseline load and the price elastic coefficient.
Table A10. The membership function parameters of the baseline load and the price elastic coefficient.
Parametersmxy2mxy3kxkyEx1Ex2Ex3Ex4
Data95%mxy105%mxy1%mxy1%mxy92%mxy95%mxy105%mxy108%mxy

Appendix C

Appendix C.1. Algorithm of PBDR:

Input: original load, original electricity price, price elastic coefficient matrix, electricity price difference generated in this iteration of BA, baseline load and elasticity coefficient membership parameter, uncertainty growth rate
Output: hourly load after response
Procedure:
Step 1. Obtain the membership function parameters of the electricity price elastic coefficient according to formula (6):
{ m x y 2 m x y 1 = m x y 4 m x y 3 m x y 2 m x y 1 = k x c x | Δ c x | + k y c y | Δ c y |
Step 2. Initialize the load of each period
E 0 . f = 0 ; E 0 . p = 0 ; E 0 . g = 0 ;
Step 3. Obtain the total load of peak time, valley time and flat time:
for t = 1 to T
  if t = (8 to 12) or (16 to 20), E 0 . p =   E 0 . p + P l o a d 0 t Δ T ;
  else if t = (22 to 24) or (1 to 6), E 0 . g =   E 0 . g + P l o a d 0 t Δ T ;
   else E 0 . f =   E 0 . f + P l o a d 0 t Δ T ;
  end if
end for
Step 4. Calculate the baseline load for each period when the confidence is according to Equation (23):
E ˜ 0 . f = ( 2 2 α ) E f 1 + ( 2 α 1 ) E f 2
E ˜ 0 . g = ( 2 2 α ) E g 1 + ( 2 α 1 ) E g 2
E ˜ 0 . p = ( 2 2 α ) E p 1 + ( 2 α 1 ) E p 2
Step5. Obtain the elastic coefficient when the confidence is α according to Equation (23):
for x = 1 to 3
  for y = 1 to 3
   m ˜ xy = ( 2 2 α ) m xy 1 + ( 2 α 1 ) m xy 2
  end for
end for
Step 6. Obtain the load of periods using outcome of step 4 and 5 according to Equation (1):
E ˜ TOU . f = E ˜ 0 . f + E 0 . f m ˜ ff Δ c f c f + E 0 . f m ˜ fp Δ c p c p + E 0 . f m ˜ fg Δ c g c g
E ˜ TOU . p = E ˜ 0 . p + E 0 . p m ˜ pf Δ c f c f + E 0 . p m ˜ pp Δ c p c p + E 0 . p m ˜ pg Δ c g c g
E ˜ TOU . g = E ˜ 0 . g + E 0 . g m ˜ gf Δ c f c f + E 0 . g m ˜ gp Δ c p c p + E 0 . g m ˜ gg Δ c g c g
Step7. 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)
     P ˜ L t = E ˜ TOU . p P l o a d 0 t E 0 . p
  else if t = (22 to 24) or (1 to 6)
     P ˜ L t = E ˜ TOU . f P l o a d 0 t E 0 . f
  else
     P ˜ L t = E ˜ TOU . g P l o a d 0 t E 0 . g
  end if
end for

Appendix C.2. Sketch Map of IBDR:

Figure A2. Sketch map of IBDR.
Figure A2. Sketch map of IBDR.
Energies 12 01711 g0a2

References

  1. Notice on Printing and Distributing the Measures for Power Demand Side Management. Available online: http://bgt.ndrc.gov.cn/zcfb/201011/t20101116_498818.html (accessed on 4 November 2010).
  2. Meng, X.; Cao, J.; Wang, Z.; Du, W. Construction and Energy Storage Mode Analysis of Energy Interconnected Microgrid Multi-energy Complementary System. Proc. CSEE 2018, 38, 5727–5737. [Google Scholar]
  3. Zhang, L.; Yuan, Y.; Sun, D.; Yuan, X.; Li, Q.; Su, D. Based wind farm joint operation scheduling model based on two-stage robust interval optimization. Electr. Power Automat. Equip. 2018, 38, 59–66. [Google Scholar]
  4. Zhang, Y.; Ren, S.; Yang, X.; Bao, W.; Xie, L.; Qi, J. Independent Microgrid Optimization Configuration Considering Price Demand Response. Electr. Power Automat. Equip. 2017, 37, 55–62. [Google Scholar]
  5. Cao, W.; Li, B.; Zheng, A.; Qi, B.; Yang, Z.; Su, Y. The Architecture and Technology of Demand Response within Energy Internet. In Proceedings of the 2016 China International Conference on Electricity Distribution (CICED), Xi’an, China, 10–13 August 2016. [Google Scholar]
  6. Niu, W.; Li, Y.; Wang, L. Based on risk assessment and opportunity constraints for uncertain interruptible load optimization scheduling. Electr. Power Automat. Equip. 2016, 36, 62–67. [Google Scholar]
  7. Huang, H.; Li, F.X.; Mishra, Y. Modeling Dynamic Demand Response Using Monte Carlo Simulation and Interval Mathematics for Boundary Estimation. IEEE Trans. Smart Grid 2015, 6, 2704–2713. [Google Scholar] [CrossRef]
  8. Zhu, L.; Liu, S.; Tang, Y.; Wang, J.; Huang, C. Modeling of charge and discharge uncertainty response and electric vehicle dealer’s dispatching strategy. Power Syst. Technol. 2018, 42, 3305–3317. [Google Scholar]
  9. Wang, B.; Sun, Y.; Li, Y. Application of Uncertainty Demand Response Modeling in Power Integral Incentive Decision. Automat. Electr. Power Syst. 2015, 39, 93–99. [Google Scholar]
  10. Zhou, B.; Lv, L.; Gao, H.; Tan, X.; Wu, H. virtual power plant optimization trading strategy based on two-stage stochastic programming. Power Constr. 2018, 39, 70–77. [Google Scholar]
  11. Yu, H.; Chung, C.; Wong, K.; Lee, H.; Zhang, J. Probabilistic Load Flow Evaluation with Hybrid Latin Hypercube Sampling and Cholesky Decomposition. IEEE Trans. Power Syst. 2009, 24, 661–667. [Google Scholar] [CrossRef]
  12. Meng, A.; Lin, Y.; Yin, H. Considers the uncertainty of the factor of household grid-wind-light-storage collaborative economic dispatch optimization method. Power Syst. Technol. 2018, 42, 162–173. [Google Scholar]
  13. Good, N.; Karangelos, E.; Navarro-Espinosa, A.; Mancarella, P. Optimization Under Uncertainty of Thermal Storage-Based Flexible Demand Response With Quantification of Residential Users’ Discomfort. IEEE Trans. Smart Grid 2015, 6, 2333–2342. [Google Scholar] [CrossRef]
  14. Niu, W.; Li, Y.; Wang, B. Considers uncertainty in demand response virtual power plant modeling. Proc. CSEE 2014, 34, 3630–3637. [Google Scholar]
  15. Zeng, D.; Yao, J.; Yang, S.; Wang, W.; Zhou, J.; Shi, F. Probabilistic power flow calculation with taking into account price-type load response uncertainty. Automat. Electr. Power Syst. 2015, 39, 66–71. [Google Scholar]
  16. Zhang, M.; Hu, Z.; Li, Y.; Xie, S. Based on feasibility test for unit combination robust optimization considering wind power and demand response. Proc. CSEE 2018, 38, 3184–3194. [Google Scholar]
  17. Yi, J.; Lyons, P.F.; Davison, P.J.; Wang, P.; Taylor, P.C. Robust Scheduling Scheme for Energy Storage to Facilitate High Penetration of Renewables. IEEE Trans. Sustain. Energy 2016, 7, 797–807. [Google Scholar] [CrossRef]
  18. Yi, W.; Zhang, Y.; Zhao, Z.; Huang, Y. Multiobjective Robust Scheduling for Smart Distribution Grids: Considering Renewable Energy and Demand Response Uncertainty. IEEE Access 2018, 6, 45715–45724. [Google Scholar] [CrossRef]
  19. Zeng, M.; Yang, Y.; Xiang, H.; Wang, L.; Sun, J. Power System Robust Optimization Programming Model with Taking Demand Side Response. Automat. Electr. Power Syst. 2016, 40, 137–145. [Google Scholar]
  20. Zhao, C.; Wang, J.; Watson, J.; Guan, Y. Multi-Stage Robust Unit Commitment Considering Wind and Demand Response Uncertainties. IEEE Trans. Power Syst. 2013, 28, 2708–2717. [Google Scholar] [CrossRef]
  21. Luo, C.; Li, Y.; Xu, H.; Li, L.X.; Hou, T.; Miao, S. Analysis of the Influence of Demand Response Uncertainty on Day-Time Optimization Scheduling. Automat. Electr. Power Syst. 2017, 41, 22–29. [Google Scholar]
  22. Liu, X.; Gao, B.; Li, Y. Bayesian Game-Theoretic Bidding Optimization for Aggregators Considering the Breach of Demand Response Resource. Appl. Sci. 2019, 9, 576. [Google Scholar] [CrossRef]
  23. Gao, Y.; Sun, Y.; Wang, X.; Chen, F.; Ehsan, A.; Li, H.; Li, H. Multi-Objective Optimized Aggregation of Demand Side Resources Based on a Self-organizing Map Clustering Algorithm Considering a Multi-Scenario Technique. Energies 2017, 10, 2144. [Google Scholar] [CrossRef]
  24. Liu, B.; Zhao, R.; Wang, G. Uncertain Planning and Application, 1st ed.; Tsinghua University Press: Beijing, China, 2003; pp. 185–187. [Google Scholar]
  25. Yang, X.S. A New Metaheuristic Bat-Inspired Algorithm. Comput. Knowl. Technol. 2010, 284, 65–74. [Google Scholar]
  26. Khan, K.; Sahai, A. A Comparison of BA, GA, PSO, BP and LM for Training Feed forward Neural Networks in e-Learning Context. Int. J. Intell. Syst. Appl. 2012, 4, 23–29. [Google Scholar] [CrossRef]
  27. Ghosh, S.; Kaur, M.; Bhullar, S.; Karar, V. Hybrid ABC-BAT for Solving Short-Term Hydrothermal Scheduling Problems. Energies 2019, 12, 551. [Google Scholar] [CrossRef]
  28. Dinh, B.H.; Nguyen, T.T.; Quynh, N.V.; Van Dai, L. A Novel Method for Economic Dispatch of Combined Heat and Power Generation. Energies 2018, 11, 3113. [Google Scholar] [CrossRef]
  29. Fayaz, M.; Kim, D. Energy Consumption Optimization and User Comfort Management in Residential Buildings Using a Bat Algorithm and Fuzzy Logic. Energies 2018, 11, 161. [Google Scholar] [CrossRef]
  30. Jose, J.T. In Economic load dispatch including wind power using Bat Algorithm. In Proceedings of the 2014 International Conference on Advances in Electrical Engineering (ICAEE), Vellore, India, 9–11 January 2014. [Google Scholar]
  31. Wang, J.Q.; Chen, J.; Qu, T.; Huang, G.Q.; Zhang, Y.F.; Sun, S.D. New entropy weight-based TOPSIS for evaluation of multi-objective job-shop scheduling solutions. In Proceedings of the 2012 IEEE International Conference on Industrial Engineering and Engineering Management, Hong Kong, China, 10–13 December 2012. [Google Scholar]
  32. Lu, J.; Wang, W.; Zhang, Y.; Cheng, S. Multi-Objective Optimal Design of Stand-Alone Hybrid Energy System Using Entropy Weight Method Based on HOMER. Energies 2017, 10, 1664. [Google Scholar] [CrossRef]
  33. Sun, C.; Mi, Z.; Ren, H.; Jing, Z.; Lu, J.; Watts, D. Multi-Dimensional Indexes for the Sustainability Evaluation of an Active Distribution Network. Energies 2019, 12, 369. [Google Scholar] [CrossRef]
  34. Li, H.; Li, G.; Liu, S.; Wang, Y.; Wang, Z.; Wang, J.; Zhang, N. Research on Optimal Planning of Access Location and Access Capacity of Large-Scale Integrated Wind Power Plants. Energies 2017, 10, 442. [Google Scholar] [CrossRef]
  35. Huang, W.; Xiong, W.; Yan, B.; Liu, L.F.; Liu, Z.F. Optimization scheduling strategy of virtual microgrid with at different time scales. Automat. Electr. Power Syst. 2017, 41, 12–19. [Google Scholar]
  36. Zhang, X.H.; Zhao, C.M.; Liang, J.X.; Li, K.; Zhong, J.Q. Robust Fuzzy Economic Scheduling for Power Systems Considering Bilateral Uncertainty of Power Generation. Automat. Electr. Power Syst. 2018, 17, 67–75. [Google Scholar]
Figure 1. IEEE33 node power distribution system.
Figure 1. IEEE33 node power distribution system.
Energies 12 01711 g001
Figure 2. The power of load, wind turbine, photovoltaic output and purchased power by microgrid.
Figure 2. The power of load, wind turbine, photovoltaic output and purchased power by microgrid.
Energies 12 01711 g002
Figure 3. Outputs of the intraday phase from 11:00 to 17:00: (a) Interruption of the interruptible load; (b) unplanned purchases; (c) PLR.
Figure 3. Outputs of the intraday phase from 11:00 to 17:00: (a) Interruption of the interruptible load; (b) unplanned purchases; (c) PLR.
Energies 12 01711 g003
Figure 4. The output in each intraday credibility scenario: (a) Γ = 1 ; (b) Γ = 3 ; (c) Γ = 5 .
Figure 4. The output in each intraday credibility scenario: (a) Γ = 1 ; (b) Γ = 3 ; (c) Γ = 5 .
Energies 12 01711 g004aEnergies 12 01711 g004b
Figure 5. The output in each day-ahead credibility scenario: (a) Γ = 1 ; (b) Γ = 3 ; (c) Γ = 5 .
Figure 5. The output in each day-ahead credibility scenario: (a) Γ = 1 ; (b) Γ = 3 ; (c) Γ = 5 .
Energies 12 01711 g005
Figure 6. Prices of different weight coefficients.
Figure 6. Prices of different weight coefficients.
Energies 12 01711 g006
Figure 7. Different objective function values when different weights are applied.
Figure 7. Different objective function values when different weights are applied.
Energies 12 01711 g007
Figure 8. Convergence curves of two algorithms.
Figure 8. Convergence curves of two algorithms.
Energies 12 01711 g008
Table 1. Time-of-use price period price.
Table 1. Time-of-use price period price.
PeriodTimePrice (yuan/(kW·h))
Peak8:00–12:00, 16:00–20:000.55
Flat6:00–8:00, 12:00–16:00, 20:00–22:000.52
Valley22:00–24:00, 0:00–6:000.3
Table 2. Operating costs for each scenario.
Table 2. Operating costs for each scenario.
ScenarioDay-ahead Cost (yuan)The Cost of IL (yuan)Extra Purchase Cost (yuan)Peak Load Regulation (yuan)Total Cost (yuan)
128,692953.47876.211844.632,264.28
228,9551104.38659.791305.9332,025.1
330,845894.33326.851178.4733,244.65
429,681917.4235.42904.8131,738.63
Table 3. Entropy weight of objective functions.
Table 3. Entropy weight of objective functions.
Entropy WeightsData
w10.3145
w20.6855
Table 4. Comparison of two algorithms.
Table 4. Comparison of two algorithms.
AlgorithmPSOBA
Optimal solution0.93440.9632
Standard deviation0.02790.0001
Average time consuming (s)18.4724.58
Operating costs (yuan)2802227715
User transfer coefficient0.07630.0759

Share and Cite

MDPI and ACS Style

Sheng, S.; Gu, Q. A Day-ahead and Day-in Decision Model Considering the Uncertainty of Multiple Kinds of Demand Response. Energies 2019, 12, 1711. https://doi.org/10.3390/en12091711

AMA Style

Sheng S, Gu Q. A Day-ahead and Day-in Decision Model Considering the Uncertainty of Multiple Kinds of Demand Response. Energies. 2019; 12(9):1711. https://doi.org/10.3390/en12091711

Chicago/Turabian Style

Sheng, Siqing, and Qing Gu. 2019. "A Day-ahead and Day-in Decision Model Considering the Uncertainty of Multiple Kinds of Demand Response" Energies 12, no. 9: 1711. https://doi.org/10.3390/en12091711

APA Style

Sheng, S., & Gu, Q. (2019). A Day-ahead and Day-in Decision Model Considering the Uncertainty of Multiple Kinds of Demand Response. Energies, 12(9), 1711. https://doi.org/10.3390/en12091711

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop