Demand shaping in an electrical power grid

Abstract

A method (100) for demand shaping through load shedding and shifting in an electrical smart grid.

Claims

1. A method for controlling the flow of energy in an electrical power grid, comprising: a. obtaining an optimal load shedding profile for a user of the grid; b. obtaining an optimal billing structure and an optimal load shifting profile for the user; c. iterating the steps of obtaining an optimal load shedding profile and obtaining an optimal billing structure and optimal load shifting profile until convergence is achieved; and d. controlling the energy consumption of the user based on the optimal load shedding profile, optimal billing structure, and optimal load shifting profile.

2. The method of claim 1, wherein obtaining an optimal load shedding profile utilizes a known billing structure and a known load shifting profile for each user of the grid.

3. The method of claim 2, wherein obtaining an optimal load shedding profile is performed without express information of load preferences of each user of the grid.

4. The method of claim 1, wherein each of obtaining an optimal load shedding profile and obtaining an optimal billing structure and an optimal load shifting profile is a solution to a robust convex optimization problem.

5. A method of controlling the energy demand of a user of an electrical power network, comprising: a. determining an optimal amount of user energy demand to shed and a price incentive required to achieve said optimal amount; b. determining a billing structure and an optimal amount of user energy demand to time-shift; c. repeating the first and second method steps until an equilibrium is achieved; and d. controlling the energy demand of the user based the determined optimal amount of user energy demand to shed, price incentive required to achieve said optimal amount, billing structure, and optimal amount of user energy demand to time-shift.

6. The method of claim 5, wherein determining the optimal amount of user energy demand to shed comprises formulating and solving an optimization problem using a known billing structure and a known user consumption for shiftable loads.

7. The method of claim 6, wherein formulating and solving comprises utilizing standard convex optimization techniques to solve the optimization problem in the case that energy price incentives are set by a party other than a network operator.

8. The method of claim 6, wherein formulating and solving comprises setting a price offer for the user and re-formulating and solving the re-formulated optimization problem with user preferences to load shedding in the case that energy price incentives are set by a network operator.

9. The method of claim 5, wherein determining a billing structure and an optimal amount of user energy demand to time-shift comprises formulating and solving an optimization problem based on results for calculating an optimal amount of user energy demand to shed and a price incentive required to achieve said optimal amount.

10. The method of claim 9, wherein formulating and solving comprises formulating and solving respective billing structures for each of user shiftable plus sheddable load and user non-moveable load using a fixed utility ratio.

11. The method of claim 10, wherein formulating and solving respective billing structures comprise re-formulating and solving re-formulated optimization problem to obtain an optimal amount of user energy demand to time-shift.

Description

DESCRIPTION OF THE DRAWINGS

(1) For a better understanding of the present invention, reference is made to the following description of an exemplary embodiment thereof, and to the accompanying drawings, wherein:

(2) FIG. 1 is a block diagram of a smart electrical power network (simplified);

(3) FIG. 2 is a schematic representation of a control method implemented in accordance with the present invention;

(4) FIG. 3a shows histograms of synthetic data of monthly demand distribution for a user;

(5) FIG. 3b shows histograms of synthetic data of hourly demand distribution for a user;

(6) FIG. 3c shows histograms of synthetic data of single day of the week and single hour demand distributions for multiple users;

(7) FIG. 3d shows synthetic data demand distribution comparison (same day of the week and hour) with empirical data for a user;

(8) FIG. 3e is a chart of simulated average user demand for each user cluster for one day;

(9) FIG. 3f is a chart of the shiftable, sheddable, and non-moveable demand for a user cluster;

(10) FIG. 3g shows a demand pattern comparison (shiftable, sheddable, and non-moveable demand in the same day) between two users in the same user cluster;

(11) FIG. 3h is a chart of exemplary DAM and RTM energy procurement costs for a single day;

(12) FIG. 3i is a chart of the synthetic demand pattern (shiftable, sheddable, and non-moveable demand) used by a base case in a simulation comparison with the method of FIG. 2;

(13) FIG. 3j is a chart of the DAM and RTM energy procurement costs used by a base case in a simulation comparison with the method of FIG. 2;

(14) FIG. 3k is a chart of the demand profile resulting from the use of the method of FIG. 2 in a simulation comparison with a base case; and

(15) FIG. 3l is a chart of the procurement strategy resulting from the use of the method of FIG. 2 in a simulation comparison with a base case.

DETAILED DESCRIPTION

(16) FIG. 1 is a block diagram of a typical smart power network 10 (simplified) that operates to provide intelligent control and two-way communication among the various network 10 components. The network 10 comprises a generating station(s) 12 that generates electricity/electrical power from an energy source, such as a renewable source (e.g., water, wind, etc.) or a non-renewable source (e.g., coal, natural gas, etc.). The generating station 12 may be configured in many ways and may utilize a combination of various energy sources. The generating station 12 may also utilize energy that was stored for later use. The network 10 also comprises a transmission system 14 that converts the generated electrical power into high voltage power (via a transmission substation 14a) and transfers the high voltage power over a long distance (via electrical transmission lines 14b) to a distribution system 16. The distribution system 16 (via a power substation 16a) converts the high voltage power into power at different lower voltages and splits, or distributes, these lower-voltage power flows off in multiple directions. The distribution system 16 may further modulate and regulate the lower-voltage power flows. Ultimately, the distribution system 16 transfers the distributed power (via distribution lines 16b) to electricity consumers or users 18. The users 18 are typically categorized as residential users 18a, commercial users 18b, and industrial users 18c. In many locations, the entities that generate electricity are separate from the entities that locally distribute electricity (such as a local utility company) so that the local distributors (with their associated distribution systems 16) may also be categorized as a user though not a final or end user (or customer). In some network 10 configurations, the distribution system 16 may also transfer the distributed power to sub-distribution systems (not shown) that further distribute the electricity to final or end users 18.

(17) Each of the network 10 components will comprise and utilize appropriate equipment (for example, sensors, controls, computers, automation, smart meters, etc.) that operates to implement the “smart” functionality of the network 10. The network 10 also comprises an operations center 20 that permits intelligent management and control of the various network 10 components and the overall network 10. The operations center 20 comprises a two-way communications network 22 to connect to the various network 10 components for gathering and acting on information. The operations center 20 also comprises a computer system 24 to provide, among other functions, monitoring, reporting, supervising, and intelligence functions for managing and controlling the various network 10 components and the overall network 10.

(18) The various components of the network 10 are well understood smart grid components. However, the computer system 24 of the operations center 20 is also adapted to permit the network 10 to operate and to implement methods in accordance with the present invention. FIG. 2 is a schematic representation of a control method 100 implemented in accordance with the present invention.

