Modeling and Investigation of Demand Response Uncertainty on Reliability Assessment

Abstract

Demand Response (DR) provides an opportunity for customers to reduce their loads during times of high prices and therefore to shave the peak loads. The power outputs of large-scale generator units can be predicted and controlled easily; therefore, the pricing and reliability of conventional power utilities can be assessed straightforwardly. However, the customer loads are very variable and difficult to predict and control; therefore, the integration of DR might cause uncertainty issues on pricing and reliability and is essential to be further investigated. A novel uncertainty model for load reduction is proposed in this paper. The probability intensities of load reduction are first estimated from the measured load reduction variations. A multi-state probability model is then proposed for load reduction and the Markov process is used to calculate the state probabilities. A stochastic analysis scheme using Monte Carlo simulation for pricing and reliability taking the DR uncertainty into account is then investigated. Several cases are designed to compare the effects of DR uncertainty. Simulation results show that the proposed uncertainty model can be integrated into conventional economic dispatch to precisely evaluate the DR uncertainty on system operation and reliability.

1. Introduction

In conventional power systems, customer loads can respond to predetermined electricity prices and then adjust their load profiles to obtain economic benefits. The procedure was treated as a Demand Side Management (DSM). Due to deregulation and energy efficiency, power utilities confronted many new challenges such as suppressing the peak load, increasing the off-peak load, reaching the regional power balance, improving the safety and reliability of power system operation, etc. Demand Response (DR) provides an opportunity for customers to reduce their loads, coupled with electricity, heat, and other energy, during time slots of high prices and therefore to shave the peak loads. With DR, customers can dynamically adjust their power consumption in response to real-time electricity prices and transfer the power consumption into other time periods. The benefits of DR at least include reducing the peak load, decreasing the number and costs of start-up/shut-down of thermal power plants, enhancing the usage of unstable renewable energy, defer the investments of new power plants and transmission lines, improving the reliability and security of system operation and so on. Many papers have been published to investigate the definitions, classifications, architectures, advantages, benefits, costs, etc. of different DR programs. For example, ref. [1] investigated the reliability of DR resources. The potential advantages of DR were also discussed. Ref. [2] represented the achievements and deployments of DR and advanced metering in the United States. The definition and classification of DR and the overview of DR in the electricity market were presented in [3]. Ref. [4] reviewed and analyzed the trend of DR from six different topics and the characterization of the future DR was obtained. A new performance-based method for characterizing and assessing the resilience of multi-functional DR with the integrated energy system was proposed in [5]. Energy management and planning in smart cities integrating DR was proposed in [6]. Base on analysis and collection of DR and load scheduling, the possible realization form of load scheduling in the future was designed in [7]. Taking the traditional DR in a power system as a starting point, the fundamental theory, framework design, and potential estimation of integrated DR in power systems were introduced in [8]. The current theoretical research and application of integrated DR were also assessed. Ref. [9] explored the main industry drivers of smart grid and DR and discussed different DR programs. A comprehensive review of several DR techniques with specific focuses on pricing signals, optimization, appliance scheduling was compared in [10]. The various applications of DR were explored in [11]. The implementation and the potential opportunities for DR through environmental protection were also investigated. The benefits and challenges to each market entity using DR-as-a-service were represented in [12]. A case study was performed using the proposed DR model and the results showed that DR curtailment were the most influencing factors that impact the benefit of using DR-as-a-service. Ref. [13] presented a probabilistic modeling approach to minimize network’s overhead lines aging and maximize its reliability at uncertainties by utilizing DR. Ref. [14] proposed the mathematical models and solution schemes for different DR programs. The benefits of DR were also investigated. DR can be divided into incentive-based and price-based categories [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. Demand-side bidding is one of the main schemes of an incentive-based category that individual customers, including residential, commercial, and industrial users, can be aggregated to reduce their power consumption during a settled economic transaction and share the benefits. The price-based category includes the schemes of real-time price, critical-peak price, time-of-use price, etc. that customers reduce their power consumption with respect to the formally predetermined price variations. Hybrid DR are composed of incentive-based and price-based categories.

Some research focused on the technologies and challenges of integrating DR into commercial and industrial customers [15,16]. The integration of renewable generation resources, providing more flexibility for DR, were studied in [17,18]. An economic model of DR interacting with renewable resources in markets was proposed in [17]. A dynamic market mechanism that reaches the market equilibrium through continuous negotiations between key market players and renewables was represented in [18]. DR can achieve many advantages and make energy usage more effective and efficient; however, a higher level of uncertainty in customer load behaviors should be investigated. The characteristics of renewable generation resources are quite different from conventional generation units, they also cause more uncertainty in DR; therefore, it is vital to study the pricing and reliability of DR integrating with renewable generation resources. Only a few works paid attention to DR with uncertainty [19,20,21,22,23,24]. Ref. [19] discussed the stochastic unit commitment with uncertain DR. A stochastic DR representation was studied by scenarios and each scenario corresponded to a price-elastic demand curve. The loss of load probability was then calculated and compared. Ref. [20] proposed a computational framework to take the clearing process uncertainties in the day-ahead electricity market with DR providers into account. An economic dispatch based on the scenario approach was also represented. The worst-case conditional value-at-risk theory with minimizing the worst-case cost caused by a disaster for a community energy system considering DR uncertainty was proposed in [21]. Different risk preferences of natural hazards were also studied. Ref. [21] focused on the community energy system and did not specifically investigate the uncertainty of renewable energy resources in DR. A stochastic hybrid system model to capture DR continuous dynamics and discrete events that arise from failures and repairs was proposed in [22]. The DR probability that can support a certain amount of power during a time period can be estimated by a probability-capacity-duration contour; however, the reliability assessment of DR uncertainty did not study. Refs. [20,21,22] did not fully investigate the DR uncertainty in system operation and reliability especially for the DR with renewable energy resources. Refs. [23,24] studied the uncertainty effects of DR in distribution systems.

The power outputs of large-scale generator units can be predicted and controlled easily; therefore, the pricing and reliability of conventional power utility can be assessed straightforwardly. Many probabilistic evaluation techniques considering failure and repair probabilities can be used to calculate the probability distribution of generation and then be utilized to evaluate system reliability [25,26]. The customer loads are very variable and difficult to predict and control especially the integration of renewable generation resources; therefore, the utilization of DR might cause higher uncertainty on pricing and reliability and is essential to be further investigated. A novel uncertainty model for load reduction is proposed in this paper. The probability intensities of load reduction are first estimated from the measured load reduction variations. A multi-state probability model is then proposed for load reduction and the Markov process is used to calculate the state probabilities. Monte Carlo simulation is used to investigate the DR uncertainty on pricing and reliability. Several cases are designed to demonstrate the validity of the proposed method. The integration of DR can effectively improve the system’s reliability. However, if the DR uncertainty cannot be effectively considered, the results may be too optimistic and may cause reliability issues and additional interruption costs. The proposed method can effectively integrate the DR uncertainty into account and accurately evaluate the system’s reliability.

2. Basic Concepts and Mathematical Models of DR

The cost of supplying system load without DR is the solution of economic dispatch based on the marginal costs of generator units. The cost function of generator unit i is usually assumed to be a quadratic polynomial and can be expressed as [14]

$$G{C}_i(P_G^i)=\frac{1}{2}{a}_i^{}{(P_G^i)}^2+b_i^{}{P}_G^i+c_i^{},$$

(1)

where $G{C}_i(•)$ is the cost function of generator unit i. $P_G^i$ is the power output of generator unit i. $a_i^{}$, $b_i^{}$ and $c_i^{}$ are coefficients of the quadratic polynomial.

The minimum cost for economic dispatch without DR can be obtained by solving the following problem.

$$\begin{array}{l}T{C}_G(P_G^{})=\min {\displaystyle \sum _{i=1}^{N_G}G{C}_i(P_G^i)}\\ s.t.\\ {\displaystyle \sum _{i=1}^{N_G}{P}_G^i}-P_{loss}-P_{SL}=0\\ P_G^{i,\min }\le P_G^i\le P_G^{i,\max }\hspace{1em}i=1\cdots N_G\end{array},$$

(2)

where $T{C}_G(•)$ is the total economic dispatch cost without DR. $N_G$ is the number of interconnected generator units. $P_{SL}$ and $P_{loss}$ are the system load and loss, respectively. $P_G^{i,\min }$ and $P_G^{i,\max }$ are the allowable minimum and maximum power output of generator unit i, respectively. If necessary, the line flow and bus voltage constraints can be integrated into (2) and be treated as a security-constrained economic dispatch.

The cost function of load reduction is more difficult to be assessed; however, it can also be considered as a quadratic function for simplification [14]. Therefore, the cost function for load reduction at bus i can be written as

$$D{C}_i(P_D^i)=\frac{1}{2}{d}_i^{}{(P_D^i)}^2+e_i^{}{P}_D^i+f_i^{},$$

(3)

where $D{C}_i(•)$ is the cost function for load reduction at bus i. $P_D^i$ is the load reduction at bus i. $d_i^{}$, $e_i^{}$, and $f_i^{}$ are coefficients of the quadratic polynomial.

For incentive-based DR such as demand-side bidding, the minimum cost for serving system load with DR can be obtained by solving the following problem.

$$\begin{array}{l}T{C}_{DR}^{IB}(P_G^{},P_D^{})=\min \left({\displaystyle \sum _{i=1}^{N_G}G{C}_i(P_G^i)}+{\displaystyle \sum _{i=1}^{N_D}D{C}_i(P_D^i)}\right)\\ s.t.\\ {\displaystyle \sum _{i=1}^{N_G}{P}_G^i}-{\displaystyle \sum _{i=1}^{N_D}{P}_D^i}-P_{loss}-P_{SL}=0\\ P_G^{i,\min }\le P_G^i\le P_G^{i,\max }\hspace{1em}i=1\cdots N_G\\ P_D^{i,\min }\le P_D^i\le P_D^{i,\max }\hspace{1em}i=1\cdots N_D\end{array},$$

(4)

where $T{C}_{DR}^{IB}(•,•)$ is the total economic dispatch cost with incentive-based DR. $N_D$ is the number of DRs. $P_D^{i,\min }$ and $P_D^{i,\max }$ are the allowable minimum and maximum load reductions at bus i, respectively.

The marginal costs for generation unit and load reduction can be calculated by

$$\begin{array}{l}{\lambda }_G^i=\frac{\partial \left(G{C}_i(P_G^i)\right)}{\partial P_G^i}=a_i^{}{P}_G^i+b_i^{}\\ {\lambda }_D^i=\frac{\partial \left(D{C}_i(P_D^i)\right)}{\partial P_D^i}=d_i^{}{P}_D^i+e_i^{}\end{array},$$

(5)

where ${\lambda }_G^i$ and ${\lambda }_D^i$ are the marginal costs for generation unit i and load reduction at bus i, respectively.

The solution of system marginal cost, power output of generator unit i and load reduction at bus i for demand-side bidding without taking the system loss into account can be calculated by

$$\begin{array}{l}{\lambda }_S^{}=\frac{{\displaystyle \sum _{i=1}^{N_G}\frac{b_i^{}}{a_i}}+{\displaystyle \sum _{i=1}^{N_D}\frac{e_i^{}}{d_i}}+P_{SL}}{{\displaystyle \sum _{i=1}^{N_G}\frac{1}{a_i}+{\displaystyle \sum _{i=1}^{N_D}\frac{1}{d_i}}}}\\ P_G^i=\frac{{\lambda }_S-b_i}{a_i}\hspace{1em}i=1\cdots N_G\\ P_D^i=\frac{{\lambda }_S-e_i}{d_i}\hspace{1em}i=1\cdots N_D\end{array},$$

(6)

where ${\lambda }_S^{}$ is the system marginal cost.

For the price-based DR, the price for load reduction should be determined by Independent System Operator (ISO) first, the objective can be expressed as

$$\begin{array}{l}T{C}_{DR}^{PB}(P_G^{},P{R}_{DR}^{opt})=\min \left({\displaystyle \sum _{i=1}^{N_G}G{C}_i(P_G^i)}+{\displaystyle \sum _{i=1}^{N_D}P{R}_{DR}^{opt}\ast P_D^i}\right)\\ s.t.\\ {\displaystyle \sum _{i=1}^{N_G}{P}_G^i}-{\displaystyle \sum _{i=1}^{N_D}{P}_D^i}-P_{loss}-P_{SL}=0\\ P_G^{i,\min }\le P_G^i\le P_G^{i,\max }\hspace{1em}i=1\cdots N_G\\ P_D^{i,\min }\le P_D^i\le P_D^{i,\max }\hspace{1em}i=1\cdots N_D\end{array},$$

(7)

where $T{C}_{DR}^{PB}(•,•)$ is the total economic dispatch cost with priced-based DR. $P{R}_{DR}^{opt}$ is the optimal price for price-based DR determined by ISO.

The solution of the optimal price for price-based DR, system marginal cost, power output of generator unit i and load reduction at bus i without taking the system loss into account can be calculated by

$$\begin{array}{l}P{R}_{DR}^{opt}=\frac{{\lambda }_S({\displaystyle \sum _{i=1}^{N_D}\frac{1}{d_i}})+{\displaystyle \sum _{i=1}^{N_D}\frac{e_i^{}}{d_i}}}{{\displaystyle \sum _{i=1}^{N_D}\frac{1}{d_i}}}\\ {\lambda }_S^{}=\frac{{\displaystyle \sum _{i=1}^{N_G}\frac{b_i^{}}{a_i}}+\frac{1}{2}{\displaystyle \sum _{i=1}^{N_D}\frac{e_i^{}}{d_i}}+P_{SL}}{{\displaystyle \sum _{i=1}^{N_G}\frac{1}{a_i}+\frac{1}{2}{\displaystyle \sum _{i=1}^{N_D}\frac{1}{d_i}}}}\\ P_G^i=\frac{{\lambda }_S-b_i}{a_i}\hspace{1em}i=1\cdots N_G\\ P_D^i=\frac{P_{DR}^{opt}-e_i}{d_i}\hspace{1em}i=1\cdots N_D\end{array},$$

(8)

The benefits of incentive-based and price-based DRs can be expressed by

$$\begin{array}{l}B{C}_{DR}^{IB}= T{C}_{DR}^{IB}(P_G^{},P_D^{})-T{C}_G(P_G^{})\\ B{C}_{DR}^{PB}= T{C}_{DR}^{PB}(P_G^{},P{R}_{DR}^{opt})-T{C}_G(P_G^{})\end{array},$$

(9)

where $B{C}_{DR}^{IB}$ and $B{C}_{DR}^{PB}$ are the benefits of incentive-based and price-based DRs, respectively.

Note that the coefficients of quadratic functions for incentive-based and price-based DRs might be different and can be estimated and predicted from market information and/or customer interruption costs. In the following, only incentive-based DR is discussed and simulated due to limited space.

3. Proposed Uncertainty Model for DR

The outputs of large-scale utility generators can be predicted and controlled more effortlessly; therefore, the pricing and reliability in conventional power utilities can be easily assessed. Traditionally, a two-state probability model as shown in Figure 1 is used for most of the reliability assessments where ${\lambda }_G^i$ and ${\mu }_G^i$ are the failure rate and repair rate of generator unit i, respectively. The Mean Time To Failure (MTTF) and Mean Time To Repair (MTTR) can be expressed as

$$\begin{array}{l}MTT{F}_G^i=\frac{1}{{\lambda }_G^i}\\ MTT{R}_G^i=\frac{1}{{\mu }_G^i}\end{array},$$

(10)

where $MTT{F}_G^i$ and $MTT{R}_G^i$ are the MTTF and MTTR of generator unit i, respectively.

The availability and unavailability of generator unit i can be expressed as

$$\begin{array}{l}{A}_G^i=\frac{{\mu }_G^i}{{\mu }_G^i+{\lambda }_G^i}\\ U{A}_G^i=1-A_G^i\end{array},$$

(11)

where $A_G^i$ and $U{A}_G^i$ are the availability and unavailability of generator unit i, respectively.

The probability density function for a two-state model of generator unit i is shown in Figure 2. After those parameters as illustrated in (10) and (11) were obtained, the reliability indices of a power utility can be assessed effortlessly by capacity outage table with convolution, effective load approach, Monte-Carlo simulation, etc. [25,26].

The customer loads are very variable and difficult to predict and control; therefore, DR especially the integration of renewable generation resources might cause higher uncertainty on pricing and reliability. A novel multi-state probability model based on the real measurements for load reduction is proposed in this paper. Due to the higher variation of load reduction and renewable generation resources, the transitions between different states of load reduction as illustrated in Figure 3 might occur frequently. Generally, the DR contract is confirmed if the measured power is larger than the set power and is denoted as “Normal State” in Figure 3. The DR contract is failed and the customer will be charged a penalty when the measured power is lower than a certain ratio of the set power and is denoted as “Fail State” in Figure 3. Due to the load reduction variations, several derated load reduction states denoted as “Derated State” may be occurred. ${\lambda }_{ij}^i$ in Figure 3 is the transition intensity from state i to state j. If the transition intensities between the different states can be estimated, the probability of each state can be solved by a Markov-based state-space method.

The definition of transition intensity between the states can be expressed as

$${\lambda }_{ij}=\underset{\Delta t\to 0}{\lim }\frac{1}{\Delta t}p[x(t+\Delta t)=j|x(t)=i)],$$

