Patent Publication Number: US-8983811-B2

Title: Simulating customer behavior for demand response

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This application claims the benefit of U.S. Provisional Application No. 61/553,364, filed on Oct. 31, 2011, the disclosure which is incorporated herein by reference in its entirety. 
    
    
     BACKGROUND OF THE INVENTION 
     The invention relates generally to Demand Response. More specifically, the invention relates to a method that generates load patterns from historical customer consumption data, and generates customer loads from the generated load patterns and price responses from given economic parameters. 
     Demand Response (DR) has become an important topic in recent years. Soaring fuel prices coupled with the environmental need to reduce fuel consumption has directed the players in the electricity sector to search for new ways to manage the generation, transmission and distribution process more efficiently. With the introduction of smart power metering systems to a mass number of customers, real-time monitoring of customer loads is made possible. 
     In order to employ DR mechanisms, methods such as machine learning algorithms are needed to interpret results. Yet, these methods require real world data to be calibrated, which is not readily available. Moreover, obtaining such data requires serious investment, pilot DR studies and in-depth econometric analyses. 
     What is desired is a method to effectively and efficiently interpret customer response and manage it, so that the overall process from electricity generation to consumption is more desirable for both the utility companies and the consumers. The utilities benefit as they can avoid costly levels of peak electrical load, while consumers enjoy lowered bills as a result of their participation. 
     SUMMARY OF THE INVENTION 
     The inventors have discovered that it would be desirable to provide simulation methods that generate load patterns from historical customer consumption data, and generate customer loads from the generated load patterns and price responses from given economic parameters. The historical customer consumption data and price responses are exogenous variables that are provided to the simulator. The historical customer consumption data may be a sample of real world customer load curves which are typically the power consumption of a customer throughout a day, and tend to change depending on if the day is a weekday, or on a weekend. Embodiments use a two-step method that first simulates the power demand of a customer over a time period, and second, generates resulting load levels. Embodiments offer a realistic substitute for a real world system. 
     Embodiments fulfill this need for data during the development stages of such methods. A second benefit is policy evaluation. Embodiments may be used to assess the results of DR policies before they are tested in the real world. It provides valuable data which can be used to develop better pricing strategies for utilities. 
     One aspect of the invention provides a customer electrical load curve simulation that uses historical customer load patterns to generate customer load curves. Methods according to this aspect of the invention include inputting pattern parameters {circumflex over (μ)} 0 , {circumflex over (σ)} 0 , {circumflex over (μ)}=[{circumflex over (μ)} 0+Δt , . . . {circumflex over (μ)} 24 ] and {circumflex over (σ)}=[{circumflex over (σ)} 0+Δt , . . . {circumflex over (σ)} 24 ], and simulation parameters P {tilde over (t)}     start   , n sim  and [{tilde over (t)} start ,{tilde over (t)} end ,Δ{tilde over (t)}] wherein P {tilde over (t)}     start    is the desired starting load level, n sim  is the number of simulated load curves that will be generated, {tilde over (t)} start  is the desired starting time of the simulation, {tilde over (t)} end  is the desired ending time of the simulation and Δ{tilde over (t)} is the desired time resolution for the simulation, {circumflex over (μ)} 0  and {circumflex over (σ)} 0  are the mean and standard deviation of the starting load level of estimated pattern parameters and {circumflex over (μ)} and {circumflex over (σ)} are vectors which represent the mean and the standard deviation of the load changes between each load sample Δt, determining a scaling parameter α=P {tilde over (t)}     start   /μ {tilde over (t)}     start   , where μ {tilde over (t)}     start    is a desired starting time, and multiplying the mean {circumflex over (μ)} and variance a {circumflex over (σ)} parameters by the scaling parameter α, generating random numbers that are distributed according to a standard noinial distribution, wherein generating u {tilde over (t)}  that are uniformly distributed in [0,1] and converting the u {tilde over (t)}  to standard normal random numbers by inverse transform sampling Φ −1 (u)), calculating the load P j,{tilde over (t)} =P j,{tilde over (t)}−Δ{tilde over (t)} +μ {tilde over (t)} +Φ −1 (u {tilde over (t)} )σ {tilde over (t)}  for each simulation replication jε{1, . . . , n sim } and time point {tilde over (t)}ε[{tilde over (t)} start ,{tilde over (t)} end ], and calculating a load matrix 
             P   =       [     P     j   ,     t   ~         ]     ∈     R       n   sim     ×     (             t   ~     end     -       t   ~     start         Δ   ⁢           ⁢     t   ~         +   1     )                 