(19) As will be explained in more detail below, the control method 100 is a two step iterative process. In the first step, the load shedding problem is solved by assuming a known billing structure and a load shifting profile for each user/customer. This is done without knowing explicitly the preferences of each user, but still reaching a solution that is optimal for them. The second step uses this result to find the optimal billing structure and load shifting profile. The method 100 iterates between these two steps, reducing the cost of the solution at each step and thus reaching optimality. Instead of solving a complex optimization problem, the method 100 breaks the problem into two simpler robust convex optimization problems that can be even solved in a distributed fashion, thus reducing data transmission and processing requirements.

(20) The framework for methods and systems of the present invention, including the above control method 100, is described as follows. The objective is to build a mechanism where a utility company has direct control of part of the energy consumption of their final customers given the energy they can procure from generators in the market. The present invention considers that the utility company has two tools to control the energy demand: demand shifting and demand shedding, and, further, that it procures the energy it delivers from a two-stage market with a day-ahead market (DAM) and a real-time market (RTM). Currently, there are no methods or systems that involve two different mechanisms simultaneously to shape the demand, nor one that accounts for the energy procurement in a two-level market.

(21) In formulating a demand model, the present invention assumes that the utility company can monitor and control each individual customer, or cluster of customers, and that these customers can categorize their energy consumption in three different groups: loads that can be shed, those that can be shifted, and those that can neither be shed nor shifted, i.e., the non-movable loads. Although it sounds reasonable that the consumption of each of the different client's appliances are assigned to a single load category, it is important to note that there are many cases in which this is not the case. For example, take the energy consumption of an air conditioning (A/C) unit. The energy consumption level required for the A/C unit to give a minimum level of comfort could be categorized as a non-moveable load, whereas the additional load required to get to the current temperature desired by the customer can be categorized as a sheddable one. The present invention also assumes that the time window the utility company wants to optimize can be divided into T equal-sized time steps, and, without loss of generality, it assumes T to be one day.

(22) Continuing with the model, let q.sub.it, with i={1, . . . , n} and t={1, . . . , T} denote the total load of user i at time-step t; and thus q.sub.i={q.sub.i1 . . . q.sub.iT} denotes the total demand profile for user i for that one day. At each time-step t, the demand of user i is composed by a sheddable demand q.sub.it.sup.(d), a shiftable demand q.sub.it.sup.(s), and a non-sheddable, non-shiftable demand q.sub.it.sup.(n), i.e., q.sub.it=q.sub.it.sup.(d)+q.sub.it.sup.(s)+q.sub.it.sup.(n).

(23) As noted above, the sheddable demand is the part of the user's demand that the utility company can shed, making the user incur a personal loss of value. For example, reducing the number of lights used to light an area or increasing the target temperature of the A/C unit will drop the total consumption, bringing certain levels of discomfort to the user. In all these cases, the demand is not shifted/postponed to a different time.

(24) Thus, let y.sub.it denote the amount of load user i sheds at time t, and y.sub.i denote the shedding profile vector for user i. Given that at time-step t user i has a sheddable load of q.sub.it.sup.(d) (i.e., the maximum amount of sheddable load) then it follows that y.sub.it≦q.sub.it.sup.(d). In addition, let D.sub.i(y.sub.it) denote the discomfort cost for user i for shedding a load y.sub.it, with D.sub.i:.sub.30.fwdarw..sub.+. To incentivize this shedding, the utility company gives a compensation of p.sub.it.sup.(d) to user i at time-step t per unit of load shed. Given the incentive p.sub.it.sup.(d) per unit of load, user i will decide to shed an amount y.sub.it=S.sub.i(p.sub.it.sup.(d)), with S.sub.i:.sub.+.fwdarw..sub.+. In the article by L. Chen, N. Li, S. H. Low and J. C. Doyle referenced above, the authors use a similar model where they assumed that the discomfort function D.sub.i(y.sub.it) is a continuous, increasing, strictly convex function, with D.sub.i(0)=0, and that the shedding function is of the form y.sub.it=S.sub.i(p.sub.it.sup.(d))=a.sub.ip.sub.it.sup.(d) where a.sub.i is a user-dependent parameter.

(25) As noted above, the shiftable demand is the part of the client's demand that can be shifted during the day. Items like washing machines, dish washers, and charging electric vehicles fall in this category since there is a window of time for using these machines and not a single specific point in the day. All these appliances can be used in a certain interval within the T time-step window. Since the shiftable demand cannot be shed, the total shiftable demand for user i has a fixed value Σ.sub.tq.sub.it.sup.(s)={circumflex over (q)}.sub.i.sup.(s). The shiftable demand is composed by the consumption requirements of m different appliances, where the consumption of appliance j of user i at time t is given by x.sub.ijt and the time profile of the appliance is denoted by x.sub.ij. Thus, at any given point in time the total consumption of user i is given by q.sub.it.sup.(s)=Σ.sub.jx.sub.ijt.

(26) The total consumption of an appliance will be a given value denoted by Σ.sub.tx.sub.ijt={circumflex over (x)}.sub.ij, with a maximum value of x.sub.ij and a minimum of x.sub.ij at any given time-step t. For every appliance j, user i sets a time window [α.sub.ij, β.sub.ij] in which the user i requires the appliance to run. For example, an electric vehicle could charge for 3 hours at any point between 10:00 p.m. and 7:00 a.m. This also implies that x.sub.ijt=0, ∀t.Math.└α.sub.ij, β.sub.ij┘. In the article by A.-H. Mohsenian-Rad, et al. referenced above, the authors use a similar model for shifting user loads to achieve demand response. They also add a minimum stand-by requirement for x.sub.ijt, but since that can be considered as a non-moveable consumption the present framework adds that to q.sub.ijt.sup.(n) and keeps only the shiftable portion in x.sub.ijt.

(27) The non-moveable demand q.sub.it.sup.(n), as its name suggests, is the part of the demand that cannot be modified or has not been programmed ahead. The user requires that amount of energy at that specific time. Examples of this are the demand of a refrigerator at the lowest viable temperature setting or the on position in the A/C unit. Given that this energy consumption is not necessarily known ahead, q.sub.i.sup.(n) can be considered as a stochastic process. Further, it can assumed that the user, or at least its smart meter, has some forecast {tilde over (q)}.sub.i.sup.(n) of the non-moveable demand for the whole time-window at the beginning of the period. Using the ideas in the article by S. Meyn, M. Negrete-Pincetic, G. Wang, A. Kowli, and E. Hafieepoorfard referenced above, the error between the meter's forecast and the real static consumption {circumflex over (q)}.sub.i.sup.(n), i.e., ε.sub.i.sup.(n)={tilde over (q)}.sub.i.sup.(n)−{circumflex over (q)}.sub.i.sup.(n), can be modeled as a driftless Brownian motion with variance σ.sub.i.sup.(n).