(12)

where $x(t)$ is the random variable of system state at time t. $p(•)$ is a probability function.

In most power system applications, the transition intensities don’t depend on the time t; therefore, the transition intensity can be rewritten as

$${\lambda }_{ij}\Delta t=p[x(t+\Delta t)=j|x(t)=i)],$$

(13)

From (13), it can be seen that if the conditional probability between different states can be calculated, then the transition intensity can be calculated accordingly. Figure 4 shows the load reduction of a customer with wind generation resources. In Figure 4, the measured powers are actual measurements and the set DR powers are the assumed values. From Figure 3, the variation of load reduction can be significantly observed; therefore, the probabilities between the measured and set load reduction values can be identified. If the uncertainty of load reduction can be quantified, the pricing and reliability of DR can be effectively analyzed. Using the actual measurements and time sequences, the transitions between different states can be counted, and then the transition intensities can be calculated. Using the enlarged part of Figure 4 as an example and shown in Figure 5, the set power at each time period is treated as 100%. In this example, the normal state is the measured power larger than 90% of the set power; the derated states 1–3 indicate that the measured power is in 90% to 80%, 80% to 70%, and 70% to 60% of the set power, respectively; and the fail state 4 means that the measured power is less than 60% of the set power. Based on the time-sequence process of Figure 5, the percentages and state transitions of load reduction are listed in