as output.
 
     Another aspect of the invention provides a method that calculates pattern parameter data from historical customer load data. Methods according to this aspect of the invention include inputting historical electrical load data {circumflex over (P)} ti  for each time interval tε[0,24] and observation i=1, 2, . . . , N, calculating the mean {circumflex over (μ)} 0 Σ i=1   N {circumflex over (P)} 0i /N and the standard deviation {circumflex over (σ)} 0 =√{square root over (Σ i=1   N ({circumflex over (P)} 0i −{circumflex over (μ)} 0 ) 2 /(N−1))} of starting load level P 0 , and calculating the mean {circumflex over (μ)} t =Σ i=1   N Δ{circumflex over (P)} ti /N and the standard deviation {circumflex over (σ)} t =√{square root over (Σ i=1   N (Δ{circumflex over (P)} ti −{circumflex over (μ)} t ) 2 /(N−1))} of the load change ΔP t  between time t−Δt and t. 
     Another aspect of the invention provides a customer response simulation that divides customer electrical demand into categories with different price elasticities and then aggregates the results once a price change occurs to simulate the customer&#39;s response to the changing electricity price. Methods according to this aspect of the invention include inputting parameters P 0 , r 0 , r, v, ε, s and f type  wherein P 0  is the corresponding electrical load level, r 0  is the current price of electricity, r is the proposed price of electricity, the vector v=[v 1 , . . . , v n ] shows the proportion of each end use of electricity within the current consumption for each end use k=1, . . . , n, the vector ε=[ε 1 , . . . ε n ] represents the price elasticity of each end use, the vector s=[s 1 , . . . , s n ] which is the coefficient of variation, the ratio of the standard deviation of the response to the mean load level and f type  is the type of the function that is used in the characterization of the customer&#39;s response, examining the function type f type  comprising if the function type f type  is constant elasticity, calculating the customer&#39;s expected load level for end use k as 
               P   k     =       P   0     ⁢         v   k     ⁡     (     r     r   0       )         ε   k               
where P 0 v k  is the customer&#39;s current power consumption for end use k, and if the function type f type  is linear, estimating the demand curve by calculating the slope m k  of the individual demand curve of each end use k at the current price r 0  and load level P 0  as
 
                 m   k     =       r   0         P   0     ⁢     v   k     ⁢     ε   k           ,         
and calculating the customer&#39;s expected load level for each end use k as
 
               P   k     =     max   ⁢     {     0   ,         P   0     ⁢     v   k       +       r   -     r   0         m   k           }             