(28) The demand model also addresses the billing system. At time-step t, user i is billed an amount B.sub.it.sup.S(q.sub.it) by the utility company, where B.sub.it.sup.S(q.sub.it):.sub.+.fwdarw..sub.+. In order to incentivize customers to set part of their demand as shiftable or sheddable demand, these loads are charged with a different billing structure B.sub.it.sup.f(q.sub.it), where B.sub.it.sup.f(q.sub.it):.sub.+.fwdarw..sub.+. Since the utility company will be able to control these loads directly, it is assumed that the billing structure does not depend on t, i.e., B.sub.it.sup.f(q.sub.it)=B.sub.t.sup.f(q.sub.it), ∀t (∀ being a universal quantifier signifying “given any” or “for all”).

(29) In formulating the procurement model, the utility company must be able to satisfy the total demand Q.sub.t=Σ.sub.iq.sub.it at every time-step t, and thus deliver a demand profile Q=[Q.sub.1 . . . Q.sub.T]. The utility company satisfies this demand by procuring this energy from a two-stage market composed by a day-ahead market (DAM) whose purchases are closed a day in advance of the actual demand, and a real-time market (RTM) that clears close to real-time. The prices in these two markets can differ significantly. Let C.sub.t.sup.d(Q.sub.t), with C.sub.t.sup.d(Q.sub.t):.sup.2.sub.+.fwdarw..sub.+ denote the total cost of procuring energy from the DAM and delivering to satisfy a demand of Q.sub.t at time t. Similarly, let C.sub.t.sup.S(Q.sub.t), with C.sub.t.sup.S(Q.sub.t):.sup.2.sub.+.fwdarw..sub.+ denote the total cost of procuring energy from the RTM and delivering to satisfy a demand of Q.sub.t at time t. In general, we can assume that C.sub.t.sup.d(Q.sub.t)≦C.sub.t.sup.S(Q.sub.t) with C.sub.t.sup.d(Q.sub.t) being known by the utility company and C.sub.t.sup.S(Q.sub.t) being an unknown stochastic process that represents the variation in prices that the RTM shows.

(30) Within the framework described above, there are several different problems of interest for a utility company that can be addressed. The first is the cost/revenue problem. Given the cost of procuring and delivering load, C.sub.t.sup.d(Q.sub.t) and C.sub.t.sup.S(Q.sub.t), a billing structure B.sub.it.sup.f(q.sub.it) and B.sub.it.sup.S(q.sub.it), ∀i, t and incentive prices per unit of load p.sub.i.sup.(d), ∀i, what is the optimal load delivery that maximizes the utility company's revenue and satisfies the demand. This is similar to the problems posed in the articles by A.-H. Mohsenian-Rad et al. and L. Chen, N. Li, S. H. Low and J. C. Doyle, where the billing structure B.sub.it.sup.f(q.sub.it) is given and the authors are interested in computing the demand pattern Q that maximizes the revenues.

(31) Another is the PAR problem that seeks to minimize the asset requirements to deliver peak-energy, which can be achieved by minimizing the peak-to-average ratio (PAR). This problem can be written as follows (equation 1):