Table 1. Percentages and state transitions of load reduction in Figure 5.

Measured Point	Power Reduction (%)	From State	To State	Measured Point	Power Reduction (%)	From State	To State
1	100	0	0	11	100	1	0
2	100	0	0	12	100	0	0
3	98	0	0	13	51	0	4
4	89	0	1	14	98	4	0
5	100	1	0	15	56	0	4
6	100	0	0	16	56	4	4
7	100	0	0	17	22	4	4
8	100	0	0	18	36	4	4
9	100	0	0	19	78	4	2
10	85	0	1	20	100	2	0

. The number of state transitions is shown in Figure 6. From Table 1 and Figure 6, it can be observed that there are 8 transitions from state 0 to state 0, 2 transitions from state 0 to state 1, 2 transitions from state 0 to 4, etc. Therefore, the number of state transitions for Figure 6 can be expressed in matrix form as

$$N_{ST}=\left[\begin{array}{ccccc}8& 2& 1& 0& 2\\ 2& 0& 0& 0& 0\\ 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 1& 1& 1& 0& 3\end{array}\right],$$

(14)

where $N_{ST}$ is the number of state transition matrix.

Using 4320 actual measurements as an example, the number of state transitions can be calculated similarly as shown in Table 1 and Figure 6 and be written as

