M.Sc Degree project
An evaluation of a new Pricing
technique to integrate Wind
energy using two Time scales
scheduling
Author: Mutaz Tuffaha Date: 2012-02-08
An evaluation of a new Pricing technique to integrate Wind
energy using two Time scales scheduling
Mutaz Tuffaha February 8, 2012
Abstract
1
Introduction
Energy is one of the most fundamental needs of the human being. The electrical networks were designed and implemented to deliver the energy from the power plants to the end users. At the be-ginning, the researchers and engineers were so focused on increasing the efficiency and reliability of these networks. Then, it was so clear that depending on fossil fuels only in producing energy was- and still is- the biggest threat to our planet and environment. Besides, these fossil fuels are not everlasting. So, the attention now in the research and treatises has been drawn to study the exploitation of the renewable resources of energy such as sun and wind. However, these kinds of energy are intermittent, unpredictable, and they have very low density.
In addition to the above, the power grids have been suffering from the problem of peak hours. The demands of the domestic end users on this network are not constant; they change stochas-tically according to daytime, season, nature of the consumer, and many other factors. In fact, it can be noticed that the demands on the network increase drastically in certain hours of the day, especially in the evening where the people are usually in their houses using various appliances, lights, heaters or AC’s, etc. Such hours of the day are called peak-hours. To tackle this problem, power plants usually use stand by generators working on gas. These stand by generators should not take too long time to operate -unlike the ordinary generators and that is why they are costly.
So, how can we integrate the renewable sources of energy like wind and sun in the power net-work given that these sources are of stochastic nature? How can we increase the reliability and efficiency of such networks? How can we deal with the problem of peak hours? How can we decrease the dependence on the fossil fuels as much as possible?
The researchers, engineers and economists spurred by these questions came up with the idea of smart grids. That is to say, an intelligent network which integrates all possible kinds of energy resources, provides some kind of interaction between the end user and the energy producer via a two-way communication channel and tries to give some solutions to the problems mentioned above.
Of course, one can study these smart grids from different perspectives. One of the most im-portant and recent approaches is the so-called Demand Response or Dynamic response. Demand response is a technique used to reduce the demands in peak hours by raising the price of energy. Thus, the consumer is encouraged to shift his demands away from the peak-hours to off-peak hours. It is believed that this technique would work better if the distribution and procurement of energy are managed by a separate unit (company). Such a unit is called in some articles: Energy Management System(EMS) [1] while others call it Independent System Operator (ISO) [2]. This technique can be formulated as a pricing problem which is, basically, an optimization problem. Considering the stochastic nature of the renewable sources of energy and the uncertainties incurred by the end user demands make this optimization problem very complicated and worthy of study-ing. I have to mention here that to my best knowledge, I have not found any other article that addresses this very problem.
In this thesis, I present the model suggested in [1] with some detailed proofs of the main re-sults. Then, I demonstrate my simulations of this model, the experiments I carried out, and their results. Then, I explain the defect I encountered in this model together with the suggested solution. Finally, I discuss the applicability of this technique in energy markets by pointing out some of its advantages and drawbacks.
In the next section I present this model with some detailed proofs of the main results. In sec-tion three, I explain different methods to generate sample paths or time series of wind energy and traditional users’ demands to pave the way for the simulations. In section four, I present and an-alyze the results of my experiments including the defect I found. In the last section, I draw some conclusions about the advantages and drawbacks of this model.
2
Mathematical Models
M. He et al. in their article [1], suggested a complicated mathematical model of dynamic response that can be used to integrate the wind energy in the power grid. Their model is based on combining the real-time with day-ahead schedules. In this section I try to explain some terminology that is commonly used in the arena of smart grids. I also present this model with all its components. In addition, I try to explain and motivate the mathematical proofs of the solutions suggested by filling in some detailed steps.
2.1 Main elements of the model
Before going through the details of the mathematical model some technical terminology must be explained. In energy markets, one usually think of two different time scales; Day-Ahead schedul-ing in which prices are declared one day ahead based on predictions of the wind energy and user demands. The second is the Real-time scheduling which is "used to tune the balance between the energy amount scheduled day-ahead and the real-time load, i.e., within minutes" [1].
The demand response technique explained in the introduction necessitates the use of smart ap-pliances or meters which, as the name suggests, enable two-way communication between the consumer and the EMS. In simple words, these meters or appliances can receive the prices of en-ergy and feed back the demands to the EMS.
The authors in [1] divided the energy resources into two main categories: The conventional ther-mal energy and wind energy. The wind energy is of stochastic nature, and it is assumed to have zero running cost. The power plants that exploit thermal energy regardless of the source (fossil fuel or nuclear), in turn is divided into base-load and peaking generators. The former represents the ordinary and slow-starting generators with cost c1per unit. These generators have high start-up
cost cp. Hence; these base-load generators are assumed to be working on the day-ahead schedule.
While the peaking generators are fast-starting generators with cost c2per unit that is higher than
the cost of the base-load generators, i.e., c2 is greater than c1. These peaking or stand by
gener-ators are used upon request (in emergency), hence they must start up quickly, and they should be scheduled on the real-time scale. Usually, they work on gas.
To combine the two schedules (day-ahead and real-time) the authors in [1] suggested two time scales. The day is divided into M slots (e.g. 24) each of length T1. Then each T1-slot is divided
into K slots (e.g. 60) each of length T2.
In their model, they also classified the users under two main categories: The traditional user who does not use smart meters or appliances and the opportunistic user who uses them. The tradi-tional user receives a price u per unit corresponding to each T1-slot one day ahead, and this is the
day-ahead pricing. Whereas, the opportunistic user receives a price v per unit at the beginning of each T2-slot. According to their reaction to the retail real-time prices v the opportunistic users are
Figure 1: A flow chart demonstrating the main components of the model.
who leaves the market if he does not accept the price. Both prices (u and v) have upper and lower bounds. So,[1]:
u∈ [umin, ucap], and v ∈ [vmin, vcap]. (1)
Fig.1 depicts the main components of this model.
Before going any further, I list herein the nomenclature used by the authors in [1]:
D= Dt+ Do: The aggregate energy demand from both traditional Dt and opportunistic Do
users in a T2-slot.
W: The wind energy generated in a T2-slot.
L= D −W : The net demand for conventional thermal energy in a T2-slot.
S: The base-load energy generation scheduled by the EMS for a T1-slot of the next day.
s=KS: The base-load generation scheduled for each T2-slot within the T1-slot.
u: The day-ahead retail price corresponding to a T1-slot.
v: The real-time retail price corresponding to a T2-slot.
c1: The cost per unit of the base-load generation.
c2: The cost per unit of the peaking (stand-by) generators.
cp: The start-up cost of the base-load generation.
Let us define the quantity Y as:
Y = s +W − Dt. (2)
Then, the model is supposed to work in the way explained in the flow chart in fig.2. If L is less than s, or equivalently Y is greater than Dobecause:
L< s =⇒ D−W < s =⇒ Dt+ Do−W < s
=⇒ Y > Do,
then the scheduled generation is larger than the demands, which means that the EMS has ordered a surplus energy. So, he loses the cost of the generation of the excess which is equal to (c1−
Figure 2: A flow chart demonstrating how the EMS is supposed to work.
losses. However, this case is unlikely to take place. Finally, if L is greater than s (or equivalently Y is less than Do) then there is a shortage of energy, which means that the EMS needs to order
(dispatch) an amount of L − s (or Do−Y ) from the peaking (stand by) generators at the high cost
of c2. So, he loses an amount of c2(Do− Y ). Thus, we can see how crucial the decision of the
EMS is, and this makes this problem a tempting optimization problem.
2.2 Wind model
The authors in [1] assumed that the wind energy W can be modeled as a non-stationary Gaussian random process across T1-slots. So, in the kth T2-slot of the mth T1-slot the wind energy is given
by [1]:
Wk,m∼N (θm, σ2), (3)
where θmand σ2represent the mean and variance respectively.
2.3 Demand model
The demand D is divided into Dt and Dowhich represent the demands of the traditional users and
the opportunistic users respectively. The traditional users demand Dt can be predicted accurately
and it is modeled as a r.v which has a mean dependent on the retail price u [1]:
Dt= E[Dt] + εt, (4)
with [1]
E[Dt] = αtuγt, (5)
where
αt: a normalizing constant
γt: the price elasticity at the corresponding T1-slot and it is defined by [1]
γt=dE[D t]/E[Dt] du/u = u E[Dt] dE[Dt] du . (6)
Simply put, the price elasticity represents the ratio of the incremental change of the demand to that of the price. This implies that it is a measure of how much the demands change when the price changes. In other words, it infers how sensitive (elastic) the user is towards the price. Worth mentioning that, the price elasticity is usually negative, i.e., γt is less than zero, since the demand
decreases when the price increases in energy markets.
The opportunistic users, on the other hand, are assumed to arrive according to Poisson process with mean κ1equal to λoT2. So, if the total number of opportunistic users is denoted by N then:
N∼Po(κ1).
The rate λois "assumed to be constant over one T1-slot but it may vary from one T1-slot to another
[1]." So, it can be said that:
E[N] = κ1
Var[N] = κ1. (7)
The opportunistic users are assumed also to arrive independently and randomly. Finally, each opportunistic user is assumed to choose randomly and independently a specific acceptance level Vi, if he receives a real-time price v less than or equal to Vihe will accept, otherwise he will reject.
The total number of active opportunistic users Nain a T2-slot, thus, is [1]:
Na= N
∑
i=1χ{v≤Vi}, (8)
with χ is an indicator function, i.e.:
χ{v≤Vi}=
1 if v≤ Vi
0 otherwise
The authors in [1] suggest that the number of active users Na also follows a Poisson
distribu-tion with mean κ1P(vi). This can be understood by considering the number of active users Na as
a compound Poisson process1. According to the theory, if Y (t) is a compound Poisson process defined as: Y(t) = N
∑
i=1 Di, (9)with N is a Poisson process with rate λ and Diare i.i.d, then:
E[Y (t)] = λ tE[D]
Var[Y (t)] = λ tE[D2]. (10)
Applying the result from equation (10) onto the process Nawe get:
E[Na] = λoT2E[χ{v≤Vi}] = κ1P(v ≤ Vi)
Var[Na] = λoT2E[χ{v≤V2 i}] = κ1P(v ≤ Vi). (11)
Now, assuming that each opportunistic user consumes a constant amount of energy of Eo in a
T2-slot, the total opportunistic demand is given by [1]:
Do= NaEo, (12)
hence, the price elasticity of the opportunistic user γocan be defined, just as equation (6), by [1]:
γo= v E[Do] dE[Do] dv = v E[Na] dE[Na] dv . (13)
The opportunistic users are classified according to the value of γo. "If γois between zero and -1 he
is said to be relatively inelastic, otherwise he is relatively elastic [1]". The solution of the above equation is simply given by:
E[Na] = kvγo, (14)
where k is a constant that can be found from boundary conditions. As a boundary condition let us define the value "vmin as the highest price acceptable to all opportunistic energy users [1]".
Obviously, when v is equal to vmin, E[Na] will be equal to E[N] so:
E[N] = kvγmino ⇐⇒ k = E[N] vγo min , (15) and hence: E[Na] =E[N] vγo min vγo = E[N]α ovγo, (16)
where, αois a normalizing constant given by v−γmino. Now, insert equation (7) in the above equation
and compare with equation (11), we can see that:
E[Na] = κ1αovγo
P(v ≤ Vi) = αovγo
⇒ Var[Na] = κ1P(v ≤ Vi) = κ1αovγo. (17)
Furthermore, if we calculate the mean and variance of Doin equation (12) using the results from
the precedented equation, we get:
E[Do] = κ1αovγoEo
Var[Do] = κ1αovγoEo2. (18)
Finally, knowing that κ1 is large enough one can approximate the demand of the opportunistic
users Doby a Gaussian distribution of mean qo(v), and standard deviation σo(v), i.e., [1]:
Do∼N (qo(v), σo2(v))
qo(v) = κ1αovγoEo
σo2(v) = κ1αovγoEo2, (19)
2.4 Problem formulation and solution
To combine the real-time and day-ahead schedules the authors in [1] suggested the following general problem: P : max π M
∑
m=1 Rum(π), (20)and the real-time retail price v. Furthermore, Ru
m(π) is the total profit in a T1-slot and it is given by
[1]: Rum(π) = K
∑
k=1 Eπψl k,m Rlk,m(ψk,ml , π), (21)where Rlk,mis the net profit in the (k, m)th slot i.e., the mth T1-slot of the kth T2-slot. Recall that
Mand K represent the numbers of the T1-slots and T2-slots respectively. While, ψk,ml represents
the system state in the (k, m)th slot and it varies according to the type of the opportunistic user, persistent or non-persistent.
For simplicity, I consider the case of the non-persistent user only in which ψk,ml is equal to {W, Dt}.
The authors in [1] split the optimization problem in (20) into two problems; one to deal with the real-time schedule and the second to deal with the day-ahead schedule.
2.4.1 Real-time optimization problem
The real-time optimization problem is to find the retail real-time price v that maximizes the profit of the EMS in each T2-slot given the wind energy forecast, energy procurement schedule s and the
price u one day ahead.
The Problem:
This problem can be formulated as [1]: PRT non-pst: max v R l(ψl, s, u, v), (22) with [1] Rl(ψl, s, u, v) = uDt+ EvDo[vDo+ χ{A}(−cps) + χ{B}(−cpε − c1(s − ε)) +χ{C}(−c1s+ c2ε )], (23)
where EvDo[.] denotes the expectation with respect to Dogiven v, and
ε = W + s − (Dt+ Do), (24)
furthermore, χ is an indicator of the events A,B and C given by [1]:
χ{A}= 1 if W≥ Dt+ Do 0 otherwise , (25) χ{B}= 1 if χ {A} = 0 and ε ≥ 0 0 otherwise (26) and χ{C}= 1 if ε < 0 0 otherwise (27) The Solution:
To solve this problem the following assumption had to be made: "Wind generation is not suffi-cient to meet the total energy demand in the system [1]". Under this assumption χ{A} is zero,
because W is less than Do+ Dt, and hence χ{B}+ χ{C} is equal to one. Thus the solution would
be [1]:
when the opportunistic user is relatively inelastic, i.e., γois between zero and -1. And [1]: v∗= γo(c1−cp) 1+γo if Y> qo( γo(c1−cp) 1+γo ) γoc2 1+γo if Y< qo( γoc2 1+γo) q−1o (Y ) otherwise , (29)
when the opportunistic user is relatively elastic, i.e., γois less than -1.
The Proof:
Before going through the details of the formal mathematical proof, perhaps it is a good idea to motivate the result shown in equation (29) by explaining why it should look like this. The whole idea behind this result is the price elasticity γo. If the users are relatively inelastic (insensitive),
i.e., γois greater than -1, the users do not care for the price since they need the energy urgently.
So, the EMS can raise the price to its maximum value vcapto ensure maximum profit.
While, when the users are elastic (sensitive), the EMS has to lower the prices to encourage them to consume energy, this explains the function γo
1+γo which is an increasing function as will be shown
later and it approaches unity as γoapproaches −∞. Hence, the coefficient of this function
repre-sents the minimum value reached. Keeping in mind the flow chart shown in fig.2, we can notice that when there is surplus energy (Y is greater than Do) the EMS loses an amount at the cost of
c1− cp, so this is the minimum price he can declare to avoid losses. Whereas, when there is a
shortage of energy (Y is less than Do), the EMS loses an amount at the cost of c2, so this is the
minimum price he can declare to avoid losses.
Now, let us show the mathematical proof in details. The authors in [1] suggested that equation (23) can be simplified to:
Rl(ψl, s, u, v) = uDt+ EvDo[(v − c2)Do− cχ{B}(Y − Do)] − c1s+ c2Y, (30)
I merely try to fill in the detailed steps. Obviously, equation (23) can be simplified using the condition mentioned above to:
Rl(ψl, s, u, v) = uDt+ EvDo[vDo+ χ{B}((c1− cp)ε − c1s) + χ{C}(−c1s+ c2ε )], (31)
which can be rewritten as:
Rl(ψl, s, u, v) = uDt+ EvDo[vDo+ χ{B}((c1− cp)ε − c1s− c2ε + c2ε ) + χ{C}(−c1s+ c2ε )]. (32)
Define:
c= c2+ cp− c1, (33)
and notice from equations (24) and (2) that ε is equal to Y − Do, the last equation can be written
as:
Rl(ψl, s, u, v) = uDt+ EvDo[vDo+ χ{B}(−c(Y − Do) − c1s+ c2(Y − Do)) + χ{C}(−c1s+ c2(Y − Do))]
= uDt+ EvDo[vDo+ χ{B}(−c(Y − Do)) − c1s(χ{B}+ χ{C}) + c2(Y − Do)(χ{B}+ χ{C})]
= uDt+ EvDo[(v − c2)Do− cχ{B}(Y − Do)] − c1s+ c2Y, (34)
where in the second and third lines of the above, the facts that χ{B}+ χ{C}is equal to one under
the assumption mentioned earlier, EvDo[Y ] is equal to Y , and E
v
Do[s] is equal to s were used.
distributed with mean qo(v), and standard deviation σo(v), which means that the density function
fDo(x) of the r.v. Docan be written as:
fDo(x) = 1 σo √ 2πe −(x−qo(v))2 2σ 2o ,
where, x denotes the real value taken by the r.v. Do instead of do to resolve any confusion that
may appear in integrations later. Having said this, the expectation EDvo in equation (34) can be calculated as follows:
EvDo[(v − c2)Do− cχ{B}(Y − Do)] = (v − c2)qo(v) − cE
v
Do[χ{B}(Y − Do)]. (35)
Keeping in mind that χ{B}is equal to unity when ε is greater than or equal to zero (i.e., Y − Dois
greater than or equal to zero, or Dois less than or equal to Y from equations (24) and (2)) the inner
expectation can be calculated as follows:
EvDo[χ{B}(Y − Do)] = Z Y −∞(Y − Do)dFDo(x) = Z Y −∞ (Y − x) σo √ 2πe −(x−qo(v))2 2σ 2o dx = Y Z Y −∞ 1 σo √ 2πe −(x−qo(v))2 2σ 2o dx− Z Y −∞ x σo √ 2πe −(x−qo(v))2 2σ 2o dx. (36)
The first integral in the above equation corresponds to the CDF FDo(Y ). Letting y be equal to
Y−qo(v)
σo , and letting Φ denote the CDF of the normalized normal distribution, FDo(Y ) can be written
as Φ(y). Whereas, the second integral can be determined by the substitution of u as (x−qo(v))2
2σ2 o to get: Z Y −∞ x σo √ 2πe −(x−qo(v))2 2σ 2o dx = Z Y −∞ x− qo(v) σo √ 2π e −(x−qo(v))2 2σ 2o dx+ Z Y −∞ qo(v) σo √ 2πe −(x−qo(v))2 2σ 2o dx = ( u = (x−qo(v))2 2σ2 o du dx = x−qo(v) σo2 ) = Z y2/2 ∞ σo √ 2πe −udx+Z Y −∞ qo(v) σo √ 2πe −(x−qo(v))2 2σ 2o dx = −√σo 2πe −u u=y2/2 u=∞ + qo(v)Φ(y) = −√σo 2πe −y2/2 + qo(v)Φ(y). (37)
Inserting the result from the above equation in equations (36), then in (35) and (34) we finally get:
Rl(ψl, s, u, v) = uDt− c1s+ c2Y+ (v − c2)qo(v) − cσo √ 2πe −y2/2 − c(Y − qo(v))Φ(y). (38)
The final stage to prove the solutions given in equations (28) and (29) is to find the derivative of the profit function Rl(ψl, s, u, v) in equation (38). However, the authors in [1] distinguish between
two cases: The case when the opportunistic user is relatively inelastic (−1 ≤ γo< 0) and the case
when the opportunistic user is relatively elastic (γo< −1). The idea behind this classification is to
approximate the profit Rl(ψl, s, u, v) with a less complicated function. This approximation depends
on observing the fact that σo2(v) is less than σ2+ σε2t where σ
2and σ2
εt represent the variances of
that of the first case. In a nutshell, the authors in [1] proved that the profit Rl(ψl, s, u, v) can be
approximated by [1]: ˜
Rl(ψl, s, u, v) = (v − c2)qo(v) + uDt− c1s+ c2Y (39)
when γois between zero and (-1). And [1]:
˜
Rl(ψl, s, u, v) = (v − c2)qo(v) + uDt− c1s+ c2Y− c(Y − qo(v))+ (40)
when γois less than -1.
The case when −1 ≤ γo< 0:
When γo is between zero and -1, the profit ˜Rl(ψl, s, u, v) is a linear increasing function of the
retail real-time price v. Since it was assumed from the beginning that the prices have a cap de-noted by vcap, the maximum of the profit function is attained at this value.
The case when γo< −1:
When γois less than -1, let us write the profit ˜Rl(ψl, s, u, v) as:
˜
Rl(ψl, s, u, v) = (v − c2)qo(v) + uDt− c1s+ c2Y− c(Y − qo(v))χ{Y −qo(v)≥0} (41)
to account for the positive part of Y − qo(v). Then, using the result of qo(v) from equation (19) the
derivative is obtained to be:
∂ ˜Rl ∂ v (ψ
l, s, u, v) =
κ1αovγoEo+ (v − c2)κ1αoγovγo−1Eo
+ c(κ1αoγovγo−1Eoχ{Y −qo(v)≥0}− (Y − qo(v))δ{Y =qo(v)}), (42)
where χ and δ corresponds to the indicator and Dirac function respectively. Denote the value of v at which the profit attains its maximum by v∗. Then, to find the zeros of ∂ ˜Rl
∂ v we consider the three
cases below:
I. When (Y − qo(v∗) > 0): The indicator χ{Y −qo(v)≥0}is equal to unity, hence and after dividing
by κ1αoEowe get:
v∗γo+ (v∗− c
2)γov∗γo−1+ cγov∗γo−1= 0
⇐⇒ v∗γo−1(v∗+ v∗
γo− c2γo+ cγo) = 0.
Obviously, the solution v∗= 0 does not count because it corresponds to the minimum. So: v∗+ v∗γo− c2γo+ cγo= 0
⇐⇒ v∗(1 + γ
o) = (c2− c)γo= (c1− cp)γo
⇐⇒ v∗=(c1−cp)γo
(1+γo) .
II. When (Y − qo(v∗) < 0): The indicator χ{Y −qo(v)≥0}is zero, hence we get:
v∗+ v∗γo− c2γo= 0
⇐⇒ v∗(1 + γo) = c2γo
⇐⇒ v∗= (c2)γo
(1+γo).
III. When (Y − qo(v∗) = 0): In this case the Dirac function dominates so:
2.4.2 Day-Ahead optimization problem
While the other optimization problem deals with the day-ahead schedule. So, the target is to find the day-ahead policy or decision (i.e., day-ahead price u and energy procurement schedule S of each T1-slot next day) that maximizes the profit of the EMS.
The Problem:
This day-ahead problem can be formulated as [1]:
PDA non-pst: max S,u K
∑
k=1 EWkE u Dtkmax vk [Rl(ψkl, s, u, vk)], (43)where the inner maximum corresponds to the real-time optimization problem in (22). The target is to obtain the day-ahead policy that maximizes the profit over one T1-slot which consists of K
T2-slots, that explains the sum. However, the maximum of this problem can also be achieved by
considering the maximum for one T2-slot when the opportunistic users are non-persistent, so the
above problem can be reformulated as:
PDA non-pst: max S,u EWE u Dt[R l (ψl, s, u, v)]. (44) The Solution:
The policy of the day-ahead optimization problem in (43) ,that contains the decision of S and u, was found to be [1]: u∗ = ucap if − 1 ≤ γt< 0 γt 1+γtc1 if γt < −1 S∗ = arg max S {(c2− c1)S + KE v W,D[(v − c2)Do −cχ{B}(S − K(αtu∗γt+ εt+ Do−W ))]}, (45)
which will be proved later to be specifically [1]:
S∗= K(κ1Eoαov∗γo+ αtu∗γt− FZ−1(cp/c)), (46)
where FZ−1(.) denotes the inverse of the CDF of Z which is equal to W − εt.
The Proof:
Surely, I can motivate the results in equation (45) in the same way I motivated the results of the day-ahead problem, but I do not want to repeat myself, so let us start with the proof straight away.
The authors in [1] define:
for simplicity. Thus, Y in equation (2) would be equal to s0+W −εtand the function EWEuDt[R
l(ψl, s0, u, v)]
can be rewritten using equation (34)as:
EWEuDt[R l(ψl, s0, u, v)] = EWEuDt[uDt+ E v Do[(v − c2)Do− cχ{B}(Y − Do)] − c1s+ c2Y] = EWEuDt[uDt+ E v Do[(v − c2)Do− cχ{B}(s 0 +W − εt− Do)] − c1s +c1EuDt[Dt] − c1E u Dt[Dt] + c2(s 0 +W − εt)] = (u − c1)EuDt[Dt] + EWE v Do[(v − c2)Do− cχ{B}(s 0+W − ε t− Do)] +(c2− c1)s0+ c2EW[W ] = (u − c1)αtuγt+ EWEvDo[(v − c2)Do− cχ{B}(s 0 +W − εt− Do)] +(c2− c1)s0+ c2EW[W ]. (48)
Intuitively, the expectation above can be split into three separate functions, as suggested by the authors [1], namely f1(u), f2(s0) and c2EW[W ] with
f1(u) = (u − c1)αtuγt
f2(s0) = EWEvDo[(v − c2)Do− cχ{B}(s
0+W − ε
t− Do)] + (c2− c1)s0, (49)
so, the maximum of the expectation in equation (48) can be obtained by maximizing each function in equation (49) separately.
I. maxuf1(u): The derivative of f1(u) is obtained to be:
d f1(u)
du = αt(u
γt+ γ
tuγt−1(u − c1)).
Denote the value of u at which the profit attains its maximum by u∗, solve for the zeros of
d f1(u)
du and divide by αtwe get:
u∗γt+ γ
tu∗γt−1(u∗− c1) = 0
⇐⇒ u∗γt−1(u∗+ γ
tu∗− γtc1) = 0.
Needless to say that the trivial solution (u∗= 0) has to be excluded, hence:
u∗(1 + γt) = γtc1⇐⇒ u∗=
γtc1
1 + γt
.
However, we can not accept this solution when γt is between zero and -1, because the price
would be negative. Instead, the cap value of the price ucapwould be adopted.
II. maxs0 f2(s0): Actually, by finding the function f2(s0) the result shown in equation (45) is
proved keeping in mind the fact that S is equal to Ks and the definition of s0 from equa-tion (47). All I need to prove now is the specific result of S∗given in equation (46).
First, define Z by W − εt, then ε in the event B which corresponds to the case when ε which
is given by W + s − (Dt+ Do) is greater than or equal to zero, can be rewritten using this
definition and equations (47), (4) and (5) as:
ε = W + s − (Dt+ Do) = W + s − αtuγt− εt− Do
So, the event B corresponds to the case when Z + s0− Dois greater than or equal to zero. To
find the maximum of f2(s0) we need to find its derivative, but we can say:
d f2(s0) ds0 = d ds0(EWE v Do[(v − c2)Do− cχ{B}(s 0+W − ε t− Do)] + (c2− c1)s0) = EWEvDo[ d ds0((v − c2)Do− cχ{Do≤Z+s0}(s 0 + Z − Do))] + (c2− c1) = EWEvDo[−cχ{Do≤Z+s0}− c(s 0+ Z − D o)δ (s0+ Z − Do)] + (c2− c1) = −cEWEvDo[χ{Do≤Z+s0}] − cEWEDo[(Z + s 0− D o)δ (s0+ Z − Do)] +(c2− c1). (50)
Of course, the Dirac function makes all point zero except when s0+ Z − Do is zero, so the
second term in the last line of the precedent equation would be zero. While the first term is more complicated. Bear in mind that the r.v. Do is normally distributed with mean qo(v),
and variance σo2(v). Remember also that Z is equal to W − εt, hence it is a r.v. normally
distributed as:
Z∼N (θm, σε2t+ σ2).
Given that, σ2is much larger than σ2
o(v) + σε2t, the r.v. Do− Z can be approximated by:
(Do− Z) ∼N ((qo(v) − θm), σ2).
Hence, I suggest proceeding as follows:
EWEvDo[χ{Do≤Z+s0}] = EWE v Do[χ{Do−Z≤s0}] = Z s0 −∞ dFDo−Z(x) = Z s0 −∞ fDo−Z(x)dx = Z s0 −∞ 1 σ √ 2πe −(x−(qo(v)−E[Z]))2 2σ 2 dx= u = qo(v) − x du dx = −1 = Z qo(v)−s0 ∞ 1 σ √ 2πe −(u−E[Z])2 2σ 2 (−du) = Z ∞ qo(v)−s0 1 σ √ 2πe −(u−E[Z])2 2σ 2 du = 1 − Z qo(v)−s0 −∞ 1 σ √ 2πe −(u−E[Z])2 2σ 2 du = 1 − PZ(Z ≤ qo(v) − s0), (51)
where, in the second line I use the substitution shown. Insert the result of the last equation in equation (50) and solve for d f2(s0)
From the definition of s0 in equation (47) and definition of qo(v) in equation (19), and by
taking the values of the retail prices v∗and u∗which maximize the profit, the result in equa-tion (46) is proved.
In summary of this section I presented the mathematical model of pricing and scheduling electric-ity suggested by the authors in [1] to integrate the wind energy. I also tried to motivate and explain the main results by filling in some detailed steps.
3
Generating Sample paths of the next 24 hours
Designing a simulation for such a model is a cumbersome task. Besides, I do not think that any one can simulate the human behavior because it is so fuzzy and unpredictable especially, when it comes to money. However, generating a time series of the wind energy for a certain day can be done if one has information on the wind energy one day-ahead which can be obtained form weather forecast as mentioned earlier. While, according to this model the traditional users’ demands depend on the prices and price elasticity. Hence, a time series can be generated if we know or assume these prices and elasticities. In this section, I describe different methods to generate time series of wind energy and traditional users’ demands to establish a base to carry out some experiments later.
3.1 Wind model
This model depends on two time scales as mentioned earlier. The day was divided into 24 hours and each hour was divided into 60 minutes, so M is equal to 24, and K is equal to 60. The first step in my simulation was to create a stochastic process for the wind. As mentioned above, the wind energy in the (m, k)th slot was modeled as a non-stationary Gaussian random (stochastic) process as shown in equation (3), in which, θmand σ2 represent the mean and variance respectively. The
mean and variance of the wind energy in a T1-slot are assumed to be available one day ahead. It is
assumed that weather forecast centers can predict the wind speed for the next 24 hours with good accuracy, then the energy manufacturer can process this data to predict the wind energy available over the next 24 hours.
Generating a stochastic process is not an easy task. My target is to generate a time series represent-ing the wind energy expected for the next 24 hours. I commenced by assumrepresent-ing M (which is equal to 24) values for the means of the wind energy in each hour of the next day, I called this vector of means Θ. To make it look more realistic I assumed that the wind energy during midday would be less than that in the early morning and at night. My assumption was based on a simple calculation. Assuming that the energy consumption of an ordinary user-whether traditional or opportunistic-is ∼ 0.25KWh in an hour on average. Assuming also that we want to consider a town of 20,000 households, the total demands would be in the range of ∼ 5MWh in an hour. It is not unrealistic to assume that the wind generation can reach up to one-third of this total demand. So, in the units of MWh, let this vector be:
Θ = ( 1.5 1.4 1.4 1.3 1.6 1.6 1.6 1.4 1.2 1.0 1.0 0.9
0.8 1.0 1.1 1.1 1.4 1.5 1.5 1.6 1.7 1.7 1.7 1.6 ). (53)
Actually, I meant by each element of this vector the mean wind energy in a T1-slot (an hour), so the
mean wind energy in each T2-slot (a minute) was obtained by dividing each element by K (=60).
While the variance is a completely different story. In fact, the authors in [1] proposed that the standard deviation of the wind energy is usually in the same order of its mean. Hence, I assumed values for the variances depending on the mean wind energy of each T1-slot as I show later.
matrix would be diagonal (with zero off-diagonal elements), otherwise we must think of the cor-relation coefficients. Of course it looked much easier to consider them uncorrelated. Anyway, I talk about that covariance matrix after explaining the techniques I used to generate this time series.
I used three techniques to generate the wind generation time series. The generation of this time series is not the target of my study here I just want to see the influence of the pricing algorithm explained above on it. Hence, the time series by itself is not of great importance, however I want my simulation to be not far from reality. So, I explain these three techniques in the following subsections.
3.1.1 Random Number Generator
This technique is easy and straight-forward. In fact, MATLAB includes the function randn that generates random numbers normally distributed with mean zero and variance one. One can modify the mean and variance easily in MATLAB. Let X be a one-dimensional r.v. with mean zero and variance one, such that:
X ∼N (0,1). If Y is a one-dimensional r.v. with mean θ and variance σ2,i.e.,
Y ∼N (θ,σ2).
A sample of Y can be easily generated in MATLAB using the following command:
X= randn(.,.) Y = θ + σ ? X .
However, when the r.v. Y is M-dimensional with mean vector Θ and covariance matrix Σ we need to be careful. Let X be an RM-valued normally distributed r.v. with mean zero and covariance
matrix equal to I the identity M × M matrix. Then a sample of Y can be obtained by [4]:
Y = Θ + β X , (54)
where β is a lower triangular matrix obtained from Cholesky decomposition. The idea behind this decomposition can be summarized as follows. If Σ is a symmetric positive definite matrix, a triangular matrix β such that:
Σ = β0β , (55)
can be obtained by the following equation [4]:
βi j=
σi j− ∑r=1j−1βirβjr
q
σj j− ∑r=1j−1β2jr
. (56)
Fortunately, MATLAB contains the function chol which calculates this matrix β without the need to write an algorithm to execute equation (56). In fact, the choice of the covariance matrix Σ plays a crucial role here. If the elements of the process Y are assumed to be uncorrelated then the co-variance matrix Σ will be diagonal as mentioned earlier. Under this assumption, there is no need to use Cholesky decomposition because the matrix β would be a diagonal matrix with diagonal elements equal to the square roots of the variances.
Figure 3: The sample path of the wind process using normal random generator when the covari-ance matrix is diagonal.
to be equal to the mean wind energy in a T2-slot, i.e., in one minute divided by 2 to reduce the
probability of getting negative values for the energy, hence I set σ to θ /(2K). Finally, the stochas-tic process W (t) can be formed by combining M (=24) rows, each of length K elements, in one long time series horizontally called ts_wind. So, the following M-file was written in MATLAB to generate the process and the result is shown in fig.3.
M=24; K=60; theta=[1.5 1.4 1.4 1.3 1.6 1.6 1.6 1.4 1.2 1.0 1.0 0.9 0.8 1.0 1.1 1.1 1.4 1.5 1.5 1.6 1.7 1.7 1.7 1.6]; Beta=diag(theta’/(2∗K)); Theta=repmat(theta’/K,1,K); X=Theta+Beta∗randn(M,K); W=X(1,:); fori=2:M
W=horzcat(W,X(i,:));Combining all the M random vectors horizontally end
ts_wind = timeseries(W(1,:), 1:1440,’name’, ’Wind Energy’); ts_wind.DataInfo.Units = ’MWh’;
ts_wind.TimeInfo.Units = ’minutes’; ts_wind.TimeInfo.StartDate=’00:00’ ; ts_wind.TimeInfo.Format=’HH:MM’ ; figure(1)%The figure depicted in fig. 3 plot(ts_wind)
I have to admit that fig.3 does not look so reasonable because the energy reaches negative val-ues at some moments, but this is a result of the model which assumed that the standard deviation is in the range of the mean!
While, if the elements of Wm,kare correlated then the covariance matrix will not be diagonal. To
decide about that, I think that we have to go back to the origin of the wind energy and electricity generation to understand this process more deeply and then we may know something about the correlation coefficient. Anyway2, I replicated the values of the standard deviations over an M × K matrix, applied the random number generator in equation (54) and multiplied element wise. The details are given in the following M-file and the result is shown in fig. 4.
Figure 4: The sample path of the wind process using normal random generator when the covari-ance matrix is not diagonal.
M=24; K=60; theta=[1.5 1.4 1.4 1.3 1.6 1.6 1.6 1.4 1.2 1.0 1.0 0.9 0.8 1.0 1.1 1.1 1.4 1.5 1.5 1.6 1.7 1.7 1.7 1.6]; Beta=repmat(theta’/(2∗ K),1,K); Theta=repmat(theta’/K,1,K); X=Theta+Beta.∗ randn(M,K); W=X(1,:); for i=2:M W=horzcat(W,X(i,:)); end
ts_wind = timeseries(W(1,:), 1:1440,’name’, ’Wind Energy’); ts_wind.DataInfo.Units = ’MWh’;
ts_wind.TimeInfo.Units = ’minutes’; ts_wind.TimeInfo.StartDate=’00:00’ ; ts_wind.TimeInfo.Format=’HH:MM’ ; figure(2)%The figure depicted in fig. 4 plot(ts_wind)
Comparing fig.4 with fig.3 it can be noticed that the sample paths are not so different from each other.
3.1.2 Monte Carlo Method using Fourier series
M. Grigoriu in [4] suggested the use of Fourier series to generate a sample of RM-valued non-stationary Gaussian process X (t) by Monte Carlo method. The process X (t) in this algorithm is supposed to have zero mean and covariance function ci, j(t, s) given by [4]:
ci, j(t, s) = E[Xi(t)Xj(s)] for i, j = 1, · · · , M.
He also proved that "if the Fourier series of ci, j converges to the covariance function ci, j of X at
every interior point of the rectangle [0, τ] × [0, τ], then the process X(ni)with coordinates X(ni)
becomes a version of the Gaussian process X (t) as ni→ ∞ for i = 1, · · · , M [4]." Where: X(ni) i (t) = Ai,0 2 + ni
∑
k=1[Ai,kcos(νkt) + Bi,ksin(νkt)], (57)
with Ai,k= 2 τ Z τ 0
Xi(u) cos(νku)du, k= 0, 1, · · · , ni,
Bi,k0 =2
τ Z τ
0
Xi(u) sin(νk0u)du, k0= 1, · · · , ni, (58)
and νk= kν1= 2πτ k. So, his suggested algorithm was [4]:
a Generate samples of the correlated Gaussian random variables (Ai,k, Bi,k0) for i = 1, · · · , M,
k= 0, 1, · · · , niand k0= 1, · · · , ni.
b Calculate X(ni)(t) from (57).
The covariance matrix (Σ) was assumed in this algorithm to be diagonal as justified before. The samples of Ai,kand Bi,k0 in equation (57) had to be generated. It looked intuitive to assume that Ai
is normally distributed with mean zero, and covariance matrix equal to the identity matrix I, where Iis (ni+ 1) × (ni+ 1). The same can be said about (Bi) except that ni+ 1 had to be replace with
ni because k0= 1, · · · , ni where k = 0, 1, · · · , ni as stated in equation (58). I wrote the following
Figure 5: The sample path of the wind process using the Monte Carlo algorithm when ni= 200.
W=X(1,:); for i=2:24
W=horzcat(W,X(i,:)); end
ts_wind = timeseries(W(1,:), 1:1440,’name’, ’Wind Energy’); ts_wind.DataInfo.Units = ’MWh’;
ts_wind.TimeInfo.Units = ’minutes’; ts_wind.TimeInfo.StartDate=’00:00’ ; ts_wind.TimeInfo.Format=’HH:MM’ ; figure(3)%The figure depicted in fig.5 plot(ts_wind)
Again, in the above algorithm I did not have to use the chol function to obtain β , I simply assumed it to be a diagonal matrix of the standard deviations which are equal to the mean in a T2-slot divided by 2 as done in the previous algorithm. To make the sample more consistent with
results obtained by the first algorithm I had to divide the β matrix by the variance of the sample Xsin the M-file because the sample was made of standard normally distributed random numbers
Ai,k0 and Bi,k multiplied by a series of sine and cosine functions, so the sample Xs did not have
unity variance necessarily. Hence, I created this sample y which contained the reciprocals of the standard deviations of Xs. The result of the above M-file is shown in fig. 5 below.
3.1.3 Modified Brownian Motion
The other method I found to generate a sample of the wind process is the function wfbm. It is a function in the wavelet toolbox in MATLAB designed to generate a fractional Brownian motion depending on the parameter H called Hurst parameter with:
0 < H < 1.
This H parameter controls the sample according to the following relation [5]:
Var(X (t) − X (s)) = σ |t − s|2H, (59)
Figure 6: The sample path of the wind process using the modified BM when ni= 200.
increments will be positively correlated. The following M-file was written and the result is shown in fig. 6. M=24; K=60; Ni=200; theta=[1.5 1.4 1.4 1.3 1.6 1.6 1.6 1.4 1.2 1.0 1.0 0.9 0.8 1.0 1.1 1.1 1.4 1.5 1.5 1.6 1.7 1.7 1.7 1.6]; A=zeros(M,Ni); for i=1:M H =0.1; GAi=wfbm(H,Ni); theta_Ai=theta(1,i)/K*ones(1,Ni); A(i,:)=theta_Ai+(theta(1,i)/(2*K))*GAi/var(GAi); end W=A(1,:); fori=2:24 W=horzcat(W,A(i,:)); end
ts_wind = timeseries(W(1,:), 1:length(W),’name’, ’Wind Energy’); ts_wind.DataInfo.Units = ’MWh’;
ts_wind.TimeInfo.Units = ’minutes’; ts_wind.TimeInfo.StartDate=’00:00’ ; ts_wind.TimeInfo.Format=’HH:MM’ ; figure(4)%The figure depicted in fig.6 plot(ts_wind)
From fig.6 it can be noticed that the modified Brownian motion is far from the results of the other techniques regardless of how much the parameter H can be changed. Anyway, I chose the second technique in generating the wind process (shown in fig.4) throughout my simulations due to its simplicity and straight-forwardness.
3.2 Traditional user and Day-Ahead price
The traditional energy users’ demand Dt model explained in equations (4) and (5) was easier to
(4). Let us remember here that this uncertainty can be considered as a white noise, i.e.,:
εt∼N (0,σε2t).
Therefore, it was clear that a sample of this uncertainty should be generated as the wind process sample was generated in precedented subsection. Thus, εt ∈ RMwith each random vector εtm
con-sisting of K elements. Then, the process of εt was formed by combining the M random vectors of
εti horizontally. I assumed the variance σ
2
εt to be less than (let us say one fifth) the variance of the
wind energy σ2 which was assumed to be equal to half of the mean wind energy in each T2-slot,
hence I set σε2t to the mean of the vector Θ divided by 10K. The following M-file was used in MATLAB to generate this sample:
epsilon_t=mean(theta)/(10*k); X_epsilon=epsilon_t*randn(M,K); epsilon_t=X_epsilon(1,:);
for i=2:M
epsilon_t=horzcat(epsilon,X_epsilon(i,:)); end
As mentioned earlier, the traditional energy users’ demand can be expected with good accuracy depending on the day-ahead price as shown in equation (5). Hence, I had to find the optimal retail day-ahead price from equation (45). Obviously, the value of optimal day-ahead price u∗depends mainly on γt and the other parameters. Firstly, I had to assume values for the other parameters as
follows:
ucap = 100
c1 = 60 c2 = 75
cp = 40, (60)
which can be thought of as prices ine/MWh (or as percentages) with the condition: cp< c1< c2.
However, I have to mention here that one of the most ambiguous issues here is the maximum price ucap, and the normalizing constant αt. The authors in [1] did not suggest explicitly a method to
choose these values. But they state that the price elasticity of the traditional energy user γt may
vary from one T1-slot to another and they assumed a value of −0.5 for γt in their experiments.
So, to create a situation closer to reality I assumed a vector of values of γt for the 24 T1-slots as
follows:
γt = ( −1.1 −1.09 −1.1 −1.09 −1.1 −1.0 −0.9 −0.8 −0.7 −0.8 −0.7 −0.6
−0.5 −0.5 −0.3 −0.2 −0.2 −0.3 −0.4 −0.3 −0.4 −0.6 −1.15 −1.1 ).
My strategy in assuming these values was based on the fact that in early morning hours the de-mands are low so the users are sensitive (elastic) to prices, while in midday hours they are less sensitive because they need the electricity more during this period.
The normalizing constant αt was a little bit trickier. In equation (53), I assumed the means of
the wind generation amount for 24 hours, the mean traditional energy users’ demand should be comparable to those means. Remembering that εt was assumed to be white noise of mean zero, I
had to choose a value of αt that could increase the mean of Dt and keep it insensitive to the price
elasticity γt. At the same time I did not want the mean of Dt to be always close to unity. So, I
Figure 7: The sample path of the traditional users’ demands of the next day.
Actually, the elasticity αt may vary from one T1-slot to another as well, but I considered this
vari-ation dependent only on the elasticity of the user.
Finally, the value of u∗ was calculated from equation (45). The following M-file was written in MATLAB and the results are depicted in fig.7.
u_cap=100; c1=60; c2=75; cp=40; u_opt=zeros(1,M); gamma_t=[-1.1 -1.09 -1.1 -1.09 -1.1 -1.0 -0.9 -0.8 -0.7 -0.8 -0.7 -0.6 -0.5 -0.5 -0.3 -0.2 -0.2 -0.3 -0.4 -0.3 -0.4 -0.6 -1.15 -1.1]; alpha_t=zeros(1,M); for m=1:M
if gamma_t(1,m)<0 && gamma_t(1,m)>=-1 u_opt(1,m)=u_cap; alpha_t(1,m)=0.1*(1/c1)^gamma_t(1,m); elseif gamma_t(1,m)<-1 u_opt(1,m)=gamma_t(1,m)*c1/(1+gamma_t(1,m)); alpha_t(1,m)=(1.5/c1)gamma_t(1,m); end end U_opt=alpha_t(1,1)*(u_opt(1,1))^gamma_t(1,1)*ones(1,60); for m=2:M U_opt=horzcat(U_opt,alpha_t(1,m)*(u_opt(1,m))^gamma_t(1,m)*ones(1,60)); end D_t=U_opt+epsilon;
ts_Dt = timeseries(D_t(1,:), 1:1440,’name’, ’Traditional user Demand’); ts_Dt.DataInfo.Units = ’MWh’;
ts_Dt.TimeInfo.Units = ’minutes’; ts_Dt.TimeInfo.StartDate=’00:00’; ts_Dt.TimeInfo.Format=’HH:MM’ ; figure(7)%The figure depicted in fig.7 plot(ts_Dt)
adequate to describe the demands of domestic traditional users in a normal working day in winter. I tried to make the demands increase in the peak hours which spread over afternoon and evening in normal days.
In this section, I explained different methods to generate the stochastic process or time series of the wind energy for a whole day according to the model suggested. in addition I generated a time series of the traditional users’ demands depending on the expected prices and price elasticities.
4
Experiments and Evaluation
Having generated the time series of the wind energy and traditional users’ demands in the previous section, several experiments and measurements can be carried out now. One of the most interesting measurements that can be found is the profit. In this section, I calculate the scheduled energy, real-time prices and profits under several circumstances. Then I try to clarify the influence of the price elasticity and wind energy on the profits to eventually study the dynamics of this model. In addition, I explain the defect of this model which I found during my experiments. Then, I suggest a solution for this defect.
4.1 Scheduled Energy, Opportunistic users and Profits
From equation (45) it can be noticed that the scheduled dispatch issued by the (EMC) depends on:
a) The expected traditional users’ demand which depends on the day-ahead price u and its elasticity γt and we have this from the precedent section.
b) The expected wind energy Wt and we have it from the time series of the wind energy in the
previous section.
c) The expected opportunistic users’s demand which is given by κ1αov∗γoEo and we do not
have this yet.
Due to this expected opportunistic users’ demand calculating the scheduled dispatch was a tedious task. On the one hand, I could assume values like γo= −2 and hence v∗would be equal to vcap. On
the other hand, the opportunistic price elasticity γovaries from one T1-slot to another. In addition,
letting γobe constant would make the expected opportunistic users’ demand constant as well, and
thus, the model would not be realistic. Another problem is that the optimal scheduled energy S∗ must be dispatched one day-ahead and it depends on the optimal real-time price v∗, which in turn, depends on the elasticity γt. So, if one wants to determine the scheduled energy from 02:00 till
03:00 on December 21st, for example, he would use the price elasticity γo(and hence the price v∗)
from 02:00 till 03:00 on December 20th.
Anyway, I decided to use the same value of the price elasticity of the traditional users γt for the
price elasticity of the opportunistic user γt because my main target was to see the influence of the
wind energy on this model. I set the other parameters to the following values:
Eo= 1 × 10−3(MWh)
λo= 50, T2= 1(min) ⇐⇒ κ1= 50
vmin= c2⇐⇒ αo= (1/c2)γo
Then, FZ−1(.) had to be determined and that was easy. Recall that:
Then, the r.v’s W and εt both were assumed to be normally distributed, and they both ∈ RM.
Assuming they are independent, the r.v. Z would be also normally distributed with mean equal to the difference between their means and variance equal to the sum of their variances, i.e.,:
Z∼N (θZ, σZ2),
with:
θZ= θm− θεt
σZ2= σ2+ σε2t. (62)
Fortunately, MATLAB has the function norminv(· · · ) which finds the value of the normally dis-tributed r.v. corresponding to a certain probability. To find the r.v. Z, we need to remember the standardization of random variables. If z denotes the standard r.v. then:
z= Z− θZ σZ
. (63)
Taking the value of z to be norminv(cP/c)and using equation (62) we get:
Z= norminv(cP/c)
q
σ2+ σε2t+ θm− θεt. (64)
Of course, θεt was assumed to be zero, and I had the variance which was called epsilon_t in
the previous M-file.
Having all the elements, I could calculate the optimal scheduled energy S∗ from equation (45) using the value of γo0 for a certain T1-slot of the next day. Then, as mentioned above, the value of
γowas used to determine the optimal real-time price v∗from equation (29) for each T2-slot in that
T1-slot.
The last step was to determine the profits of the EMS in that specific T1-slot. Simply put, the profit
can be calculated according to this equation:
R= Gt+ Go−C, (65)
where, Gt, Goand C denote the gross income from the traditional users in a T1-slot, gross income
from the opportunistic users in a T1-slot and the aggregate cost in a T1-slot respectively, and they
are given by:
Gt= K
∑
k=1 u∗kDtk, (66) Go= K∑
k=1 v∗kDok (67) and C= K∑
k=1 c1s∗ if Wk+ s∗− (Dtk+ Dok) ≥ 0 c1s∗+ c2(Dtk+ Dok− (Wk+ s ∗)) if W k+ s∗− (Dtk+ Dok) < 0 , (68)which can be rewritten using equation (2) as:
C= K
∑
k=1 c1s∗ if Yk− Dok≥ 0 c1s∗+ c2(Dtk−Yk) if Yk− Dok< 0 . (69)Actually, this equation says that when there is surplus energy (Wk+ s∗− (Dtk+ Dok) ≥ 0) the
loss is due to the cost of the scheduled energy, whereas when there is a shortage of energy (Wk+ s∗− (Dtk+ Dok) < 0) the EMS has to procure the difference from the stand by (fast start-up
and costly) generators at a cost of c2 per unit energy. The following M-file was written in
lambda_o=50; kappa_1=lambda_o*1; E_o=e-3; v_cap=100; x=3; gamma_op=-1.1; gamma_o=-1.1; S_opt(1,m)=K*(kappa_1*E_o*(v_cap/c2)^(-2)+alpha_t(1,x)*(u_opt(1,x))^gamma_t(1,x) -(norminv(cp/c)*sqrt(Beta(x,1)^2+epsilon_t^2)+theta(1,x)/K)); end Y=s_opt+W(1,(x-1)*60+1:x*60)-D_t(1,(x-1)*60+1:x*60); v_opt=zeros(1,K); alpha_o=(1/c2)^gamma_o; v1=gamma_o*(c1-cp)/(1+gamma_o); v2=gamma_o*c2/(1+gamma_o); q_o=kappa_1*alpha_o*E_o; sigma_o=sqrt(kappa_1*alpha_o*E_o^2); for k=1:K
if gamma_o<0 && gamma_o>=-1 v_opt(1,k)=v_cap; elseif gamma_o<-1 if Y(1,k)>=q_o*v1^gamma_o v_opt(1,k)=v1; elseif Y(1,k)<q_o*v2^gamma_o v_opt(1,k)=v2; else v_opt(1,k)=(Y(1,k)/q_o)^(1/gamma_o); end end end D_o=q_o*v_opt.^gamma_o+sigma_o*sqrt(v_opt.^gamma_o).*randn(1,60); G_t=u_opt(1,x)*sum(D_t(1,(x-1)*60+1:x*60)); G_o=v_opt*D_o’; C=0; for k=1:K if (Y(1,k)-D_o(1,k)) >= 0 C=C+c1*s_opt(1,k); else C=C+c1*s_opt(1,k)+c2*(D_o(1,k)-Y(1,k)); end end R=G_t+G_o-C;
In the above M-file, x denotes the number of the T1-slots during which the calculations will be
made. While, gamma_op and gamma_o denote γo0 and γorespectively.
4.2 Results and Analysis
I wanted first to see the effect of wind energy on the dynamics of this model. To ensure fair comparison I kept the price elasticities of the traditional and opportunistic users constant and I increased the mean wind generation by a step of 0.2MWh and recorded the scheduled energy and the profits in a T1-slot (an hour). I carried out this experiment for two different times of the day:
1) 02:00-03:00 during which the users can be considered elastic (or sensitive) to prices because it can be considered as an off-peak hour.
2) 19:00-20:00 during which the users can be considered inelastic (or insensitive) to prices because it can be considered as a peak hour.
2:00-3:00 (γo= −1.1, γt= −1.1) 19:00-20:00 (γo= −0.3, γt = −0.3) W(MWh) S∗(MWh) R(e) W(MWh) S∗(MWh) R(e) 0.5 2.999 1807.1 0.5 7.2329 333.4 0.55 2.9360 1827.4 0.55 7.1690 338.3 0.6 2.8719 1828.3 0.6 7.1049 330.7 0.7 2.7433 1833.1 0.7 6.9764 346.7 0.8 2.6143 1845.8 0.8 6.8475 352.0 0.9 2.4850 1859.1 0.9 6.7183 358.5 1.0 2.3556 1845.1 1.0 6.5889 367.0 1.1 2.2260 1858.5 1.1 6.4593 372.8 1.3 1.9665 1862.1 1.3 6.1999 376.1 1.5 1.7068 1856.9 1.5 5.9402 390.6 1.7 1.4469 1903.2 1.7 5.6804 397.6 1.9 1.1869 1903.5 1.9 5.4204 411.5 2.1 0.9268 1908.5 2.1 5.1603 419.0
Table 1: The effect of increasing wind energy on the profits during two one-hour intervals corre-sponding to off-peak hour and peak hour respectively.
4.2.1 Price volatility
This problem of price volatility was not clear to me from the beginning. To understand the reason of this phenomenon we have to look at the price function given in equation (45). During the time interval from 02:00-03:00 the price elasticity was assumed to be -1.1. When the price elasticity γt
is less than -1, the optimal price u∗is given by c11+γγtt which is depicted in fig.8 below.
Obviously the price function has a horizontal asymptote at u = ucap and this means that as γt
approaches −∞ the price u∗ approaches c1. That is logical, because when the user is extremely
sensitive to the price, the EMS has to lower the price as much as possible and the lowest price allowed uminwas assumed to be c1. The price function has also a vertical asymptote at γt= −1. In
the region where γtis larger than -1 we have no problem as the optimal price is set to the maximum
price ucap. While in the region where γtis close to -1 from the left we have a problem, because the
price u will exceed tremendously the maximum price ucap. For example, if γt is equal to -1.1, the
price will be
u∗= c1
−1.1
1 − 1.1 = 11c1.
As a suggestion to solve this problem, I found the price elasticity corresponding to the value of the maximum price ucap, and I denote it here by γtcap as shown in fig.8:
ucap= c1 γtcap
1+γtcap
⇐⇒ ucap+ γtcapucap= c1γtcap
⇐⇒ γtcap(c1− ucap) = ucap
⇐⇒ γtcap =
ucap
c1−ucap
So, I suggest to modify the optimal price function to the following:
u∗ = ucap if γt≥ ucap c1−ucap γt 1+γtc1 if γt< ucap c1−ucap (70)
Figure 8: The behavior of the price function u = c11+γγt t.
4.2.2 The influence of wind generation on the profits
Going back to the results shown in table 1 we can notice, as a rule of thumb, that the profits increase when the wind generation increases. Which makes sense since the wind generation was assumed to be cost free. We wanted to describe this relation quantitatively by considering the rate of change of the profits with respect to the wind generation. However, finding the rate of change of the profits with respect to wind generation ∆R
∆W is not easy due to the stochasticity of the wind
energy generation and opportunistic users’ demands. Instead, I suggested an empirical formula to describe the incremental change of the profits. This formula was based on the results presented in table 1 and on observing the equations (65) to (69). So, assuming that the wind energy in each T2-slot changes by ∆θm,k, I suggested the following formula for the incremental change of the
profits: ∆R = K
∑
k=1 v∗∆Dok+ c1norminv(cp/c) K∑
k=1 ∆σZ,k+ c1 K∑
k=1 ∆θm,k, (71)because, the gross income Gt of the traditional users would not change by the change of wind
energy except for some change due to the uncertainty εt. While, the gross income Goof the
oppor-tunistic users are supposed to be more random since it is set each T2-slot, and since it was modeled
as a stochastic process. The last two terms in the above equation come from the definition of the FZ−1and the r.v. Z in equation (64).
To test this theory, I considered the values obtained for the time interval 19:00-20:00 since in this interval the price elasticity γo is greater than -1, i.e., the optimal price v∗ would be equal to
vcapfor all T2-slots. So, the first term of the above equation can be rewritten as: K
∑
k=1 v∗∆Dok= v ∗ K∑
k=1 ∆Dok, (72)while in the time interval 02:00-03:00 the price elasticity γo is less than -1, i.e., the optimal price
v∗would change from one T2-slot to another, thus the calculations would be more difficult. Then,
I neglected the change of the gross income of the traditional users’ demands Gt for the reason
About ∆σZ, remember that I assumed in my model σ to be2Kθm, and σεt to be θm 10K so: ∆σZ= K q σ22+ σε2 t,2− K q σ12+ σε2 t,1 = K r (θm,2 2K ) 2+ (θm,2 10K) 2− K r (θm,1 2K ) 2+ (θm,1 10K) 2 = Kθm,2 r ( 1 2K) 2+ ( 1 10K) 2− Kθ m,1 r ( 1 2K) 2+ ( 1 10K) 2 = K∆θm r ( 1 120) 2+ ( 1 600) 2≈ 0.51∆θ m. (73)
Finally, recalling that I assumed θm,kto be θm/K in my model of the wind process, equation (71)
can be simplified to:
∆R ≈ v∗
K
∑
k=1∆Dok+ 0.0085c1norminv(cp/c)∆θm+ c1∆θm. (74)
I used MATLAB to find ∑Kk=1∆Dok of a T1-slot for the values of the mean wind energy given in
table 1 and I calculated the change of profits from table 1 and using equation (74) and the results are shown in table 2 below.
In table 2, I considered the value of wind mean energy of 0.5MWh as a reference, so all the
19:00-20:00 (γo= −0.3, γt= −0.3) ∆θm(MWh) ∑Kk=1∆Dok(MW h) v∗∑Kk=1∆Dok+ c1∆θm(0.51norminv(cp/c) + 1) ∆R(e) 0.05 -0.066 -2.7 4.9 0.1 -0.1635 -8.55 -2.7 0.2 -0.0442 11.2 13.3 0.3 -0.0159 21.8 18.6 0.4 -0.0202 29.2 25.1 0.5 -0.0510 33.9 33.6 0.6 -0.0581 41.0 39.4 0.8 -0.0427 58.1 42.7 1.0 -0.0206 75.9 57.2 1.2 -0.1185 81.7 64.2 1.4 -0.0178 107.4 78.1 1.6 -0.0492 119.9 85.6
Table 2: Calculation of the incremental profits.
changes are calculated from that value. It can be noticed from the results in table 2 that the incre-ment of the profits ∆R found in MATLAB is close to the value calculated according to equation (74) which makes my theory work except, of course, for some difference which I think is due to the uncertainty of the traditional users’ demands εtsince I did not consider that change in equation
(74).
Anyway, this result is very significant since it shows the advantages of this model. Firstly, we can say that the profits increase with the wind generation, with a rate proportional to the cost of the base-load generators c1. Actually, that is logical because the more wind energy we have, the
less energy we need to dispatch from the energy production plants. Secondly, the incremental profits depend also on the change of the opportunistic users’ demands and i think this is the main merit of this model. Although the change of the opportunistic users’ demands Doare negative in
(a) (b)
Figure 9: The time series of the incremental change of the profits ∆Rkin the time interval
19:00-20:00 when the mean wind energy changes by (a) 0.5MWh, and (b) 0.05MWh.
can be achieved if we can increase γo, i.e., if we can make the opportunistic users more insensitive
to the prices when the wind energy increases and this is the idea of the whole model in my opinion.
Out of curiosity, we wanted to see how this incremental change of the profits varies in each T2-slot.
So, instead of adding up the changes of the profits ∆R, I found the change of the profits ∆Rk in
each T2-slot then I put them in one time series. The results are shown in fig.9. which confirms our
main result that the profits increase with the wind energy generation.
This can be seen because the curve in fig.9(a), which corresponds to a change of 0.5MWh of wind energy, is higher than the curve in fig.9(b) that corresponds to a change of 0.05MWh only. Besides, we can notice the randomness of this signal, which emphasizes the stochasticity of the incremental change of the profits. Which, in turn, explains the difficulty of finding the rate of the change of the profits with respect to the wind.
In summary of this section, I used my simulation of the wind process and the traditional users’ demand time series to carry out simple experiments of this model. In addition, I elucidated and performed the calculations of the scheduled energy, real-time price and profits under different cir-cumstances. Then, I showed that this model suffers from the defect of the price volatility when the price elasticity γo,t lies in the proximity of -1. I also suggested a solution for that defect by the
al-ternative equation given in (70). Furthermore, I showed how the profits of the EMS increase when the wind energy increases and when the demands of the opportunistic users increase. Finally, I tried to describe this change of the profits quantitatively by equation (74).
5
Conclusion
I tried to elucidate this model with its proofs by filling in some detailed steps in the second sec-tion. Then in the third section I generated a time series of the wind process and demands of the traditional users by different methods to provide a good material for my experiments. While, In the fourth section, I carried out simple experiments, in which I increased the mean wind energy in small steps, to study the influence of the wind generation on the profits using this model. During my experiments, I found that the prices will increase drastically when the price elasticity is in the proximity of -1. I explained the reason of this defect and I suggested a solution for it. Finally, I showed that, with this model, the profits increase with good rates, not only when the wind gener-ation increases, but also when the users’ demands increase.
In conclusion, I can say that this model can serve as a rich field for research and study in smart grids. The reason for that conclusion is that this model, to my best knowledge, is the first to dis-cuss integrating wind energy in energy markets. On the one hand, from the results of my simple experiment, I noticed the defect of the price volatility, and I suggested a solution for this problem by modifying the model given in equation (45) to equation (70). On the other hand, I noticed that the profits of the EMS in this model would increase if the wind generation increases and if the price elasticity increases. In other words, if the opportunistic users are encouraged strongly to use the energy more when the wind generation increases, the profits would increase too.
At the end, I would like to say that energy companies like NORDPOOL use in general the day-ahead schedule. Actually, we held a useful meeting with BIXIA in their premises in Växjö. In that meeting they explained to us how they buy the energy from NORDPOOL and how they offer their bids to sell extra energy. All those contracts are made one day-ahead. While the model sug-gested in [1] uses two timescales: day-ahead and real-time. I think that the model deserves more studying. Perhaps for future work it is a good idea to analyze some data from energy companies like BIXIA and see the influence of this model on this data.
6
Acknowledgment
I would like to thank my supervisor Röger Pettersson for his continuous support and encour-agement, his enriching comments and thoughtful discussions definitely made my research easier. Then, I need to thank my examiner Börje Nilsson for directing me to put the final touches to this thesis. I would like also to express my gratitude to the staff of DFM in Linnæus university-Växjö, especially, Astrid Hilbert for introducing me to the world of stochastic processes.
References
[1] M. He, S. Murugesan and J. Zhang. Multiple Timescale Dispatch and Scheduling for Stochas-tic Reliability in Smart Grids with Wind Generation Integration. arXiv:1008.3932v2 [cs.SY], 13 Sep 2010.
[2] M. Roozbehani, M. Dahleh and S. Mitter. Dynamic Pricing and Stabilization of Supply and Demand in Modern Electric Power Grids. IEEE 978-1-4244-6511-8, 2005.
[3] S. Caron and G. Kesidis. Incentive-based Energy Consumption Scheduling Algorithms for the Smart Grid. IEEE 978-1-4244-6511-8, 2010.
[4] M. Grigoriu. Stochastic Calculus: Applications in Science and Engineering. Boston: Birkhäuser, 2002.
SE-391 82 Kalmar / SE-351 95 Växjö Tel +46 (0)772-28 80 00