(32) $\begin{array}{cc}\underset{X,Y}{\min }\underset{t}{\max }\frac{{\mathrm{TQ}}_t}{\underset{t}{.Math.}{Q}_t}s.t.\text{:}\begin{array}{cc}{Q}_t=\underset{i}{.Math.}{q}_{\mathrm{it}},& \forall t,\\ q_{\mathrm{it}}=\left(q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}}\right)+q_{\mathrm{it}}^{\left(s\right)}+{\stackrel{\_}{q}}_{\mathrm{it}}^{\left(n\right)},& \forall t,\forall i\\ q_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\forall i\\ \underset{t}{.Math.}{x}_{\mathrm{ijt}}={\hat{x}}_{\mathrm{ij}},& \forall j,\forall i\\ x_{\mathrm{ijt}}=0,& \forall j,\forall i,\forall t.Math.\left[{\alpha }_{\mathrm{ij}},{\beta }_{\mathrm{ij}}\right]\\ y_{\mathrm{it}}\ge 0,& \forall t,\forall i\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall t,\forall i\\ x_{\mathrm{ijt}}\ge {\underset{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\\ x_{\mathrm{ijt}}\le {\stackrel{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i,\end{array}& \left(1\right)\end{array}$
where {tilde over (q)}.sub.i.sup.(n) the third line of equation 1) denotes the forecast the utility company has of the non-moveable demand for the next day. This problem could also be posed as a robust optimization problem, considering the error ε.sub.it.sup.(n) of the forecast for each user.

(33) A more sophisticated problem is, given energy costs, billing structures, and incentive prices, to compute a solution that achieves a demand profile Ψ=[Ψ.sub.1 . . . Ψ.sub.T] with a minimum cost. This is the demand shaping problem. For a given error parameter Δ, this problem can be written as follows (equation 2):

(34) $\begin{array}{cc}\underset{X,Y}{\min }\underset{t}{.Math.}{c}_t\left(Q_t\right)+\underset{i}{.Math.}\underset{t}{.Math.}{p}^{\left(d\right)}{y}_{\mathrm{it}}-\underset{i}{.Math.}\underset{t}{.Math.}{B}_t\left(q_{\mathrm{it}}\right)+\Delta {.Math.Q-\Psi .Math.}_{p}s.t.\text{:}\begin{array}{cc}{Q}_t=\underset{i}{.Math.}{q}_{\mathrm{it}},& \forall t,\\ q_{\mathrm{it}}=\left(q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}}\right)+q_{\mathrm{it}}^{\left(s\right)}+{\stackrel{\_}{q}}_{\mathrm{it}}^{\left(n\right)},& \forall t,\forall i\\ q_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\forall i\\ \underset{t}{.Math.}{x}_{\mathrm{ijt}}={\hat{x}}_{\mathrm{ij}},& \forall j,\forall i\\ x_{\mathrm{ijt}}=0,& \forall j,\forall i,\forall t.Math.\left[{\alpha }_{\mathrm{ij}},{\beta }_{\mathrm{ij}}\right]\\ y_{\mathrm{it}}\ge 0,& \forall t,\forall i\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall t,\forall i\\ x_{\mathrm{ijt}}\ge {\underset{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\\ x_{\mathrm{ijt}}\le {\stackrel{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\end{array}& \left(2\right)\end{array}$

(35) The present invention primarily focuses on the cost/revenue problem. The objective is to compute an optimal schedule for the shiftable demand and the sheddable demand, such that the revenues are maximized given cost, billing, and price incentive structures. The present invention assumes several requirements for this problem. One requirement is fairness. Specifically, the billing structure must be fair. If two customers consume the same amount at the same time, then the billing they should receive shouldn't differ even though their usage utility functions might differ. This means B.sub.i.sup.f(q.sub.it)=B.sup.f(q.sub.it) and B.sub.it.sup.S(q.sub.it)=B.sub.t.sup.S(q.sub.it), ∀i. Similarly, if two users are requested to shed the same amount of demand at the same time period t, then the payment they receive should be the same, that is p.sub.it.sup.(d)=p.sub.t.sup.(d), ∀i. Another requirement is proportionality. The present invention assumes the more energy a customer consumes, the larger the billing price should be, that is, B.sup.f(q.sub.it1)≧B.sup.f(q.sub.it2) if q.sub.it1≧q.sub.it2. Furthermore, the present invention assumes proportionality, that is, B.sup.f(q.sub.it1)/B.sup.f(q.sub.it2)=q.sub.it1/q.sub.it2. This is certainly justified for small consumers/households, whereas it might not be true for larger industrial consumers that pay additional fees for peak consumption, load factors, and other factors. The same holds true for B.sub.t.sup.S(q.sub.it).

(36) Another requirement is positive revenues. The present invention assumes that the revenues for delivering the required demand Q are positive, that is, Σ.sub.tB.sup.f(Q.sub.t.sup.1)+B.sub.t.sup.S(Q.sub.t.sup.2)≧Σ.sub.tC.sub.t.sup.d(Q.sub.t.sup.1)+C.sub.t.sup.S(Q.sub.t.sup.2). Another requirement relates to procurement. Let {tilde over (q)}.sub.i.sup.(n)=[{tilde over (q)}.sub.i1.sup.(n) . . . {tilde over (q)}.sub.iT.sup.(n)] denote an optimal procurement of the non-moveable demand of user i in the DAM. The present invention assumes that energy procurement is done in the following way. Since the utility company knows the shiftable demand and the sheddable demand a priori (q.sub.i.sup.(S) and q.sub.i.sup.(d)) and an optimal procurement of the non-moveable demand {tilde over (q)}.sub.i.sup.(n), it is assumed that it purchases a total demand {tilde over (Q)}=[{tilde over (Q)}.sub.1 . . . {tilde over (Q)}.sub.T]{tilde over (Q)}.sub.t=Σ.sub.t(q.sub.it.sup.(d)+q.sub.it.sup.(s)+{tilde over (q)}.sub.it.sup.(n) in the DAM, and it has a cost of C.sub.t.sup.d({tilde over (Q)}.sub.t) to deliver it to the final users. The difference between the purchased energy {tilde over (Q)}.sub.t.sup.(n) and the real non-moveable demand Q.sub.t.sup.(n) is purchased in the RTM and it has a total delivery cost of C.sub.t.sup.S((Q.sub.t.sup.(n)−{tilde over (Q)}.sub.t.sup.(n)).sup.+). Note that this framework can also account for the spinning reserve (i.e., the generation capacity that is online but unloaded and that can respond quickly to compensate for generation or transmission outages). This is done by purchasing in the DAM not only {tilde over (Q)}.sub.t, but (1+ε){tilde over (Q)}.sub.t where ε denotes the additional fraction purchased as spinning reserve. Then, the total cost of procuring the energy will be given by C.sub.t.sup.d((1+ε){tilde over (Q)}.sub.t)+C.sub.t.sup.S((Q.sub.t.sup.(n)−(1+ε){tilde over (Q)}.sub.t.sup.(n)−εQ.sub.t.sup.(d)−εQ.sub.t.sup.(s)).sup.+).

(37) Another requirement relates to the shedding function. Most of the current literature defines a simple shedding function S.sub.i(p.sub.it.sup.(d)) since that simplifies the analysis significantly. For example, in the article by A.-H. Mohsenian-Rad, et al., the authors use y.sub.it=S.sub.i(p.sub.it.sup.(d))=a.sub.i(p.sub.it.sup.(d)), where a.sub.i is a user-dependent parameter, and then proceed to compute an a.sub.i that maximizes the user's revenues. There are two main issues with this approach. First, the user is changing its preference every time a price is given and, second, it is already known that the user has a discomfort cost D.sub.i(y.sub.it) which can be used to find the optimal load y.sub.it that will maximize the user's revenues. Hence, instead of using a simplification for the shedding function, the present invention solves the user's optimization problem to find the optimal load to shed.

(38) With these requirements assumed, the present invention addresses the cost/revenue problem in the following way. For simplicity, let B.sub.t(q.sub.it.sup.1, q.sub.it.sup.2)=B.sup.f(q.sub.it.sup.1)+B.sub.t.sup.S(q.sub.it.sup.2) and C.sub.t(Q.sub.t.sup.1, Q.sub.t.sup.2)=C.sub.t.sup.d(Q.sub.t.sup.1)+C.sub.t.sup.S(Q.sub.t.sup.2). Similarly, let X=[x′.sub.11 . . . x′.sub.nm]′ and Y=[y′.sub.1 . . . y′.sub.n]′. With this, the optimization problem for the utility company becomes as follows (equation 3):

(39) $\begin{array}{cc}\underset{X,Y,{\tilde{Q}}_t^{\left(n\right)}}{\max }\underset{t}{.Math.}\underset{i}{.Math.}\left[B_t\left(q_{\mathrm{it}}^{\left(s\right)}+q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}},q_{\mathrm{it}}^{\left(n\right)}\right)-p_t^{\left(d\right)}{y}_{\mathrm{it}}\right]-C_t\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)},Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right)s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\forall i\\ \underset{t}{.Math.}{x}_{\mathrm{ijt}}={\hat{x}}_{\mathrm{ij}},& \forall j,\forall i\\ x_{\mathrm{ijt}}=0,& \forall j,\forall i,\forall t.Math.\left[{\alpha }_{\mathrm{ij}},{\beta }_{\mathrm{ij}}\right]\\ y_{\mathrm{it}}\ge 0,& \forall t,\forall i\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall t,\forall i\\ x_{\mathrm{ijt}}\ge {\underset{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\\ x_{\mathrm{ijt}}\le {\stackrel{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\end{array}& \left(3\right)\end{array}$

(40) This optimization problem assumes that the actual procurement cost in the RTM C.sub.t.sup.S(Q.sub.t) and the actual non-moveable demand Q.sub.t.sup.(n) are known a priori. Since that is seldom the case, a robust optimization approach is used to compute an optimal solution. The idea behind a robust optimization formulation is to consider not a single expected value of a stochastic variable, but a whole set of possible values from which the approach needs to cover the solution. The amount of robustness in the solution is then controlled by the size of the uncertainty set in which the stochastic variables are considered to be included.

(41) Given that there are two main stochastic variables in the above optimization problem (equation 3), then the optimization problem can have robustness in two different ways. First, there can be robustness in procurement costs. The RTM prices can be considered as a stochastic process, hence, the problem becomes how much energy to purchase in the DAM given that there is uncertainty in the prices in the RTM. In this case, the uncertainty set (symbolized as calligraphic C)
C.sub.t, for t={1, . . . , T}
will represent the possible cost functions the utility company will get in the RTM at time step t during the next day. Second, there can be robustness in consumption level. The second robust optimization formulation is achieved by considering the stochastic process behind the non-moveable demand. In this case
Q.sub.t.sup.(n), for t={1, . . . , T}
denotes the uncertainty set (symbolized as calligraphic Q) for the non-moveable demand Q.sub.t.sup.(n). The constraint that demand is within a certain convex set is added and an optimal solution for all demands within that set is computed.

(42) Both settings can be combined and optimal solutions can be computed when the procurement cost and the consumption are unknown. This might lead to very conservative solutions since the formulation accounts for two sources of uncertainty. The full robust problem formulation is then given by (equation 4):

(43) $\begin{array}{cc}\underset{X,Y,{\tilde{Q}}_t^{\left(n\right)}}{\max }\underset{C_t,Q_t^{\left(n\right)}}{\min }\underset{t}{.Math.}\underset{i}{.Math.}\left[B_t\left(q_{\mathrm{it}}^{\left(s\right)}+q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}},q_{\mathrm{it}}^{\left(n\right)}\right)-p_t^{\left(d\right)}{y}_{\mathrm{it}}\right]-C_t\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)},Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right)s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\forall i\\ \underset{t}{.Math.}{x}_{\mathrm{ijt}}={\hat{x}}_{\mathrm{ij}},& \forall j,\forall i\\ x_{\mathrm{ijt}}=0,& \forall j,\forall i,\forall t.Math.\left[{\alpha }_{\mathrm{ij}},{\beta }_{\mathrm{ij}}\right]\\ C_t\in C_t,& \forall t,\\ Q_t^{\left(n\right)}\in Q_t^{\left(n\right)},& \forall t,\\ y_{\mathrm{it}}\ge 0,& \forall t,\forall i\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall t,\forall i\\ x_{\mathrm{ijt}}\ge {\underset{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\\ x_{\mathrm{ijt}}\le {\stackrel{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\end{array}& \left(4\right)\end{array}$

(44) The level of robustness in the above problem (equation 4), and thus how conservative is the solution, will depend on the uncertainty sets
C.sub.t and Q.sub.t.sup.(n).

(45) In general, these uncertainty sets can be either boxes, which basically add a set of box constraints to the optimization model, or ellipsoids. The main advantage of box constraints is its ease of implementation since most of the time that is translated to interval constraints for each variable. The drawback is that the corners of the box make the solution overly conservative since they are much further away from the centre of the box compared to the sides. On the other hand, ellipsoid sets give better solutions and don't suffer from being overly-conservative. The disadvantage is that they are harder to implement since they are not linear and can only be approximated if a linear model is desired.

(46) Considering the cost/revenue problem detailed above, there is also a need to design a mechanism that would determine the billing structure B.sub.t(q.sub.it.sup.1, q.sub.it.sup.2) and the price incentive structure p.sub.t.sup.(d)), such that the whole system achieves equilibria. That is, given the structure computed by this mechanism, the optimal solution obtained by the utility company will give no incentives to the final customers to change their solutions or cheat.

(47) As a separate issue, to understand what will be the objective of the final customer, the optimization problem they see must be posed. At time-step t, user i is billed an amount B.sub.t(q.sub.it.sup.(s)+q.sub.it.sup.(d)+q.sub.it.sup.(n)) by the utility company, and is paid p.sub.it.sup.(d) per unit of load shed, hence the user will solve the following optimization problem (equation 5),

(48) $\begin{array}{cc}\underset{x_i,y_i}{\min }\underset{t=1}{\overset{T}{.Math.}}{B}_t\left(q_{\mathrm{it}}^{\left(s\right)}+q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}},q_{\mathrm{it}}^{\left(n\right)}\right)+D_{\mathrm{it}}\left(y_{\mathrm{it}}\right)-p_{\mathrm{it}}^{\left(d\right)}{y}_{\mathrm{it}}s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\\ \underset{t}{.Math.}{x}_{\mathrm{ijt}}={\hat{x}}_{\mathrm{ij}},& \forall j,\\ x_{\mathrm{ijt}}=0,& \forall j,\forall t.Math.\left[{\alpha }_{\mathrm{ij}},{\beta }_{\mathrm{ij}}\right]\\ y_{\mathrm{it}}\ge 0,& \forall t,\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall t,\\ x_{\mathrm{ijt}}\ge {\underset{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\\ x_{\mathrm{ijt}}\le {\stackrel{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\end{array}& \left(5\right)\end{array}$

(49) In summary, the cost/revenue problem has two basic parts that need to be addressed. One part is that prices and billing structures that achieve some equilibrium must be computed, and, in a second part, the allocation of demand for each user given these prices and billing structure must be computed. The control method 100 solves this problem by performing a two-step process which iterates between computing a billing structure and computing the price incentive to achieve an equilibrium.

(50) To analyze how to solve the general version of this problem posed in equation 4, additional structure is pointed out. First, it is assumed that the general billing structure is linear, that is, B.sup.f(q.sub.it)=b.sup.fq.sub.it and B.sub.t.sup.S(q.sub.it)=b.sub.t.sup.Sq.sub.it, for all t, where b.sup.f and b.sub.t.sup.S and are constants. Although this is not the case for large industrial consumers, most residential households do have a linear billing cost. As noted above, since the utility company has direct control over the shiftable demand and the sheddable demand, it is assumed that the billing structure for these demands is constant throughout the day, i.e., b.sub.t.sup.f=b.sup.f, for all t.

(51) The method 100 comprises a first step 102 that addresses price incentive and shedding calculation. The objective of this first step 102 is to compute the optimal amount of demand to shed and the price incentive required to do so. It is assumed that the billing structure B.sub.t(Q.sub.t.sup.1, Q.sub.t.sup.2) is known for all t, as well as the solution X for the shiftable loads. With these assumptions, and given p.sup.(d), the utility company's problem (given by equation 4) simplifies to the following (equation 6):

(52) $\begin{array}{cc}\underset{Y,{\tilde{Q}}_t^{\left(n\right)}}{\max }\underset{C_t,Q_t^{\left(n\right)}}{\min }\underset{t}{.Math.}\underset{i}{.Math.}\left[b^f\left(q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}}\right)+b_t^s{q}_{\mathrm{it}}^{\left(n\right)}-p_t^{\left(d\right)}{y}_{\mathrm{it}}\right]-C_t\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)},Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right)s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\forall i\\ C_t\in C_t,& \forall t,\\ Q_t^{\left(n\right)}\in Q_t^{\left(n\right)},& \forall t,\\ y_{\mathrm{it}}\ge 0,& \forall t,\forall i\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall t,\forall i,\end{array}& \left(6\right)\end{array}$
which is separable in t; hence, for every time step t the following optimization problem needs to be solved (step 104) (equation 7):