$$N_{ST}=\left[\begin{array}{ccccc}3043& 202& 135& 79& 44\\ 146& 79& 25& 30& 22\\ 107& 23& 66& 28& 10\\ 59& 27& 18& 46& 22\\ 30& 13& 20& 15& 31\end{array}\right],$$

(15)

The numbers of transitions started from states 0 to 4 are 3503, 302, 234, 172, and 109, respectively. According to (12) and (13), the transition intensity can be calculated by dividing the number of transitions by the number of transitions started from the corresponding state. For example, the transition intensity from state 0 to 1 can be calculated by dividing 202 by 3503. Therefore, the transition intensity matrix can be expressed as

$${\lambda }_{ij}=\left[\begin{array}{ccccc}0.8687& 0.0577& 0.0385& 0.0226& 0.0126\\ 0.4834& 0.2616& 0.0828& 0.0993& 0.0728\\ 0.4573& 0.0983& 0.2821& 0.1197& 0.0427\\ 0.3430& 0.1570& 0.1047& 0.2674& 0.1279\\ 0.2752& 0.1193& 0.1835& 0.1376& 0.2844\end{array}\right],$$

(16)

where ${\lambda }_{ij}$ is the transition intensity matrix.

After the transition intensities were estimated from the historical measurement data, a Markov-process-based state space method [25,26] as expressed in (17) can be used to calculate the state probabilities. Due to limited space, the detail derivations of the Markov-process-based state space method are not shown here.