wherein the load cannot be negative, and generating random numbers that are distributed according to a standard normal distribution, wherein generating u k  that are uniformly distributed in [0,1] and converting the u k  to standard normal random numbers by inverse transform sampling Φ −1 (u k ), generating the customer&#39;s load level (consumption) for end use k, which is a random number that is normally distributed with mean P k  and standard deviation P k s k , and calculating the customer&#39;s final load level which is the sum of the load levels of each individual end use as P=Σ k=1   n P k +Φ −1 (u k )P k s k .
 
     The details of one or more embodiments of the invention are set forth in the accompanying drawings and the description below. Other features, objects, and advantages of the invention will be apparent from the description and drawings, and from the claims. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  is a plot of exemplary load curves that compare load to customer type. 
         FIG. 2  is a plot of exemplary load curves that compare load to time of year. 
         FIG. 3  is a plot of exemplary load curves that compare load to day of the week. 
         FIG. 4  shows an exemplary curve that illustrates a Gaussian random walk. 
         FIGS. 5 ,  5 A,  5 B and  5 C are an exemplary customer load curve simulation method. 
         FIG. 6  is a plot of load simulation curves that compare load to the time of the day. 
         FIG. 7  is a plot of an exemplary load curve that compares load to the price of electricity (rate). 
         FIG. 8  is a plot showing an exemplary aggregation of demand. 
         FIG. 9  is an exemplary customer response to changing electricity price simulation method. 
         FIG. 10  is a plot showing exemplary linear and constant elasticity comparisons of load to the price of electricity (rate). 
         FIG. 11  is a plot showing exemplary response levels with randomization. 
         FIG. 12  is a 3-dimensional plot showing an exemplary deterministic price response comparing load with the time of day and rate. 
         FIG. 13  is a 3-dimensional plot showing an exemplary randomized price response comparing load with the time of day and rate. 
     
    
    
     DETAILED DESCRIPTION 
     Embodiments of the invention will be described with reference to the accompanying drawing figures wherein like numbers represent like elements throughout. Before embodiments of the invention are explained in detail, it is to be understood that the invention is not limited in its application to the details of the examples set forth in the following description or illustrated in the figures. The invention is capable of other embodiments and of being practiced or carried out in a variety of applications and in various ways. Also, it is to be understood that the phraseology and terminology used herein is for the purpose of description and should not be regarded as limiting. The use of “including,” “comprising,” or “having,” and variations thereof herein is meant to encompass the items listed thereafter and equivalents thereof as well as additional items. 
     The terms “connected” and “coupled” are used broadly and encompass both direct and indirect connecting, and coupling. Further, “connected” and “coupled” are not restricted to physical or mechanical connections or couplings. 
     It should be noted that the invention is not limited to any particular software language described or that is implied in the figures. One of ordinary skill in the art will understand that a variety of alternative software languages may be used for implementation of the invention. It should also be understood that some of the components and items are illustrated and described as if they were hardware elements, as is common practice within the art. However, one of ordinary skill in the art, and based on a reading of this detailed description, would understand that, in at least one embodiment, components in the method and system may be implemented in software or hardware. 
     Embodiments of the invention provide methods, system frameworks, and a computer-usable medium storing computer-readable instructions that provide two simulators to be used with Demand Response (DR) management systems. The first simulator generates customer load patterns from observed, historical customer consumption load data, and then uses the generated load patterns to generate customer load curves. The second simulator uses an econometric characterization of customer response to changing electricity prices. It uses the customer load curves output from the first simulator in conjunction with current electricity prices to generate a customer response for proposed electricity prices. Aspects of the methods may be distributed and executed at and on a plurality of computing devices. For example, a PC, a Mac, a tablet computer, a laptop computer, etc., each employing a mouse, touchpad or touchscreen as a pointing device and I/O for data input and exchange. Embodiments may be deployed as software as an application program tangibly embodied on a non-transitory computer readable program storage device. The application code for execution can reside on a plurality of different types of computer readable media known to those skilled in the art. 
     Embodiments provide input to other methods that study DR mechanisms where real data is not available. Moreover, they provide policy evaluation where one can test potential contract parameters and iteratively improve the contract. 
     Embodiments simulate the electricity load curve of a given customer over a period of time. The load curve reflects the customer&#39;s power consumption within the day and has a distinctive shape depending on customer type (commercial, residential, industrial, etc.) ( FIG. 1 ), time of year, which is a proxy for seasonal weather (spring, summer, fall and winter) ( FIG. 2 ) and day type (weekday, weekend) ( FIG. 3 ). Embodiments use these parameters to generate separate load patterns for each customer type, time of year and day type. For the simulation results to be representative, the shape of the generated load patterns should match the distinctive shape of the target tuple—customer type, time of year and day of week. 
     Electrical load is time dependent, as a result of human activity that varies within the day. Embodiments treat load changes within a day as a time series, without relating them to any underlying cause other than those already isolated above before generating the load patterns. To closely represent the daily movement, historical customer load data with a fine time resolution is needed in order not to lose important information such as peak load and time. Fine time intervals may be made coarser, but the converse is not true. 
     Embodiments use historical customer load data segregated into customer type (commercial, residential and industrial) and sampled in a fine time resolution (e.g., every hour) for proper identification of load patterns. 
     Embodiments use a Gaussian random walk to simulate customer load curves.  FIG. 4  shows a Gaussian random walk. 
     The change in load in each time step is presumed to be an independently distributed normal random variable, such that the change in load level is
 
 P   t   −P   t−Δt   =ΔP   t =μ t +σ t Φ −1 ( u   t )  (1)
 
     where P t  is the load in kW at a time tε[0,24] hours and refers to the simulated value during simulation, μ t  is the mean and σ t  is the standard deviation of the load change from t−Δt to t, where Δt is the time resolution. Φ − (u t ) is the inverse cumulative distribution function of a standard normal distribution for probability u t ε[0,1]. μ 0  is defined as the mean and σ b  is defined as the standard deviation of the load at time t=0 (i.e., the starting load level P 0 ). Note that, Φ −1 (u t ) is a random variable that is distributed with a standard noLmal distribution for a given u t . During the simulation, u t  is randomly generated. Therefore, the right hand size of (1) will convert a standard normal random variable to a normal random variable that is distributed with mean μ t  and standard deviation σ t . 
       FIGS. 5 ,  5 A,  5 B and  5 C show the customer load curve simulation method. Customer load pattern parameters are calculated from observed, historical customer load data, which are then used to generate customer load curves. Once calculated, the customer load pattern parameters can be stored in memory and used in future simulations directly without recalculation. 
     Customer load pattern parameters are calculated from historical load data of a given group of customers. The load data comprises observations of daily load (e.g., through 30 days) at a desired sampling frequency (e.g., every 15 minutes). The historical load data may be acquired using, for example, smart electricity meters measuring power (e.g., kW) and energy (e.g., kWh), and stored in memory for input to a network or computer. 
     The historical load data and generated load pattern parameters (data) are denoted using “a hat” over the unit vector variables. {circumflex over (P)} ti  is historical load data obtained at time t in observation i, where i=1, 2, . . . , N, with N the number of load level observations (step  501 ). 
     After historical load data acquisition, load pattern parameters are generated. The starting load level P 0  has a mean {circumflex over (μ)} 0 =Σ i=1   N {circumflex over (P)} 0i /N and a standard deviation {circumflex over (σ)} 0 =√{square root over (Σ i=1   N ({circumflex over (P)} 0i −{circumflex over (μ)} 0 ) 2 /(N−1))} (step  503 ). Similarly, the load change ΔP t  from time t−Δt to time t has a mean {circumflex over (μ)} t =Σ i=1   N Δ{circumflex over (P)} ti /N and a standard deviation {circumflex over (σ)} t =√{square root over (Σ i=1   N (Δ{circumflex over (P)} ti −{circumflex over (μ)} t ) 2 /(N−1))}. {circumflex over (μ)} t  and {circumflex over (σ)} t  are calculated for each t=Δt, 2Δt, . . . , 24 and are represented, respectively, as vectors {circumflex over (μ)}=[{circumflex over (μ)} 0+Δt , . . . {circumflex over (μ)} 24 ] and {circumflex over (σ)}=[{circumflex over (σ)} 0+Δt , . . . {circumflex over (σ)} 24 ] (step  505 ). 
     Load pattern parameters {circumflex over (μ)} 0  and {circumflex over (σ)} 0 , vectors {circumflex over (μ)}=[{circumflex over (μ)} 0+Δt , . . . {circumflex over (μ)} 24 ] and {circumflex over (σ)}=[{circumflex over (σ)} 0+Δt , . . . {circumflex over (σ)} 24 ] from {circumflex over (μ)} t  and {circumflex over (σ)} t , and simulation parameters P {tilde over (t)}     start   , n sim  and [{tilde over (t)} start ,{tilde over (t)} end ,Δ{tilde over (t)}] are input (step  507 ). 
     P {tilde over (t)}     start    is the desired starting load level (e.g., Watts), n sim  is the number of simulated load curves that will be generated, {tilde over (t)} start  is the desired starting time of the simulation (typical default is 0), {tilde over (t)} end  is the desired ending time of the simulation (typical default is 24) and Δ{tilde over (t)} at is the desired time resolution for the simulation (typical default is 15 minutes). 
     The method presumes that future electrical load follows the pattern parameters with appropriate scaling of the calculated parameters and randomness of the calculated mean values. {circumflex over (μ)} 0  and {circumflex over (σ)} 0  are the mean and standard deviation of the starting load level of the pattern parameter data and {circumflex over (μ)} and {circumflex over (σ)} are vectors which represent the mean and the standard deviation of the load changes between each time period Δt of the pattern parameter data. {tilde over (t)} start  and {tilde over (t)} end  may be different from the pattern parameter data. For example, the pattern parameters can be estimated using daily data between hours 0 and 24, but the patterns and the load levels that are generated using them will be needed the most for the peak period between 12 PM to 4 PM. 
     The time frame (default is [0,24]) and time resolution Δt of the pattern parameter data may be converted to a desired time frame [{tilde over (t)} start ,{tilde over (t)} end ] and time resolution Δ{tilde over (t)} (steps  509 ,  511 ,  513 ). This operation is optional and is used if the time frame and/or time resolution of the historical load data (and therefore the pattern parameter data) is different from the desired time frame [{tilde over (t)} start ,{tilde over (t)} end ] and time resolution Δ{tilde over (t)}. For example, if the pattern parameter data was sampled every 15 minutes in a 24 hour period, but the desired time resolution for the simulation will be every 10 minutes between 12 PM to 4 PM, the pattern parameter time resolution must be converted to the desired time resolution. Using the existing pattern parameter data, a starting load at hour 12 PM is found and 10 minute samples for the desired time resolution are calculated using interpolation. Cumulative means Σ l=0   t {circumflex over (μ)} l  and variances Σ l=0   t {circumflex over (σ)} l   2  are found for each time point t in the pattern parameter data (step  509 ). Note that variance is the square of the standard deviation by definition. 
     The cumulative means Σ i=0   {tilde over (t)} {circumflex over (μ)} i  Σ i=0   {tilde over (t)} {circumflex over (σ)} i   2  of the pattern parameter data are interpolated to obtain the cumulative means {tilde over (μ)} and the cumulative standard deviation {tilde over (σ)} vectors for each discrete time point {tilde over (t)} value in the simulation interval {tilde over (t)}ε[{tilde over (t)} start ,{tilde over (t)} end ] with increments of Δ{tilde over (t)} (step  511 ). 
     For example, the mean of the starting load level at the beginning of the ({tilde over (t)} start ) is interpolated as {tilde over (μ)} t     start   =Σ i=0   t     1   {circumflex over (μ)} i +{circumflex over (μ)} t     2   ({tilde over (t)} start −t 1 )/(t 2 −t 1 ), where t 1  is the closest time point in the pattern parameter data that is smaller than {tilde over (t)} start , and t 2  is the closest time point in the pattern parameter data that is larger than {tilde over (t)} start . In other words, t 1 ≦{tilde over (t)} start ≦t 2  and t 2 =t 1 +Δt. The standard deviation of the starting load level is similarly found as σ {tilde over (t)}     start   =√{square root over (Σ i=0   t     1   {circumflex over (σ)} i   2 +{circumflex over (σ)} t     2   ({tilde over (t)} start −t 1 )/(t 2 −t 1 ))}. The same method is used to calculate each {tilde over (t)} in [{tilde over (t)} start ,{tilde over (t)} end ] with increments Δ{tilde over (t)}. 
     First differences are calculated to obtain the means μ {tilde over (t)} ={tilde over (μ)} {tilde over (t)} −{tilde over (μ)} {tilde over (t)}−Δ{tilde over (t)}  and standard deviations σ {tilde over (t)} =√{square root over ({tilde over (σ)} {tilde over (t)}   2 −{tilde over (σ)} {tilde over (t)}−Δ{tilde over (t)}   2 )} for each time increment Δ{tilde over (t)} of the simulation starting from {tilde over (t)}={tilde over (t)} start +Δ{tilde over (t)} to {tilde over (t)} end  (step  513 ). Therefore, the conversion from the pattern parameter data&#39;s time resolution Δt to the desired time resolution Δ{tilde over (t)} is complete, yielding the mean of the starting load level μ {tilde over (t)}     start   , standard deviation of the starting load level σ {tilde over (t)}     start    and vectors denoting the means of the increments [μ {tilde over (t)}     start+Δ{tilde over (t)}   , . . . , μ {tilde over (t)}   end ] and the standard deviation of the increments [σ {tilde over (t)}     start+Δ{tilde over (t)}   , . . . , σ {tilde over (t)}   end ]. 
     If P {tilde over (t)}     start   ≠μ {tilde over (t)}     start   , the pattern parameter data are scaled using scaling factor α=P {tilde over (t)}     start   /μ {tilde over (t)}     start    (step  515 ). This is performed in order to properly scale the parameters of the pattern and to the load levels that are to be simulated. For example, the pattern parameters may be calculated from a group of residential units that start the day, on average, with a load of 100 kW, but a desired simulation may be for a larger group which has an average starting load of 500 kW. In this case, a scaling parameter of α=5 is used. The mean μ {tilde over (t)}  and standard deviation σ {tilde over (t)}  parameters found during conversion (step  513 ) are multiplied by the scaling parameter α (step  517 ). This yields scaled parameters μ {tilde over (t)}     start    and σ {tilde over (t)}     start    for the starting load level, and vectors [μ {tilde over (t)}     start+Δ{tilde over (t)}   , . . . , μ {tilde over (t)}   end ] and [σ {tilde over (t)}     start+Δ{tilde over (t)}   , . . . , σ {tilde over (t)}   end ] showing parameters for each time increment. 
     Random numbers u {tilde over (t)} , that are uniformly distributed in [0,1] are generated. They are converted to standard normal random numbers by inverse transform sampling Φ −1 (u {tilde over (t)} ). This is the common way of generating random numbers that follow a standard normal distribution and are therefore not shown as an independent step in the method. 
     A load P j,{tilde over (t)} =P j,{tilde over (t)}−Δ{tilde over (t)} +μ {tilde over (t)} +Φ −1 (u {tilde over (t)} )σ {tilde over (t)}  is calculated for each simulation replication jε{1, . . . , n sim } and time point {tilde over (t)}ε[{tilde over (t)} start ,{tilde over (t)} end ] (step  519 ). This yields an n sim  by 
             (             t   ~     end     -       t   ~     start         Δ   ⁢           ⁢     t   ~         +   1     )         