(53) $\begin{array}{cc}\underset{Y_t,{\tilde{Q}}_t^{\left(n\right)}}{\max }\underset{C_t,Q_t^{\left(n\right)}}{\min }\underset{i}{.Math.}\left[b^f\left(q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}}\right)+b_t^s{q}_{\mathrm{it}}^{\left(n\right)}-p_t^{\left(d\right)}{y}_{\mathrm{it}}\right]-C_t^d\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}\right)-C_t^s\left(Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right)s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall i\\ C_t\in C_t,& \\ Q_t^{\left(n\right)}\in Q_t^{\left(n\right)},& \\ y_{\mathrm{it}}\ge 0,& \forall i\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall i\end{array}& \left(7\right)\end{array}$

(54) In the setting that the price incentive is set externally (e.g., by government regulations, contracts, etc.), and the utility company can shed the load directly, this optimization problem can be solved using standard convex optimization techniques (step 104a), and the method 100 can move to a second step 110 for billing and shift calculation.

(55) In the setting that the utility company has to set the price so that the final users are motivated to shed their load, the utility company has no direct control over the load, and can only set the price offer so that respective pre-programmed smart meters at a user take the decision of what amount of load to shed depending on the user's preferences and loss function (step 104b). From the user's point of view, given a billing structure, a solution for the shiftable load, and a price incentive p.sub.t.sup.(d), the optimization problem for every time step t is given by (equation 8):