$$\{\begin{array}{c}{\left({\lambda }_{ij}-I\right)}^{T}P{R}_S=0\\ {\displaystyle \sum _{i=1}^{N_S}P{R}_S[i]}=1\end{array},$$

(17)

where $I$ is an identity matrix. $P{R}_S$ is the vector of state probability. $N_S$ is the number of states.

Using (16) as an example, the state probabilities can be calculated by

$$\begin{array}{l}\left[\begin{array}{ccccc}-0.1313& 0.4834& 0.4573& 0.3420& 0.2752\\ 0.0577& -0.7384& 0.0983& 0.1570& 0.1193\\ 0.0385& 0.0828& -0.7179& 0.1047& 0.1835\\ 0.0226& 0.0993& 0.1197& -0.7326& 0.1376\\ 1& 1& 1& 1& 1\end{array}\right]\left[\begin{array}{c}P{R}_S[0]\\ P{R}_S[1]\\ P{R}_S[2]\\ P{R}_S[3]\\ P{R}_S[4]\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ 1\end{array}\right]\\ \left[\begin{array}{c}P{R}_S[0]\\ P{R}_S[1]\\ P{R}_S[2]\\ P{R}_S[3]\\ P{R}_S[4]\end{array}\right]=\left[\begin{array}{c}0.7595\\ 0.0852\\ 0.0673\\ 0.0526\\ 0.0354\end{array}\right]\end{array},$$

(18)

The probability density function of load reduction under different states is shown in Figure 7. With the actual measurements, a multi-state probability model for load reduction is proposed in this paper. Using the multi-state model, the load reduction uncertainty on DR can be investigated effectively. The detailed steps of the proposed multi-state probability model are described below:

(1)

Collect measured power and set power data and calculate the load reduction percentage of each measurement as listed in Table 1;

(2)

Determine the number of states for load reduction (e.g., 5 states are used in the above example);

(3)

Use the time-sequence process and load reduction percentage to count the number of transitions between different states as the example shown in Table 1 and Figure 6;

(4)

Build the number of state transition matrix as expressed in (15) and calculate the number of transitions started from each state;

(5)

Build the transition intensity matrix by dividing the number of transitions by the number of transitions started from the corresponding state as expressed in (16);

(6)

Use a Markov-process-based state space method as expressed in (17) to calculate the state probabilities. After the state probabilities were calculated, a probability density function of the multi-state model for a load reduction as illustrated in Figure 7 can be obtained.

4. Investigations of DR Uncertainty on Operation and Reliability

The concepts and mathematical models of DR as described in Section 2 did not consider the uncertainties of generator unit and customer load reduction. With the proposed multi-state probability model of load reduction, the load reduction uncertainty on DR can be investigated effectively and precisely. The cumulative probability density functions of the generator units and customer load reductions can be obtained by the formulas proposed in Section 3 and then the Monte Carlo simulation can be used to calculate the expected total economic dispatch cost and Loss-of Load Probability (LOLP) straightforwardly. The expected total economic dispatch cost including the customer interruption cost of unsupplied load can be estimated by

$$T{C}_S^E=E(T{C}_{DR}^{IB}(P_G^{},P_D^{})+C_{CIC}(P_G^{},P_D^{})),$$

(19)

where $T{C}_S^E$ is the expected total economic dispatch cost including customer interruption cost considering the uncertainties of generator unit and load reduction. $E(•)$ is expected function. $C_{CIC}$ is the customer interruption cost of unsupplied system load.

The $C_{CIC}(•,•)$ can be calculated by (20).

$$C_{CIC}(P_G^{},P_D^{})=P_{SL}^{UNS}\ast M{W}_{CIC},$$

(20)

where $P_{SL}^{UNS}$ is the unsupplied system load and can be calculated by $P_G^{}$ and $P_D^{}$. $M{W}_{CIC}$ is customer interruption cost per MW ($/MW).

Equation (19) cannot be calculated by a deterministic method; therefore, a stochastic method, Monte-Carlo simulation, is used in this paper to calculate the cost of (19). The procedures are as follows:

(1)

Establish the two-state probability model of each generator unit and the multi-state probability model of each load reduction as proposed in Section 3;

(2)

Establish the cumulative probability density functions for each generator unit and load reduction;

(3)

Determine the number of Monte-Carlo simulation;

(4)

For each simulation, generate the random variables for each generator unit and load reduction and determine the availability or unavailability of each generator unit and actual load reduction of each DR customer;

(5)

Solve the problem of incentive-based DR as expressed in (4);

(6)

Calculate the total generation output, total load reduction, unsupplied system loads, and so on for this simulation;

(7)

Repeat (4) to (6) until the number of Monte-Carlo simulations is completed.

(8)

Calculate (19) and the statistical data including average and standard deviation of system marginal price, total economic dispatch cost, total generation output, total load reduction, unsupplied system load, etc.

5. Test Results and Discussions

A 6-bus system as shown in Figure 8 is simulated in this paper. From Figure 8, it can be seen that buses 1–3 are generator buses and buses 4–6 are demand buses. The parameters of cost function and the lower and upper limits of each generator bus are shown in

Table 2. Parameters of cost functions for the test system.

Bus	a	b	c	Min (MW)	Max (MW)
1	0.1633	11.669	213.1	50	180
2	0.1689	10.333	200.0	40	120
3	0.1441	10.833	240.0	30	125
	d	e	f	Min (MW)	Max (MW)
4	0.80	25.5	167.5	0	30
5	-	-	-	-	-
6	0.85	25.5	173.8	0	30

. The generator units installed at each generator bus and the reliability parameters of the two-state probability model are shown in

Table 3. Reliability parameters of generation units.

Bus	Capacity (MW)	Type	MTTF (h)	MTTR (h)	UA *
1	36	Thermal	1460	45	0.03
	36	Thermal	1460	45	0.03
	36	Thermal	1460	45	0.03
	36	Thermal	1460	45	0.03
	36	Thermal	1460	45	0.03
2	30	Hydro	2920	60	0.02
	30	Hydro	2920	60	0.02
	30	Hydro	2920	60	0.02
	30	Hydro	2920	60	0.02
3	25	Thermal	1752	45	0.025
	25	Thermal	1752	45	0.025
	25	Thermal	1752	45	0.025
	25	Thermal	1752	45	0.025
	25	Thermal	1752	45	0.025

‘*’: UA is unavailability.

. From Table 3, it can be seen that there are 5, 4, and 5 generator units installed at buses 1, 2, and 3, respectively. The loads at buses 4, 5, and 6 are 120, 140, and 140 MW, respectively, therefore, the system load is 400 MW. Buses 4 and 6 participate in DR and the parameters of the cost function and the lower and upper limits of each load reduction are also shown in Table 2. The parameters of Table 2 were obtained from Ref. [14] with minor modification for test purposes. The reliability parameters of multiple-state probability models for load reductions of buses 4 and 6 are listed in

Table 4. Reliability parameters of load reduction buses.

Bus	PRS[State]
	PRS[0]	PRS[1]	PRS[2]	PRS[3]	PRS[4]
4	0.7579	0.0852	0.0673	0.0526	0.0354
6	0.8252	0.0756	0.0461	0.0308	0.0223

where the states 0–4 indicate the measured power is larger 90%, in 90% to 80%, in 80% to 70%, in 70% to 60% and less than 60% of the set power, respectively.

The number of Monte-Carlo simulation is set as 50,000 and the customer interruption cost is 140 $/MW. Cases as listed below are simulated and discussed:

• Case A: Economic dispatch without DR and uncertainty is neglected;
• Case B: Economic dispatch with DR, but the uncertainty is neglected;
• Case C: Economic dispatch with two-state generator model but without DR;
• Case D: Economic dispatch with two-state generator model and incentive-based DR, but the DR uncertainty is neglected;
• Case E: Economic dispatch with two-state generator model and incentive-based DR and the proposed DR uncertainty is considered.

For each case, the total economic dispatch cost, total customer interruption cost, LOLP, etc. are calculated. For Case A, the marginal cost calculated by (2) is 68.63 $/MW, and the power outputs for generating buses 1, 2, and 3 are 174.43 MW, 120 MW, and 125 MW, respectively. The total economic dispatch cost without DR is 14,934.65$. For Case B, the marginal cost calculated by (4) is 54.21 $/MW, and the power outputs for generating buses 1, 2, and 3 are 130.27 MW, 120 MW, and 125 MW, respectively. The load reductions for buses 4 and 6 are 20.60 MW and 19.60 MW, respectively. The total economic dispatch cost with DR is 13,912.86$. Obviously, the total economic dispatch cost and marginal cost are significantly reduced due to the integration of DR.

The cumulative failure probability density function of generator units with a two-state probability model is shown in Figure 9. From Figure 9, it can be observed the failure probability under the system load of 400 MW is 0.2383 and therefore LOLP is about 0.2383. It is a little high due to the lower reserve without DR. For Case C, the expected marginal cost is 66.34 $/MW with a standard deviation of 5.43 $/MW and the expected total generation cost is 15,142.51$ with a standard deviation of 361.58$. The expected unsupplied system load is 14.51 MW with a standard deviation of 20.53 MW and therefore the expected total customer interruption cost is 2031.99$ with a standard deviation of 2873.93 $. The expected total economic dispatch cost including customer interruption cost is 17,174.51$ with a standard deviation of 3142.39$. Comparing Cases A and C, it can be seen that the generator failure uncertainty would cause the higher total economic dispatch cost. For Case D, the failure probability is 0.0078 and therefore LOLP is about 0.0078 since the system load can be regarded as 400 MW minus the maximum load reduction of 60 MW. The expected marginal cost is 56.39 $/MW with a standard deviation of 3.70 $/MW and the expected total generation cost is 14,062.68$ with a standard deviation of 280.03$. The expected unsupplied system load is 0.59 MW with a standard deviation of 4.22 MW and therefore the expected total customer interruption cost is 41.37$ with a standard deviation of 295.34$. The expected total economic dispatch cost including customer interruption cost is 14,145.42$ with a standard deviation of 753.84$. Comparing Cases B and D, it can be observed that total economic dispatch cost will only be slightly increased if the load reduction is absolutely reliable.

The cumulative failure probability density function of load reductions with the proposed multi-state model is shown in Figure 10. The cumulative failure probability density function of generation units and load reductions is shown in Figure 11. From Figure 11, it can be observed the failure probability under the system load of 400 MW is 0.01244 and therefore LOLP is about 0.01244. For Case E, the expected marginal cost is 56.73 $/MW with a standard deviation of 3.82 $/MW and the expected total generation cost is 14,082.74 $ with a standard deviation of 294.65$. The expected unsupplied system load is 1.01 MW with a standard deviation of 5.50 MW and therefore the expected total customer interruption cost is 141.85$ with the standard deviation of 770.47$. The expected total economic dispatch cost including customer interruption cost is 14,224.59$ with a standard deviation of 944.74$. Comparing Cases D and E, it can be seen that the LOLP is changed from 0.0078 to 0.01244 due to the DR uncertainty. The difference of LOLPs with or without taking DR uncertainty into account is not very significant due to the lower system demand. If the system demand is 440 MW, the LOLPs of Cases D and E are 0.0920 and 0.1667 from Figure 9 and Figure 11, respectively. Figure 12 illustrates the LOLPs while system load is increased from 400 MW to 480 MW. Apparently, from Figure 12, it can be seen that the consideration of DR uncertainty can effectively and accurately evaluate the system reliability and analyze the effects of DR for economic dispatch.

A larger system is also simulated to further verify the validity of the proposed method. Due to the limited space, the data as shown in Table 3 and Table 4 is used and only reliability assessment is discussed. The larger system consists of 36 thermal generator units with a rated capacity of 36 MW, 26 hydro generator units with a rated capacity of 30 MW and 37 thermal generator units with a rated capacity of 25 MW. The total generation capacity and demand reductions are 2857 and 700 MW, respectively. The cumulative failure probability density function of generator units with a two-state probability model for the larger system is shown in Figure 13. The cumulative failure probability density function of load reductions with the proposed multi-state model is shown in Figure 14. The cumulative failure probability density function of generation units and load reductions is shown in Figure 15. Figure 16 illustrates the LOLPs while system load is increased from 2400 MW to 3400 MW. Once again it can be seen from Figure 16 that DR can effectively enhance the system’s reliability. For example, the LOLPs for a total load of 3000 MW are 1.00, 0.01314, and 0.026 for Cases C, D, and E, respectively. The integration of DR such as the results of Cases D and E can effectively improve the system reliability. However, if the DR uncertainty cannot be effectively considered, the results as shown in Case D may be too optimistic and may cause unexpected load shedding. The proposed method can effectively integrate the DR uncertainty into account and accurately evaluate the system’s reliability. The proposed method will be very useful for ISOs and/or power utilities since the DR uncertainty can be fully considered.

6. Conclusions

The customer loads are very variable and difficult to predict and control; therefore, the integration of DR would cause uncertainties on pricing and reliability and is essential to be further investigated. DR will definitely improve system reliability; however, its effects should be carefully evaluated. A novel multi-state probability model was proposed for load reduction based on actual measurements and then the Markov process was used to calculate the state probabilities. A stochastic analysis scheme based on Monte Carlo simulation for pricing and reliability taking the DR uncertainty into account was investigated. Based on the proposed method, the DR uncertainty on the system operation, pricing and reliability can be assessed precisely and effectively. Simulated results showed that if the DR uncertainty cannot be effectively considered, the results may be too optimistic and may cause unexpected reliability issues. The proposed uncertainty model of DR can be integrated into a conventional economic dispatch program to precisely analyze the system reliability and will be helpful for ISOs and power utilities.