load matrix P with entries for each (j,{tilde over (t)}) pair. Finally,
 
 P=[P   j,{tilde over (t)} ]  (2)
 
     is output (step  521 ). 
       FIG. 6  shows an example load matrix P output, which depicts 10 weekday replications for a commercial customer in July. There is a load matrix P for each time point {tilde over (t)} and replication j of the simulation. 
     The customer load curve simulation method has two advantages: 1) the calculated means (μ {tilde over (t)}     start   , [μ {tilde over (t)}     start+Δ{tilde over (t)}   , . . . , μ {tilde over (t)}   end ]) and standard deviations (σ {tilde over (t)}     start   , [σ {tilde over (t)}     start+Δ{tilde over (t)}   , . . . , σ {tilde over (t)}   end ]) represent the average shape of the load curve and the variation around it and 2) the Gaussian random walk provides smooth load transitions between time steps. 
     A customer&#39;s response to changes in the price of electricity is represented using the concept of price elasticity of demand. Elasticity is the measure of responsiveness in economics. The price elasticity ε of a customer&#39;s power demand is defined as 
     
       
         
           
             
               
                 
                   
                     ε 
                     = 
                     
                       
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           P 
                         
                         P 
                       
                       / 
                       
                         
                           Δ 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           r 
                         
                         r 
                       
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     where P is the load and r is the contracted rate (i.e., the price of electricity).  FIG. 7  shows an elasticity plot. The quantity shows the percentage of load change with respect to a one percent change in price. Since this measure uses ratios rather than quantities, it is dimensionless (free of units). 
     An important aspect of the customer load is that not all end uses of the electricity have the same responsiveness. For example, a consumer might use air conditioning less if the price is high, but not turn the refrigerator off. In this example, electricity consumption for air conditioning is more sensitive to price changes (i.e., elastic) than the electricity consumption for the refrigerator. In order to reflect this aspect, embodiments divide customer demand into categories with different price elasticities, and then aggregate the results from them once a price change occurs. From this point of view, this approach is similar to obtaining the aggregate demand curve from the demand curves of individuals.  FIG. 8  shows an aggregation of demand. 
     Note that price elasticity in the above setting refers to load shedding but not rescheduling of the consumption. Embodiments presume that the consumer uses less if the price of electricity increases, but not postpone the consumption to a later time in the day. For instance, a consumer may postpone using a washing machine to the evening, when the price of electricity is lower. On the contrary, it is not true for AC. This may be handled using the similar concept of cross price elasticity, in which percentage change in the consumption of one product is calculated against one percent change in the price of another product. Then, the electricity consumed at different times are considered as different products, even if they are used for the same reason. 
     Embodiments need to obtain an empirical demand curve and estimate the price elasticities. However, there is not enough reliable data for this. Embodiments use the concept with reasonably set parameters. 
       FIG. 9  shows the method for simulating customer response to changing electricity prices. The method requires inputs P 0 , r 0 , r, v, ε, f type  and s (step  901 ). P 0  is the corresponding load level, rb is the current price of electricity, r is the proposed price of electricity, the vector v=[v 1 , . . . , v n ] shows the proportion of each end use of electricity within the current consumption for each end use k=1, . . . , n, the vector ε=[ε 1 , . . . , ε n ] represents the price elasticity of each end use, the vector s=[s 1 , . . . , s n ] is the coefficient of variation which is the ratio of the standard deviation of the response to the mean load level and f type  specifies the type of the function that is used in the characterization of the customer&#39;s response. 
     The data calculated by the customer load curve simulation is input. For each entry P j,{tilde over (t)}  of the calculated load matrix P, the customer response simulation method is used to calculate a customer response. Note that if there is no price change, then P j,{tilde over (t)}  will remain unaffected. The inputs r 0 , r, v, ε also have a time dimension. For instance, there is an existing (i.e., current) electricity price for each time point {tilde over (t)}ε[{tilde over (t)} start ,{tilde over (t)} end ]. However, since the same method is going to be implemented for each j=1, . . . , n sim  and {tilde over (t)}ε[{tilde over (t)} start ,{tilde over (t)} end ] those subscripts are dropped to make the notation simpler. As the price of the electricity changes from r 0  to r, the customer responds to the price change by altering the load level P. Electricity is consumed for different end uses. If there are n uses, then the vector v=[v 1 , . . . , v n ] shows the proportion of each end use within the customer&#39;s current consumption. For example, v k =0.2 means end use kε{1, . . . , n} constitutes 20% percent of current consumption. Each end use has a different sensitivity to the changes in the price of electricity. Therefore, the vector ε=[ε 1 , . . . , ε n ] is the price elasticity of each end use. For example, from (3), ε i =0.1 means the power demand for end use k would decrease by 0.1% if the price of electricity increases by 1%. 
     Two types of demand functions are used to simulate customer response, a linear demand function or a constant elasticity demand function. The input parameter f type  specifies demand function type. 
     Randomness is introduced to the customer&#39;s response, s=[s 1 , . . . , s n ,] which is the coefficient of variation, i.e., the ratio of the standard deviation of the response to the mean load level for each end use k. 
     Embodiments presume that the end uses are not substitutable, therefore, the cross price elasticity among end uses are discarded. From a different perspective, the same end use over two different time periods (e.g., laundry at 2 PM or 10 PM) can be considered substitutable. Then, there is cross price elasticity between these two time periods. This effect is discarded because the simulation is for emergency DR and peak period loads are considered. Postponing the consumption from a time point in the peak period to another time point which is still within the peak period is negligible. Nevertheless, extending the model to these settings is straightforward. 
     Embodiments output customer load P depending on whether the function type f type  is constant elasticity (step  903 ) or linear (step  905 ). 
     A customer responds to a price change from r 0  to r, by changing P 0 v k , the initial load associated with each end use k. If the function type f type  is constant elasticity, the customer&#39;s expected (i.e., mean, average) load level for end use k is calculated as 
               P   k     =       P   0     ⁢         v   k     ⁡     (     r     r   0       )         ε   k               