(56) $\begin{array}{cc}\underset{y_i}{\min }{b}^f\left(-y_{\mathrm{it}}\right)+D_{\mathrm{it}}\left(y_{\mathrm{it}}\right)-p_t^{\left(d\right)}{y}_{\mathrm{it}},s.t.\text{:}{y}_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},y_{\mathrm{it}}\ge 0.& \left(8\right)\end{array}$
which is easily solvable once D.sub.it is known. If it is assumed D.sub.it is continuous, increasing, and differentiable, then the solution for this user optimization problem (equation 8) is given by (equation 9):
y*.sub.it=S.sub.it(p.sub.t.sup.(d))=min{q.sub.it.sup.(d),max{0,D′.sub.it.sup.−1(p.sub.t.sup.(d)+b.sup.f)}}  (9)

(57) This computes the optimal value of y.sub.it for each user i at each time step t, and thus can be considered as the shedding function for user i at time t. Let S.sub.t(p.sub.t.sup.(d))=Σ.sub.iS.sub.it(p.sub.t.sup.(d)), then the utility company's problem for each time step t is given by (step 104b) (equation 10):

(58) $\begin{array}{cc}\underset{p_t^{\left(d\right)},{\tilde{Q}}_t^{\left(n\right)}}{\max }\underset{C_t,Q_t^{\left(n\right)}}{\min }\underset{i}{.Math.}\left[b^f{q}_{\mathrm{it}}^{\left(d\right)}+b_t^s{q}_{\mathrm{it}}^{\left(n\right)}\right]-\left(p_t^{\left(d\right)}+b^f\right)Y_t-C_t\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)},Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right)s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall i\\ Y_t=S_t\left(p_t^{\left(d\right)}\right),& \\ C_t\in C_t,& \\ Q_t^{\left(n\right)}\in Q_t^{\left(n\right)},& \\ y_{\mathrm{it}}\ge 0,& \forall i\\ y_{\mathrm{it}}\le q_{\mathrm{it}}^{\left(d\right)},& \forall i.\end{array}& \left(10\right)\end{array}$

(59) If S.sub.t(p.sub.t.sup.(d)) is known then again this can be directly solved, although there is no assurance that the problem is convex any more since S.sub.t(p.sub.t.sup.(d)) might not be convex. Knowing S.sub.t(p.sub.t.sup.(d)) is a limiting assumption since it implies that either all the user loss functions D.sub.it(y.sub.it) are known (which has important privacy concerns) or that S.sub.t(p.sub.t.sup.(d)) can be recovered by sending the price information to the users' smart meters (which might be a data intensive operation but doesn't have to be done often).

(60) The second step 110 for billing and shift calculation uses the results from the previous step 102, i.e., p.sup.(d) and y.sub.i, ∀i and computes the billing structure B.sub.t, and the demand schedule for the shiftable load x.sub.ij, ∀i, j. Since p.sup.(d) and y.sub.i, ∀i, are known, the utility company's problem is simplified to the following (step 112):

(61) 0$\begin{array}{cc}\underset{X,Y,{\tilde{Q}}_t^{\left(n\right)}}{\max }\underset{C_t,Q_i^{\left(n\right)}}{\min }\underset{t}{.Math.}\underset{i}{.Math.}{B}_t\left(q_{\mathrm{it}}^{\left(s\right)}+q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}},q_{\mathrm{it}}^{\left(n\right)}\right)-C_t\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)},Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right)s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\forall i\\ \underset{t}{.Math.}{x}_{\mathrm{ijt}}={\hat{x}}_{\mathrm{ij}},& \forall j,\forall i\\ x_{\mathrm{ijt}}=0,& \forall j,\forall i,\forall t.Math.\left[{\alpha }_{\mathrm{ij}},{\beta }_{\mathrm{ij}}\right]\\ C_t\in C_t,& \forall t,\\ Q_t^{\left(n\right)}\in Q_t^{\left(n\right)},& \forall t,\\ x_{\mathrm{ijt}}\ge {\underset{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\\ x_{\mathrm{ijt}}\le {\stackrel{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\end{array}& \left(11\right)\end{array}$

(62) There are two main approaches to solve the problems as the one posed in equation 11. One is to assume a utility function for each user that will relate the energy consumption with the price of consuming that energy, and thus control the billing. The problem with this approach in practice is that it is very hard to recover or even estimate these utility functions. The other relies on the fact that since the electric market is highly regulated, in many cases the billing costs to final users are regulated as well and to be kept within certain values. The method 100 utilizes the second approach and assumes that the utility ratio Υ (symbolized as the Greek letter upsilon with hooks and as seen in Equations 12, 14 and 15 below) is fixed. The latter is similar to the article by A.-H. Mohsenian-Rad, et al., but unlike this previous work, the method 100 keeps the ratio separately for the shiftable plus sheddable load and the non-moveable load. Hence, the formulation is as follows:

(63) $\begin{array}{cc}{\rm Y}\equiv \frac{\underset{t}{.Math.}\underset{i}{.Math.}{B}^f\left(q_{\mathrm{it}}^{\left(s\right)}+q_{\mathrm{it}}^{\left(d\right)}-y_{\mathrm{it}}\right)}{\Psi \underset{t}{.Math.}{C}_t^d\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}\right)}=\frac{\underset{i}{.Math.}{B}_t^s\left(q_{\mathrm{it}}^{\left(n\right)}\right)}{{\Phi }_i{C}_t^d\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}\right)+C_t^s\left(Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right)}\ge 1,& \left(12\right)\end{array}$
where the second equality is for all t and

(64) $\begin{array}{cc}\Psi =\frac{\underset{t}{.Math.}{Q}_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}}{\underset{t}{.Math.}{Q}_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+Q_t^{\left(n\right)}};{\Phi }_i=\frac{{\tilde{Q}}_i^{\left(n\right)}}{\underset{i}{.Math.}{Q}_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}},\forall t.& \left(13\right)\end{array}$
Using equations 12 and 13, respective billing structures may be formulated as follows (step 114a):