(step  909 ).
 
     If the function type f type  is linear, the slope (m k ) of the individual demand curve of each end use k at the current price r 0  and load level P 0  is estimated, 
               m   k     =       r   0         P   0     ⁢     v   k     ⁢     ε   k               
(step  911 ). In this manner, the linear demand function for each end use is characterized as it is sufficient to characterize a linear function with a point on the function (P 0 ,r 0 ) and the slope m k . Then, the customer&#39;s expected load level for each end use k is calculated as
 
               P   k     =     max   ⁢       {     0   ,         P   0     ⁢     v   k       +       r   -     r   0         m   k           }     .             
The load cannot be negative (step  913 ).
 
     Randomness is introduced to customer response around the expected value. u k  is a random number that is uniformly distributed in [0,1] from which the load level for end use k is generated, which is a random number that is normally distributed with mean P k  and standard deviation P k s k . 
     The customer&#39;s final load level is the sum of the load levels of each individual end user, which is calculated as P=Σ k=1   n P k +Φ −1 (u k )P k s k ) (step  915 ). 
       FIG. 10  shows the estimated demand curves.  FIG. 11  shows the responses for different rates when the customer&#39;s response s k =3% for each k.  FIGS. 12 and 13  show price response over the course of a day, respectively with no randomization and with s k =3% for each k. 
     One or more embodiments of the present invention have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the invention. Accordingly, other embodiments are within the scope of the following claims.