(65) $\begin{array}{cc}{b}^f=\frac{{\rm Y}\underset{t}{.Math.}{C}_t^d\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}\right)}{\underset{t}{.Math.}{Q}_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}},\mathrm{and}& \left(14\right)\\ b_t^s=\left(\frac{{\rm Y}{\tilde{Q}}_t^{\left(n\right)}}{Q_t^{\left(n\right)}}\right)\frac{C_t^d\left(Q_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}\right)}{\underset{t}{.Math.}{Q}_t^{\left(d\right)}-Y_t+Q_t^{\left(s\right)}+Q_t^{\left(n\right)}}+\left(\frac{{\rm Y}}{Q_t^{\left(n\right)}}\right)C_t^s\left(Q_t^{\left(n\right)}-{\tilde{Q}}_t^{\left(n\right)}\right),\forall t.& \left(15\right)\end{array}$

(66) Note that the value of b.sup.f can be computed without knowing the optimal solution for X, since only Σ.sub.t Σ.sub.i Σ.sub.jx.sub.ijt is needed, which is known and constant. On the other hand, b.sub.t.sup.S cannot be computed until the realization of the non-moveable demand; but, since X is not affected by b.sub.t.sup.S one can compute the optimal shiftable demand by solving the following (step 114b):

(67) $\begin{array}{cc}\underset{X,Q_t^{\left(n\right)}}{\max }\underset{C_i,Q_t^{\left(n\right)}}{\min }\underset{t}{.Math.}\underset{i}{.Math.}{b}^f{q}_{\mathrm{it}}^{\left(s\right)}-\underset{t}{.Math.}{C}_t^d\left(Q_i^{\left(d\right)}-Y_i+Q_t^{\left(s\right)}+{\tilde{Q}}_t^{\left(n\right)}\right)-\underset{t}{.Math.}{C}_t^s\left(Q_t^{\left(n\right)}-Q_t^{\left(n\right)}\right)s.t.\text{:}\begin{array}{cc}{q}_{\mathrm{it}}^{\left(s\right)}=\underset{j}{.Math.}{x}_{\mathrm{ijt}},& \forall t,\forall i\\ \underset{t}{.Math.}{x}_{\mathrm{ijt}}={\hat{x}}_{\mathrm{ij}},& \forall j,\forall i\\ x_{\mathrm{ijt}}=0,& \forall j,\forall i,\forall t.Math.\left[{\alpha }_{\mathrm{ij}},{\beta }_{\mathrm{ij}}\right]\\ C_t\in C_t,& \forall t,\\ Q_t^{\left(n\right)}\in Q_t^{\left(n\right)},& \forall t,\\ x_{\mathrm{ijt}}\ge {\underset{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\\ x_{\mathrm{ijt}}\le {\stackrel{\_}{x}}_{\mathrm{ij}},& \forall t,\forall j,\forall i\end{array}& \left(16\right)\end{array}$

(68) With this problem solved, the method 100 returns to the price incentive and shedding step 102 and repeats all of the method 100 steps until convergence is achieved (step 130). Importantly, the method 100 solves this problem in a centralized manner (unlike the distributed manner presented by the article by A.-H. Mohsenian-Rad, et al.).

(69) From a procurement standpoint, the method 100 considers a two level market from which the utility company purchased the energy required to cover the demand. Second, unlike previous work, the method 100 takes into account that demand level and the procurement prices are indeed unknown stochastic processes. The method 100 accommodates these factors through robust optimization. Finally, the method 100 simultaneously uses demand shedding and demand shifting as a means to shape the demand according to the required objectives. Although three different problems were identified that can be tackled in the framework of the present invention, the method 100 focuses on the revenue optimization problem which is the most important one for utility companies.

(70) In simulations using randomly-generated data, the results of using the control method 100 have been compared against the results of using no control method. The test procedure is described below.

(71) In the simulations, synthetic data was randomly generated for each user/customer. To get useful and more realistic demand patterns, household consumption data was used as a seed to generate the synthetic data for the users/customers. Eighteen (18) months of hourly data for six (6) different users was used to compute seeds that would later generate random instances. In order to define the types of seeds, the data was statistically analyzed in several different ways. FIGS. 3a and 3b show two different examples. In FIG. 3a, the hourly consumption data was grouped for a single user for every month and its distribution was computed, and are presented as vertical histograms, one for each month. Similarly, FIG. 3b shows the result of grouping the data of a single user for a single day of the week (Tuesday in this case) for each hour of the day, and are presented as vertical histograms. It is noted there is a notable change in the consumption pattern depending on the hour of the day, which also varies depending on the day of the week.

(72) By grouping demand by day of the week and hour, a good balance was obtained between recovering the details of the demand distribution and having enough data points to have a meaningful empirical distribution since adding month, for example, would leave only 4 or 5 data points for each group which cannot give enough information.

(73) As FIG. 3c shows, it is also noted that there are differences among the different users in the raw data, although it is difficult to be certain given the low amount of samples. Still, it was assumed that each user represents the average or typical distribution of a certain cluster of users, and this was used as a seed to generate a population of users that has a mix between the different clusters.

(74) Using this information as a seed, the hourly demand pattern for any single day of the week can be generated for a population that contains a mix of the six (6) different users identified in the raw data. This can be done by generating uniformly distributed random variables in the following manner. Let f.sub.u,.sub.d,t(z) denote the empirical probability mass function of the demand z of users in cluster
u ε {1, . . . , 6},
on day of the week d={1, . . . , 7}, and time
t ε {0, . . . , 23},
that were computed from the raw data, and let F.sub.u,.sub.d,t(z) denote its cumulative distribution function. Then, given a day of the week D for each user/customer of cluster U, twenty-four (24) uniformly distributed random variables z.sub.0, . . . , z.sub.23 are generated and the demand of that user i at time t is set equal to q.sub.it=F.sub.U,D,t.sup.−1(z.sub.t), making q.sub.it a random variable with probability mass function f.sub.u,.sub.d,t(z).

(75) Although this method is simple to implement and allows the generation of over a million users very quickly, it has the drawback that the demand of the users will be independent, which is not the case in real life since there is a positive correlation of the energy demand among users. The process to generate correlated random variables is more complicated and requires the definition of a variance-covariance matrix for the users, limiting the number of final users that can be simulated due to the physical memory requirements in the test procedure. The current implementation can generate up to 5,000 users/customers, being limited mainly by the simulation computer's RAM. In this case, when generating an instance, a symmetric, positive definite matrix is randomly created that has only positive correlation factors among variables, S for each cluster within the population that is being simulated.

(76) It is important to note that the final generated random variables (that is, q.sub.it) will not necessarily have S as its variance-covariance matrix since F.sub.u,.sub.d,t(z) could not preserve the correlation, but it serves as a simple and good enough proxy to positively correlate the simulated demand patterns for the users of a cluster. Also, it is noted that due to the correlation, the distribution of the demand of all users of a cluster u in a certain day of the week d and hour t will not follow the corresponding distribution f.sub.u,.sub.d,t. But, as shown in FIG. 3d, if several samples of the same user for the same day of the week and hour are taken, the final distribution for that user (bottom chart) will be identical to the empirical one (top chart), just as expected.

(77) FIG. 3e shows the simulated average user demand for each cluster and for a whole day. Finally, for each user, there needs to be a simulation of what part of that demand is shiftable, sheddable, or non-moveable at every time in the day simulation. Since those patterns might differ from one user to another, only the relative weights of each type of demand were set for each cluster of users and the proportion for each final user, using the relative weights as a seed, were generated.

(78) FIG. 3f shows the total aggregated demand for cluster 4 indicating which part is considered shiftable, sheddable, or non-moveable. To show that each user within this cluster actually has its own randomly generated demand, FIG. 3g shows the daily demand pattern of two different users within the users of cluster 4 showing that although they are different they do have some correlation; whereas, the proportions of shiftable, sheddable, and non-moveable demands are completely random among users.

(79) In the simulations, another important variable to test the method 100 is the energy procurement cost. The simulations used data from the Electric Reliability Council of Texas website (http://www.ercot.com), specifically the Houston prices, to simulate the DAM and RTM prices. Two months of price data were used and the average for the corresponding day of the week in which the simulation is done was used. The DAM prices are given hourly, but given that the RTM prices are provided with 15 minutes resolution, the maximum price for that hour was used as the hourly data to capture the worst case scenarios. FIG. 3h shows an example of the DAM and RTM prices for a Friday.

(80) In the simulations, the method 100 output is compared with a base case in order to analyze how good is the solution given by the method 100. As a base case, the simulations consider the setting in which the utility company has no control over the demand of its users and will just purchase energy in the DAM using a forecast of the total energy consumption for the next day and then buy in the RTM whatever extra energy it requires over what was forecasted. Since the performance will depend on how good the forecast is, for a first set of simulations, it was assumed the non-moveable demand and the RTM prices were known at the moment the problem was solved. In other words, the robustness part was eliminated by setting the RTM prices uncertainty set and the Q.sub.t.sup.(n) uncertainty sets equal to the actual values of those quantities. For this case, FIG. 3i shows the total value of one of the randomly generated demands. Although in the base case the utility company cannot control the final user's demand, the different types of demand (shiftable, sheddable, and non-moveable) are highlighted to compare it with the method 100 output.

(81) Given this demand profile, the base case will procure all its energy on the DAM as shown in FIG. 3j. Since no forecasts are involved there is no over-purchase in the DAM.

(82) Finally, the utility ratio Υ=1.2 is set for all the simulations to keep the same income ratio in the utility company. Then, the costs and billing rates for the base case are as follows:

(83) Output of Base Case

(84) Total Procurement Cost: $ 603087.81 Par: 2.6089 Flat Rate to Consumers: $ 38.54/unit

(85) The performance of the method 100 in simulation is now described. The method 100 is started with b.sup.f and b.sub.t.sup.S equal to the billing rate computed in the base case and the initial value for X.sub.t equal to the shiftable demand Q.sub.t.sup.(s). FIG. 3k shows the result of five (5) iterations of the method 100. When compared to FIG. 3i, it is clear that the method 100 moves the shiftable demand to the hours where the energy is cheaper and, when convenient, it also sheds some demand to reduce the energy procurement. FIG. 3l shows how the energy procurement is made. Since the method 100 can decide when to purchase part of the non-moveable demand, the procurement occurs in the RTM when the prices are lower than the DAM.

(86) Finally, the following is the output of the method 100 in which there are several important things to note. First, the method 100 converges fairly quickly, in just five (5) iterations. Second, the equilibrium results in a solution that is more than 20% cheaper than the base case in terms of procurement costs while at the same time it reduces the billing rates to consumers for the sheddable and shiftable loads, as well as for most of the non-moveable demand, (except when the procurement prices are higher). This last benefit is achieved thanks to the important reduction in procurement costs. Finally, it is noted the PAR is also slightly reduced for this simulation. ioThis might not be the case every time, and it depends on how steep are the changes in prices and the amount of available shiftable demand. In any case, this at least shows that if the PAR is included in the objective it could result in a much better solution in terms of PAR.

(87) Output of Method 100

(88) Step 1—Total Cost: $ 471647.63 Step 2—Total Cost: $ 477478.68 Step 3—Total Cost: $ 477454.26 Step 4—Total Cost: $ 477454.09 Step 5—Total Cost: $ 477454.08 Total Procurement Cost: $ 477454.08 Improvement: 20.8317% New Par: 2.2083 Flat Rate to Consumers (bf): $ 32.47/unit Spot Rate to Consumers (bs): 24.8385 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 31.6698 46.0320 31.6698 84.2850 43.1370 39.2160 36.6045 35.0310 31.6698 31.6698 31.6698

(89) The present invention provides a modeling framework and a control method 100 (expressed in part as an optimization algorithm) for demand shaping through load shedding and shifting in a electrical smart grid. The present invention tackles the problem of controlling the flow of energy, for example, through demand shaping, in a much more general, and novel, framework than the current state of the art. First, unlike current mechanisms, the present invention considers two tools simultaneously for shaping the demand: shedding and shifting of load. Through direct load control the utility company will schedule/shift part of the demand and through price incentives it controls the level of load shedding. Additionally, the present invention considers the two-level market in which utility companies purchase their energy, the day-ahead market and the real-time market, which have different cost structures, and thus addresses the problem of how much energy to purchase on each to cover the demand. Finally, the different optimization sub-problems that must be solved are all posed in a robust optimization framework, allowing the present invention to accommodate the fact that both the final demand and the costs in the real-time market are unknown stochastic processes and not the deterministic values assumed by most of the literature. The final objective of the present invention is then to find the best shedding and shifting profile and how much energy to purchase in each market to satisfy the demand while minimizing costs.

(90) Other modifications are possible within the scope of the invention. For example, the smart power network 10 may be a power sub-network connected to a larger utility power network or even a private power network. Also, the power network 10 and its components have been described in a simplified fashion and may each be constructed in various well-known manners and using various well-known components. Also, the control method 100 may be utilized by energy generators as well as energy distributers since energy generators also buy energy from the energy markets but they also generate it themselves with completely different cost structures.

(91) Also, although the steps of the control method 100 have been described in a specific sequence, the order of the steps may be re-ordered in part or in whole and the steps may be modified, supplemented, or omitted as appropriate. Also, the power network 10 and the computer system 24 may use various well known algorithms and software applications to implement the steps and substeps. Further, the control method 100 may be implemented in a variety of algorithms and software applications.