- Email: [email protected]
Modeling real-time balancing power demands in wind power systems using stochastic differential equations
Electric Power Systems Research 80 (2010) 966–974
Contents lists available at ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
Modeling real-time balancing power demands in wind power systems using stochastic differential equations Magnus Olsson ∗ , Magnus Perninge, Lennart Söder Royal Institute of Technology (KTH), Electric Power Systems, Teknikringen 33, 10044 Stockholm, Sweden
a r t i c l e
i n f o
Article history: Received 7 August 2008 Received in revised form 18 February 2009 Accepted 18 February 2009 Available online 2 February 2010 Keywords: Stochastic differential equations Real-time balancing Wind power
a b s t r a c t The inclusion of wind power into power systems has a significant impact on the demand for real-time balancing power due to the stochastic nature of wind power production. The overall aim of this paper is to present probabilistic models of the impact of large-scale integration of wind power on the continuous demand in MW for real-time balancing power. This is important not only for system operators, but also for producers and consumers since they in most systems through various market solutions provide balancing power. Since there can occur situations where the wind power variations cancel out other types of deviations in the system, models on an hourly basis are not sufficient. Therefore the developed model is in continuous time and is based on stochastic differential equations (SDE). The model can be used within an analytical framework or in Monte Carlo simulations. © 2010 Elsevier B.V. All rights reserved.
1. Algebraic symbols Variables in continuous time are denoted with the time argument with parenthesis, as e.g. P(t), while discrete time is indicated with subscripts, as e.g. Pk . Stochastic variables are in uppercase Latin letters, e.g. Pk . Realizations or historical data of the same are in lowercase letters, as e.g. pk . The indicator function I¯(A) ∈ {0, 1} take the value I¯(A) = 1 if the logical expression A is true, otherwise I¯(A) = 0. Below follows a list of the variables and most important parameters. Additional parameters are defined on first occasion. D(t) P(t) L(t) Q (t) W (t) N(t) (t) ak , dk X(t) ˛, ω(t) U(t) ,
demand for balancing power [MW] base production [MW] base load [MW] controllable production [MW] uncontrollable production [MW] net base load [MW] mean function of L(t) parameters of (t) stochastic process in L(t) parameters of X(t) deterministic part of W (t) stochastic process in W (t) parameters of U(t)
∗ Corresponding author. E-mail addresses: [email protected] (M. Olsson), [email protected] (M. Perninge), [email protected] (L. Söder). 0378-7796/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.epsr.2010.01.004
(t) Y (t) rk (tk ) ˇ Z(t) (Tn , Vn )
deterministic part of Q (t) stochastic process in Q (t) segments of (t) parameter in Y (t) stochastic jump process in Y (t) stochastic pairs defining Z(t)
2. Introduction This paper addresses the problem of balance management in power systems. The aim is to describe a new probabilistic model of the time continuous balancing demand in MW in a system with large amounts of wind power. This model can, e.g. be used to estimate energy demands for primary, secondary and tertiary control. In addition to be useful within an analytical framework, the presented model can be applied in Monte Carlo simulations. Performing such analyses and simulations is of great importance not only for system operators, but also for producers (and possibly consumers) since they in most systems provide balancing power through various market solutions. The presented model provides a tool for calculation and simulation of the demand for frequency control and required ramp rates. Probabilistic load or production models are usually based on an hourly time resolution as in, e.g. [1] and [2]. However, to be able to perform precise estimations of the impact of new power sources into the power system on the continuous need for balancing power, models using an hourly time resolution are not sufficient. This is due to the fact that the production may vary significantly during the hour, as e.g. in the case of wind power.
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
The inclusion of wind power into the power system has a significant impact on the need for real-time balancing power due to the stochastic nature of the wind power production. However, the introduced unpredictable wind power production variations will also sometimes help the system in the real-time balancing process. This is in this article denoted as counter balancing. In time discrete models with, e.g. an hourly time resolution, assumptions regarding what happens during the hour become necessary to consider counter balancing. Since this can occur at any time and span a time period of arbitrary length, models in continuous time are more suitable to capture this. Concerning intra-hourly modeling, not much work has been published. In [3], an optimization model for minimization of the real-time balancing cost in systems with uncertain wind power forecasts is presented. The proposed model suggests which bids to call and when this should be performed. The model in [3] and the model presented in this paper address the same problem, but from different perspectives. The model in [3] focuses on optimizing the bids given the load and the scheduled production, while the model presented in this paper represents the continuous demand for balancing power. In [4], a simulation model of the impact of wind power production on operation of the power system is presented. However, the presented model does not have a stochastic representation. In papers [5–8] detailed models of wind turbines or farms are presented for analyses of the interaction between the power system and the wind farms. Such models are in theory possible to use for estimation of the increased demand for real-time balancing power associated with integration of wind power, but simulating a system with such accuracy for a 24 h period, for example, is not possible due to computational time. Generally, these types of simulation models concern time frames of a few seconds in length. Articles concerning the impact of wind power on real-time balancing quantities include [9] and [10]. However, these have a time resolution of one hour and thus do not include counter balancing within the hour. This paper presents stochastic models in continuous time of loads and productions in order to capture the continuous demand for balancing power. Suitable tools for modeling stochastic processes in continuous time are stochastic differential equations (SDEs) [11]. SDEs have previously been used in various applications including biology, physics, signal processing and finance [12]. 3. Background 3.1. Market structure Some of the properties that need to be considered in the modeling process are set by the market rules. These are: • The stipulated trading period length. This will affect the dispatch of power into the system since the producers tend to change the production at the borders of the trading periods where the price changes. • The possibility for actors to trade power according to, e.g. updated forecasts before the start of the actual trading period. Considering trading periods, trading is usually performed on an hourly basis as on the European Energy Exchange (EEX) [13], PJM [14], and in the Nordic system [15]. In some systems, as e.g. UK [16], half-hour periods are instead used. For convenience, this paper assumes a trading period of one hour, but the suggested model can be applied to any period length. In this paper, no specific market structure is assumed. However, one market property needs to be considered: At which time the actors no longer have the possibility to trade power ahead for a
967
specific hour. This needs to be considered in the modeling due to forecast uncertainties. In many power systems, ahead trading can be performed on the following market places: • Day-ahead spot market, where actors can trade electric energy for the various hours of the next day. Examples of day-ahead spot markets are EEX [13], Nord Pool Elspot [17] and Omel [18]. • Intra-day market, which offers possibilities to adjust the traded energy quantities from the day-ahead market according to updated forecasts or to non-accepted bids on the day-ahead market. Here it is possible to continuously trade energy for a specific hour after the day-ahead market has closed, until a few hours before the start of that hour. Examples of intra-day markets are Nord Pool Elbas [17] and APX NL [19]. Hence, in this setting, the interesting market property is when the intra-day market closes for the various hours since this represents the last possibility to trade energy ahead. Apart from the ahead markets, there can also exist a realtime balancing market where producers and consumers have the possibility to trade power with the system operator during the operational hour. This is a part of the frequency control often referred to as tertiary control. Examples of such markets are the real-time balancing market in PJM [14], regulating market in the Nordic system [15], the real-time energy market in New England [20], and frequency control ancillary services (FCAS) market in Australia [21]. A compilation of balance management in Europe can be found in [22]. The model described in this paper does not require any assumption regarding the real-time market. However, the model can, e.g. be applied in estimations or simulations of demands for real-time balancing power. 3.2. Load forecasting Load forecasting encompasses forecasts on several time scales, ranging from hourly to yearly. The models presented in this paper consider a number of hours and include load forecasting on an hourly time scale. This means that load forecasts in MWh/h are made each hour according to updated information for individual future hours. Load forecasting has been subjected to substantial amounts of research. The motivation for development of accurate forecast methods lies in the economic impact of the forecast errors [23]. Short-term load modeling and forecasting is reviewed in [24,25]. Recent work on short-term forecasting includes [26–28]. 3.3. Wind power forecasting The uncertainties associated with wind power originates from uncertainties in wind speed forecasts. This in combination with the rapid development of wind farms generates a need for better wind power production forecasting methods. The higher forecast reliability the lower the balancing costs become in the system, which in case of large-scale integration of wind power can imply substantial savings for the wind farm owners as well as better overall efficiency of the system [29,30]. Different wind power forecast techniques are presented in e.g. [31,32] and the state-of-the-art literature review [33]. 4. Definitions 4.1. Continuous demand for balancing power The overall aim of this paper is to model the continuous demand for balancing power D(t), measured in MW. To be able to achieve
968
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
4.2. Controllable and uncontrollable generation
Fig. 1. Basic system with mechanical production P(t), electrical production P e (t) and load L(t).
this, no dispatch for frequency control is included in the developed models. Hence, dispatch related to primary, secondary or tertiary control is excluded in the model of D(t). The continuous demand for balancing power is now defined as the difference between the mechanical production P(t) and the electrical load L(t): D(t) = P(t) − L(t).
(1)
Mechanical production here refers to the power from the turbines, originating from, e.g. water or steam hitting the blades. The difference between mechanical and electrical production is pointed out in the fundamental system in Fig. 1. Note that if replacing P(t) with the electrical production P e (t) in (1) will give D(t) ≡ 0 since P e (t) ≡ L(t) by the law of conservation of energy! The difference between P(t) and P e (t) is covered by changes in rotational speed of the synchronous machines in the system. The actual mechanical production, including wind power production, unforeseen shut-downs of generation units, etc. but with excluded frequency control is denoted base production throughout this paper. The corresponding electric load is denoted base load. To clarify this terminology, various schematic production and load curves are illustrated in Fig. 2. Starting with line PL in the figure, this line represents the electrical production and load including frequency control. Hence, this is the actual electrical production/consumption in the system. The other lines in the figure represent the following: The solid lines represent the planned mechanical production (black line P1) and forecasted load (grey line L1); the dashed black line P2 corresponds to the mechanical base production including e.g. unforeseen shut-downs of units and variations in mechanical wind power production; and the dashed grey line L2 is the base load. Note that in P2 and L2, neither production nor consumption changes for frequency control are included. Concerning the load, the lines L2 and PL do not coincide, implying that some of the frequency control has been performed on the load side in this figure. In terms of the lines in Fig. 2, the overall aim of this paper is to model the difference between P2 and L2, i.e. base production minus base load. The various production types in this paper consider mechanical productions, but for simplicity, the term “mechanical” is from now on omitted.
The wind power is in this paper considered as uncontrollable, implying that there are no possibilities to control the power output. This assumption can also be applied on the base load. Contrarily, except for the wind power production, the base production is controllable since it is the sum of decisions made by the producers in the power system. Thereby two different types of processes are considered: • Uncontrollable: Base load and wind power production, which are processes independent on forecasts and decisions. • Controllable: Base production with excluded wind power production. This is a result of the decisions of the producers, which are based on forecasts. Note that controllable does not imply that the production is deterministic. For example faults can occur resulting in shut-downs of production units. The uncontrollable base load and wind power production are combined into the net base load N(t), N(t) = L(t) − W (t),
(2)
where L(t) denotes the base load and W (t) the wind power production. Hence, N(t) represents the base load that is not covered by the wind power production. The corresponding for base production can be expressed as Q (t) = P(t) − W (t),
(3)
where Q (t) is the net base production and P(t) is the total base production including wind power. Q (t) thereby represents the hydrothermal base production in the system. The definitions of N(t) in (2) and Q (t) in (3) imply that N(t) is uncontrollable while Q (t) is controllable. Q (t) is dependent on decisions made by the market actors, which in turn are based on the market structure and the available forecasts of the net load. Because of the differences between (L(t), W (t)) and Q (t), the modeling is divided into two main parts: (1) Modeling L(t) and W (t): L(t) and W (t) defines N(t) according to (2). In this paper, L(t) and W (t) are considered as endogenous variables. (2) Modeling Q (t): The aim of the model of Q (t) is to simulate the sum of the actions taken by the producers, resulting in the total base production excluding wind power in the system. Q (t) is dependent on forecasts of the base load and the wind power production. The aim to model the demand of balancing power can be reached by studying the difference between P(t) and L(t) (i.e. the difference between P2 and L2 in Fig. 2). This reflects the continuous demand for balancing power, D(t), and hence D(t) = P(t) − L(t) = Q (t) − N(t).
(4)
The advantage of studying Q (t) and N(t) instead of P(t) and L(t) is the distinction between controllable and uncontrollable processes. 5. Main assumptions The proposed models in this paper assume the following:
Fig. 2. Load and production curves. The black lines correspond to productions, where P1 is the planned production and P2 is the base production P(t). The grey lines correspond to base loads where line L1 represents the forecasted base load and L2 the actual base load L(t). Line PL represents the electrical production and consumption when frequency control is included.
(i) The wind power production is independent of the base load. (ii) The wind power production does not change discontinuously. This assumption can be interpreted as that the wind power plants are 100% reliable and will never be subjected to shutdowns. Another interpretation is that each plant is of arbitrarily
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
small size. This together with the assumption that the transmission lines are 100% reliable, and that the probability of shut-downs of two units at the exact same time is zero, implies that no shut-downs will cause any discontinuities in the total wind power production. (iii) The producers (excluding wind power production) deliver a constant base production output during the hour unless the system operator orders a production change. (iv) The market structure allows for trading of electric energy in MWh/h up to a few hours before the start of the operational hour. Market actors trade power according to updated forecasts of base load and wind power production until this market closes. Hence, when the intra-day market closes the following is valid concerning the hourly quantities: • Wind power producers are in balance considering traded quantity and forecasted production. • Other producers are in balance considering traded quantity and planned production. • Consumers are in balance considering traded quantity and forecasted load. (v) Available data series are hourly base load levels and wind power production levels. Both are measured in MWh/h. 6. Modeling the net base load N(t) In this section, the models of the base load N(t) = L(t) − W (t) ∈ R,
t≥0
(5)
are presented. The ingoing parts L(t) ∈ R+ and W (t) ∈ R+ are, as stated in main assumption (i) in section 5, assumed independent and are presented separately. 6.1. Modeling L(t)
Fig. 3. Example of mean load curve. The staircase-shaped curves are the hourly load levels and the solid line represents the model of (t) given by the DCT.
Now, the process L(t) = (t) + X(t) solves the following SDE: dL(t) =
L(t) = (t) + X(t),
t ≥ 0.
(6)
d(t) dt
+ ˛ ( (t) − L(t) ) dt + dB(t).
(9)
The inclusion of d(t)/dt is necessary in order to make L(t) a meanreverting process [35]. What is now left is the mean function (t), which can be chosen in several ways. Here a model based on the discrete cosine transform (DCT) [36] is used. The available data do, as further described in Section 8, in general not consist of samples of the total base load in MW, but of energy levels during the trading periods in MWh. Thus, the time integral of the base load is modeled, which then can be differentiated to get the load itself. The integral of (t) ≥ 0 is non-decreasing in time and therefore, also a linear term is added to the model. Hence, (t) can be written
The model of the load is briefly described below. A more detailed description can be found in [34]. The base load is modeled as the sum of a function representing the mean of the load and a stochastic process:
969
t
(s)ds = 0
a0 t +
N
ak dk cos
(2k − 1)(t − 1) 2N
(10) ,
k=1
The mean (t) ∈ R+ should resemble the continuous variation of the base load between various trading hours in order to get the typical load characteristics. The stochastic process X(t) ∈ R is chosen according to the following:
where N is the number of samples, a0 is the linear coefficient, and ak and dk are given by the DCT. An example of (t) can be found in Fig. 3 where the staircase-shaped curve corresponds to historical hourly load mean values and the solid line represents the model.
• It seems reasonable to assume, in accordance with the central limit theorem, that L(t) is normally distributed since it corresponds to a sum of various loads in the system. Hence, X(t) should be normally distributed. • The square mean of L(t) cannot, e.g. increase during many successive hours, but should instead return to the mean. This implies that some sort of mean-reversion is necessary.
6.2. Modeling W (t)
Putting the above together, the deviation of the load from its mean function is modeled with an SDE of Ornstein–Uhlenbeck type [12] according to: dX(t) = −˛X(t)dt + dB(t),
˛ ∈ R+ ,
(7)
where (B(t), t ≥ 0) is a Brownian motion. The solution to this SDE can be found by applying Itô’s formula [11] and is given by
X(t) = e−˛t (X(0) +
t
e˛s dB(s)),
(8)
0
where X(0) is the initial condition for the base load deviation from the mean and the integral is taken in the Itô sense.
The following can be noted regarding wind power production: • Wind speeds do not have symmetrical distributions [37]. Thereby, the distribution of W (t) should not be symmetrical. • The square mean of W (t) cannot, e.g. increase during many successive hours. This means that some mean-reversion should be included in the model. • The wind power production W (t) do not show the same strong cyclic behavior as the load L(t) on a daily basis. • The wind power production is bounded on the interval [0, W ], ¯ is the total installed wind power production capacity. where W The wind speed is often considered to be Weibull distributed. However, in this paper the log-normal distribution is applied for W (t) due to the benefit of having normally distributed logarithms. This is an advantage from a mathematical point of view because of the possibility of using well established SDEs based on Brownian motion.
970
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
Based on the properties above and on assumption (ii) in Section 5, the wind power production is modeled as W (t) = eω(t)+U(t) ∈ R+ ,
t ≥ 0,
(11)
where ω(t) ∈ R+ is an algebraic function representing the mean of ln W (t), and U(t) ∈ R is an Ornstein–Uhlenbeck process dU(t) = −U(t)dt + dB(t),
t≥0
(12)
where B(t) is a Brownian motion. Thereby, U(t) is normally distributed and is a diffusion, which can be expressed as
t
U(t) = e−t (U(0) +
es dB(s)).
(13)
0
ω(t) should reflect cyclic patterns of the wind power production. The modeling of ω(t) depends on the wind situation in the actual power system. In this paper, ω(t) is assumed constant and will further on be denoted without time argument, i.e. ω. ¯ , The expression in (11) can result in situations where W (t) > W which is not possible in reality. However, if the wind farms are widely spatially distributed across the system, the probability that ¯ becomes very small. Therefore, this is neglected in the W (t) > W developed model. If applying the model of W (t) within an analytical framework, it might become necessary to handle the case ¯ . In Monte Carlo simulation applications, the sample W (t) > W paths can be truncated to not violate the bounds of W (t). However, ¯ is very small, this is generally not if the probability that W (t) > W a problem in such simulations.
Fig. 4. Principle behavior of (t). The time variable t ∈ [0, ∞), while tk ∈ [0, 1) with tk = 0 when t = k.
traded by the hour. This together with assumption (iii) in Section 5, results in that changes in base production tends to be performed in steps at the border of each hour. The system operators often have the possibility to order the producers to start production earlier, e.g. up to 15 min in front of and/or at the end of hour. This is to smooth out the steps in base production that otherwise will occur in the system at the start of each hour. The function (t) is therefore assumed to be flat during the middle of the hour, but will ramp up or down at the start and the end of the hour. Hence, the segments of the curve rk (tk ) looks like the following:
rk (tk ) =
7. Modeling the net base production Q (t) In this section, the models of the net base production Q (t) = P(t) − W (t) ∈ R,
t≥0
(14)
are presented. The model of W (t) is described in Section 6.2. By the definition of the net base production in (14), no wind power production is included in Q (t). This implies that Q (t) is controllable. Similarly as for the base load and wind power production, the net base production is modeled as the sum of an algebraic function and a stochastic process: Q (t) = (t) + Y (t),
t ≥ 0,
(15)
where (t) ∈ R+ is an algebraic function and Y (t) ∈ R is a stochastic process. (t) represents the planned production if all units were 100% reliable while Y (t) simulates the unforeseen events leading to shut-downs of production units. The following can be noted regarding (t): • (t) is not a mean function of Q (t). This is because the deviation of Q (t) from (t) can only be negative since Y (t) represents unforeseen production decreases. • (t) is dependent on decisions made by the producers in the system. If assuming that the producers produce according to forecasts, the hourly integral of (t) will equal the difference between the last performed forecasts of the hourly load and the hourly wind power production. This is further discussed in Section 8.3. To model (t), piecewise linear curves rk (tk ) : [0, 1) → R+ are used: (t) = rk (tk ),
t ∈ [k + 1/2, k + 3/2).
(16)
The time variable tk ∈ [0, 1) take the value 0 at the middle of hour k and asymptotically approaches the value 1 at the middle of hour k + 1. This is illustrated in Fig. 4. The motivation for choosing a piecewise linear function for modeling (t) can be found in the market structure, where electricity is
⎧ ⎪ ⎨ bk , ⎪ ⎩
tk ∈ [0, e − 1/2),
hk (tk ), bk+1 ,
tk ∈ [ e − 1/2, s + 1/2), tk
∈ [ s
(17)
+ 1/2, 1),
where hk (tk ) = k tk + mk , e denotes the time at the end of the hour that the production starts to change, and s is the corresponding for the start of the hour. Since rk (tk ) should be a continuous function, rk (0) = bk and limtk →1 rk (tk ) = bk+1 . This also gives the values for k and mk as k = (bk+1 − bk )/( e − s ) and mk = bk (1 + ( s /( e − s ))) − bk+1 ( s /( e − s )). The principles of (t) are illustrated in Fig. 4. To consider unforeseen shut-down of units which causes momentarily decreases in production, Y (t) is modeled by using a stochastic jump process [12]. This process should also include a reversion property since the production should return to (t). This is because if a power station is taken off-line, this unit will be replaced by other base production. The owner of the power plant will replace the lost generation by increasing the generation in some other units or by purchasing power from other actors. Y (t) can be described by the following SDE: dY (t) = −ˇY (t)dt − dZ(t),
(18)
where ˇ ∈ R+ is the reversion factor and Z(t) ∈ R+ denotes a pure jump process. The physical interpretation of ˇ is how fast the lost base production will be covered with other base production in relation to the size of the production level of the unit taken off-line. Z(t) is defined by the stochastic sequence (Tn ∈ R+ , Vn ∈ R+ ) as Z(t) = Z(0) +
∞
I(Tn ≤t) Vn ,
t ≥ 0.
(19)
n=1
Thus, Z(t) is the sum of all events Vn that have occurred up to time t. Thereby, dZ(t) is a series of shocks occurring at times Tn of sizes Vn . The interpretation of this in the model of Y (t) is that dZ(t) represents events where base production units are taken off line due to failures or other unpredictable events. Now, the model of the net base production in (15) solves the following SDE: dQ (t) =
d(t) dt
+ ˇ((t) − Q (t)) dt − dZ(t).
(20)
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
Concerning the process Z(t), the differences Tn+1 − Tn are called arrival times and are assumed exponentially distributed. Hence Tn+1 − Tn ∈ exp( (Q (t)). The failure rate (Q (t)) is dependent on the number of units in operation at the current time and can therefore be expressed as a function of the production level Q (t). If the number of on-line units in the system is large at all production levels,
(Q (t)) can be approximated by a constant . Vn ≥ 0 are stochastic variables with a known distribution function FV (·; Q (t)) reflecting the base production in the on-line units in the system. Hence, the distribution function is dependent of Q (t). In reality, Vn is a discrete variable since the shut-down of units can only result in a finite number of states. If the number of on-line units in the system is large and are of different sizes, FV (·; Q (t)) can be approximated with a continuous function. Further (and similarly as for the failure rate), if the number of on-line units is large for all Q (t), FV (·; Q (t)) can be approximated by FV (·).
sponding historical data consists of the sum of the mean load levels: k =
k
i ∈ R+ .
• Hourly base load levels, ln . Depending on the modeled system, this can be actual load levels (assuming that the load do not participate in the frequency control), or load levels that are adjusted to eliminate the impact of frequency control. • Hourly wind power production levels, wn . These data can be historical data from a system containing a significant amount of wind power, but can also be a simulation result. Models for generating such data are, e.g. presented in [38,39]. Let the length of the available series of hourly loads and wind power productions be denoted by N. Typically, N = 8760 corresponding to data for one year. Let K be the time horizon of the performed simulations, e.g. K = 24 hours. The historical data can then be organized in NK = N/K sequences each of length K. Now let li,j , i = 1, . . . , NK , j = 1, . . . , K and wi,j , i = 1, . . . , NK , j = 1, . . . , K denote the load and wind power production in sequence i, hour j respectively. 8.1. Estimating L(t)
(21)
i=0
What now remains are the parameters ˛ and of the stochastic process X(t) in (7). These are estimated by applying the estimation method described in [34]. Also the estimation method described in [40] can be applied. The estimation is based on studying the energy t process S(t) = 0 (s) + X(s)ds since the available data reflects this process. The variance of S(t) is given by
t
t
=
Var[S(t)]
Cov[X(u), X(v)]dudv 0
0
t
t
= 0
8. Parameter estimation In this section, possible parameter estimation methods are presented. Note that the choice of estimation method depends on the available data. Hence, in situations where other data are available, alternative estimation methods may be preferable. In the ideal situation, the available data reflects the process that one wants to estimate. However, this is not always the case, and in some situations it might even be impossible to get such data. In this paper L(t), W (t) and Q (t) are continuous stochastic processes which parameters must be estimated. However, available data consists, in accordance with assumption (v) in Section 5, of hourly integrals of these processes. This implies that information regarding the variations within the hour is lost. The estimation procedures described below therefore includes “guesstimates” where there is lack of information. In cases where data actually are available, these should be used to make a better estimate of the model parameters. In this paper, the following data series are assumed available:
971
=
0
2 −˛|u−v| dudv e 2˛
(22)
2 1 t + (e−˛t − 1) . 2˛ ˛
By estimating the variance from historical data and putting these on the left hand side of the expression above, estimates of ˛ and can be retrieved. The corresponding historical data consists of accumulated energy levels, i,k : i,k =
k
li,j ∈ R+ .
(23)
j=0
Thus, by setting 2k = Vari [ i,k ] at the left hand side of (22), estimates of ˛ and can be extracted. 8.2. Estimating W (t) Since the wind power production is assumed log-normally distributed, the logarithm of the available data can be used to estimate the parameters of ω + U(t) according to the same methods as for L(t) presented in Section 8.1. Therefore let mi,j = ln wi,j and M(t) = ln W (t) = ω + U(t) ∈ R,
t ≥ 0.
(24)
The estimate of ω is the mean of mi,j , ω=
1
NK K
NK K
mi,j .
(25)
i=1 j=1
The parameters of the Ornstein–Uhlenbeck process U(t) are estimated by the same method as for the process X(t) of the base load. t The studied energy process is R(t) = 0 ω + U(s)ds since the available data are the logarithm of the hourly productions. The variance of this process is given by Var[R(t)] =
2 1 t + (e−˛t − 1) . 2
(26)
The corresponding historical data series consists of the sum
i,k
=
k
mi,j ∈ R.
(27)
j=0
As described in Section 6.1L(t) is the mean function of the base load, (t), modeled with a sum of cosine functions according to (10). The estimation of the parameters ak and dk is performed by applying the DCT. Since (t) is the mean function, mean values of the historical data are used in the estimation. Thus, estimation data are the hourly mean load levels k ,
k = 1, . . . , K, which NK can be calculated according to k = (1/NK ) i=1 li,k . Since the expression (10) represents the accumulated energy, the corre-
By setting the left hand side of (26) to 2 = Vari [ of and can be retrieved.
k
i,k ],
estimates
8.3. Estimating Q (t) The base production model defined by (15) includes the function (t), which is modeled by linear segments as described by (17).
972
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
By applying Itô calculus, this can be expressed as Cov[S(u), S(t)] 1 1 ˛(u∧t−u∨t) 2 − e−˛(u+t) ) . = 2 (u ∧ t) − (e−˛u + e−˛t )(e˛(u∧t) − 1) + (e ˛ 2˛ ˛
(34)
Now, the hourly load for hour k is Lk = S(k + 1) − S(k) and the covariance Cov[Lk , Lu ] thereby becomes Cov[Lk , Lu ] = Cov[S(k + 1) − S(k), S(u + 1) − S(u)] Fig. 5. Forecasting the hourly load levels when n = 2. The grey area represents the previous known load levels and the dashed lines possible outcomes of the future hourly load levels. The load forecast for hour k is desired, which is conditioned on the known load levels.
= Cov[S(k + 1), S(u + 1)]
(35)
−Cov[S(k + 1), S(u)] − Cov[S(k), S(u + 1)] +Cov[S(k), S(u)].
The parameters of (t) can be estimated by studying the hourly integrals of (t), which should reflect the difference between the forecasts of the hourly load and the hourly wind power production. For convenience introduce the following integrals: Lk =
k+1
k+1
k+1 k
L(t)dt,
Wk = k W (t)dt and k = k (t)dt. Now assume, in accordance with assumption (iv) in section 5, that the production plans are set n hours before the start of the actual hour. k can then be expressed as ˆ k |W1 , . . . , Wk−n , k = Lˆ k |L1 , . . . , Lk−n − W
Hence,
2
(
the
normal
(L , . . . , Lk−n , Lk )
vector
k−n+1
k+11
have
means
(t)dt, . . . , k−n (t)dt, k (t)dt) and covariances 1 given by (34) and (35). Thereby, Lˆ kn , which is defined by (30), can be calculated by using the normal conditional distribution fLk |L1 =l1 ,...,Lk−n =lk−n (x) =
(28)
fLk ,L1 ,...,Lk−n (x, l1 , . . . , lk−n ) fL1 ,...,Lk−n (l1 , . . . , lk−n )
(36)
,
ˆ k denotes the forecasted base load and wind power where Lˆ k and W ˆ n as the forecasted production, respectively. Introducing Lˆ kn and W k hourly base load and wind power production for hour k performed n + 1 hours before hour k (hence information regarding hour n is available, but not for hour n + 1), (28) becomes
whose parameters now are known. As for the load, the forecasted wind power production is modeled as the conditional expectation of the process,
ˆ n. k = Lˆ kn − W k
Since M(t) = ln W (t) have the same model structure as L(t), the same method as for E[Lk |L1 , . . . , Lk−n ] can be applied for the calˆn= culation of E[Mk |M1 , . . . , Mk−n ]. The forecast of Wk is then W k eE[Mk |M1 ,...,Mk−n ] . ˆ n have been calculated, the When the hourly values k = Lˆ kn − W k
(29)
The load forecasting problem is depicted in Fig. 5, which illustrates an example when n = 2. The shaded area represents the known information consisting of previous hourly loads up to hour k − 2. The calculation of k require some assumption regarding the forecasts of the load and wind power production. Starting with the load, the assumption is made that the forecast equals the conditional expectation Lˆ kn = E[Lk |L1 , . . . , Lk−n ].
(30)
This can be found be studying the energy process
+
T
S(T ) =
L(s)ds = 0
0
T
+ e−˛T
(31)
(e˛T − e˛s )dB(s) . 0
Observe that s = 0 is when the measurements started. Hence, T T is very large, implying that X(0)(1 − e−˛T ) e−˛T 0 (e˛T − e˛s )dB(s). Thereby, the part originating from the initial distribution of X(0) can be neglected. Hence,
T
S(T ) ≈
(s)ds
0
+
−˛T e ˛
(32)
T
(e˛T − e˛s )dB(s). 0
(S(u), S(t)) is a normal vector with expectations E[S(t)] = Thereby t (s)ds. The covariance can be calculated as follows prerequisite 0 that the approximation in (32) holds:
Cov[S(u), S(t)] =E
−˛u e ˛
u
(e˛u − e˛s )dB(s)
0
linear segments of the curve (t) are chosen so that k , i.e.
1 1 2
rk−1 (tk−1 )dtk−1 +
1 2
rk (tk )dtk = k ,
(37)
k+1 k
(t)dt =
(38)
0
which can be expressed as
T
(s)ds
1 X(0)(1 − e−˛T ) ˛
ˆ n = E[Wk |W1 , . . . , Wk−n ]. W k
×
−˛t e ˛
t
(e˛t − e˛s )dB(s)
.
0
(33)
Abk−1 + Bbk + Cbk+1 = k ,
(39)
where A, B and C are constants depending only on s and e . Thus, by solving the system of linear equations defined by (38), values of bk can be found and (t) is defined as the piecewise curve with segments rk (tk ). The production model defined by (15) also includes the stochastic process Y (t), which in turn includes the stochastic pairs (Tn , Vn ). The distributions of the arrival times Tn+1 − Tn ∈ exp( (Q (t)) and production decreases Vn can either be estimated from historical data of the real-time operation, or by studying the actual system and the generation units within it. Getting access to information concerning all the on-line generation units might not however be possible. One might then be forced to rely on historical data, which does not necessarily reflect the total situation in the system. For example, rare events might not be represented in historical data, or recent investments in units can have changed the distributions of Tn and Vn . Finally, the mean reversion parameter ˇ also need to be set. To estimate this parameter outgoing from historical data is certainly a challenge. Instead ˇ is set heuristically so that the reversion of the deviation from the function (t) is “reasonably” quick.
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
9. Case study To show the feasibility of the developed models, a case study consisting of a Monte Carlo simulation has been performed. The case study evaluates the impact of wind power on the demand for real-time balancing power in a fictive system. To get reasonable load curves, etc., available data from the Swedish power system has been used in the estimation process. These data consists of hourly load levels ln for the whole year of 2005 and hourly wind power production levels wn from [41] for one year. Hence, N = 8760. The simulations in this case study spans over one day and thereby K = 24. This gives NK = 365. The wind power data from [41] was scaled to represent 1000 MW of installed wind power. Concerning the distributions of Tn and Vn , some information regarding the available units in the system and times between failures was available from [17]. The assumed market structure allows for trading up to one hour before the start of the operational hour. This means that information regarding the hour starting two hours before the operational hour is assumed known, while information for the hour before the operational hour is not. Thereby, n = 2 was used for the forecasts of the hourly load and wind power production levels. The developed models and estimation methods were implemented by using the software package Matlab. 9.1. Parameter estimation The model parameters were estimated according to the methods described in Section 8. The parameter values of (t) and (t) are not presented due to limitations in space but are given on request. The parameters of the stochastic process X(t) of the base load were estimated to ˛ = 0.10, = 75.
(40)
The parameters of W (t) were estimated to ω = 5.3, = 0.01,
(41)
= 0.1. For the base production, the times s and e were set to s = 0.25 and e = 0.75. This implies that the production will start to ramp up or down 15 min before the hour starts and will end to increase/decrease 15 min after the start of the hour. Since the
973
market structure allows for trading up to one hour before the operational hour starts, the hourly levels k can be expressed as ˆ 2 . In this case study, Lˆ 2 and W ˆ 2 are for simplicity calk = Lˆ k2 − W k k k ˆ 2 = E[Wk |Wk−2 ]. culated as Lˆ 2 = E[Lk |Lk−2 ] and W k
k
The intensity of Tn+1 − Tn , i.e. (Q (t)) and the distribution of Vn are assumed independent of the current production level Q (t). Thus
(Q (t)) ≡ and FV (·, Q (t)) ≡ FV (·). was set to = 24 and hence, Tn+1 − Tn ∈ Exp(24). The density function fV (·) is assumed triangular with distribution function
fV (x) =
⎧ 1 ⎪ ⎨ fV , x ∈ [100, 500] ⎪ ⎩
fV2 ,
x ∈ [500, 1000]
0,
otherwise,
(42)
where fV1 (x) = −5.56 × 10−4 + 5.56 × 10−6 x,
(43)
fV2 (x) = 4.44 × 10−3 − 4.44 × 10−6 x.
(44)
Hence, the size of the shut-downs Vn will exist on the interval [100, 1000]. Further, the reversion parameters was set to ˇ = 24. 9.2. Simulations When the model parameters of L(t), W (t) and Q (t) have been set, sample path simulations of these processes can be performed and the associated sample paths of D(t) = Q (t) − N(t) = Q (t) − (L(t) − W (t)) can be calculated. In the simulations, two different cases were run: (1) With wind power, and (2) without wind power. The simulated sample paths of the load L(t) are the same for the both cases, while W (t) only were applied in case 1. Since Q (t) is dependent of forecasts of L(t) and W (t), different simulations of Q (t) were performed for the two cases. To distinguish between these two, the denotation Q #1 (t) refers to the net base production in case 1 and Q #2 (t) for case 2. The corresponding denotation for D(t) is D#1 (t) and D#2 (t). The presented models were used to generate 100 sample paths, each covering 24 h, of L(t), W (t), Q #1 (t) and Q #2 (t). The sample paths of L(t), W (t) and Q #1 (t) are displayed in Fig. 6. Note that since the algebraic function #1 (t) is dependent on successive forecasts of Lk and Wk , the #1 (t) s are different for each sample path. Concerning the differences between the cases with and without wind power, the associated hourly energy process of D(t), defined
k+1
as Dk = k D(t)dt, was computed for the various sample paths. The histograms for case 1 and case 2 of Dk are shown in Fig. 7. The
Fig. 6. Sample paths of L(t) (left), W (t) (middle) and Q #1 (t) (right). The black lines corresponds to the algebraic functions (t), eω(t) and #1 (t). Note that for Q #1 (t) there is one algebraic function per sample path.
974
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
Fig. 7. Histograms of dk#1 (black bars) and dk#2 (white bars).
sample mean and variance of Dk were also calculated: E[Dk#1 ] = −15.7, E[Dk#2 ] = −18.9, Var[Dk#1 ] = 6558,
(45)
Var[Dk#2 ] = 6332. Hence, with wind power, the mean moves slightly towards zero while the variance increases. 10. Conclusion This paper presents models based on SDEs of base production, base consumption and wind power production. These models can be used to evaluate the impact of wind power on the real-time balancing of the power system. This is of great interest for system operators as well as for producers and consumers. Using continuous time models enables the possibility of considering counter balancing. Besides the description of the models, estimation methods for parameters are proposed and a case study is presented. The case study suggests that using SDEs is a possible approach for modeling imbalances in continuous time. Future research concerning this topic could consist of performing more detailed case studies where the developed models are applied. Acknowledgement Support from the Nordic power exchange Nord Pool considering market data is gratefully acknowledged. References [1] R. Weron, Modeling and Forecasting Electricity Loads and Prices—A Statistical Approach, John Wiley & Sons Inc., 2006. [2] N. Amjady, Short-term hourly load forecasting using time-series modeling with peak load estimation capability, IEEE Trans. Power Syst. 16 (3). [3] E. Lindgren, L. Söder, Minimizing regulating costs in multiarea systems with uncertain wind power forecasts, Wind Energy 11 (1) (2008) 97–108.
[4] H. Banakar, C. Luo, B.T. Ooi, Impacts of wind power minute-to-minute variations on power system operation, IEEE Trans. Power Syst. 23 (1) (2008). [5] K. Elkington, V. Knazkins, M. Ghandhari, On the stability of power systems containing doubly fed induction generator-based generation, Electric Power Syst. Res. 78 (9). [6] I. Norheim, K. Uhlen, J.O. Tande, T. Toftevang, M.P. Pálsson, Doubly fed induction generator model for power system simulation tools, in: Proceeding of the NordicWind Power Conference, Gothenburg, Sweden, 2004. [7] T. Petru, T. Thiringer, Modeling of wind turbines for power system studies, IEEE Trans. Power Syst. 17 (4). [8] V. Akhmatov, H. Knudsen, An aggregated model of a gridconnected, large-scale, offshore wind farm for power system stability investigation—importance of windmill mechanical system, Int. J. Electrical Power Energy Syst. 24 (9) (2002). [9] H. Holttinen, Impact of hourly wind power variations on the system operation in the Nordic countries, Wind Energy 8 (2). [10] Wind power plants and system operation in the hourly time domain, NREL, 2003. URL: http://www.nrel.gov. [11] B. Oksendal, Stochastic Differential Equations, Springer-Verlag, Berlin Heidelberg, 2005. [12] F.C. Klebaner, Introduction to Stochastic Calculus with Applications, Imperial Collage Press, 2005. [13] European Energy Exchange. URL: http://www.eex.com. [14] Scheduling operations, PJM Manual 11: PJM, 2007. URL: http://www.pjm.com. [15] Nordic grid code 2004. Nordel, 2004. URL: http://www.nordel.org. [16] The Grid Code, National Grid Electricity Transmission Plc., 2008. URL: http://www.nationalgrid.com/. [17] Nord Pool AS web site. URL: http://www.nordpool.com. [18] Operador del Mercado Ib’erico de Energia web site. URL: http://www.omel.es. [19] APX Group. URL: http://www.apxgroup.com. [20] Market Operations, ISO New England Inc., 2006. URL: http://www.iso-ne.com. [21] Operating Procedure: Frequency Control Ancilliry Services, National Electricity Market Management Company (NEMMCO), 2005. URL: http://www.nemmco.com.au. [22] Current State of Balance Management in Europe, ETSO, 2003. URL: http://www.etso-net.org. [23] D. Ranaweera, G.G. Karady, R.G. Farmer, Economic impact analysis of load forecasting, IEEE Trans. Power Syst. 12 (3) (1997) 1388. [24] M.A. Abu-El-Magd, N.K. Sinha, Short-term load demand modeling and forecasting: a review, IEEE Trans. Syst., Man, Cybern. SMC-12 (3) (1982) 370. [25] I. Moghram, S. Rahman, Analysis and evaluation of five shortterm load forecasting techniques, IEEE Trans. Power Syst. 4 (4). [26] S. Fan, L. Chen, Short-term load forecasting based on an adaptive hybrid method, IEEE Trans. Power Syst. 21 (1) (2006) 392–401. [27] N. Amjady, Short-term bus load forecasting of power systems by a new hybrid method, IEEE Trans. Power Syst. 22 (1) (2007) 333–341. [28] Z. Yun, Z. Quan, S. Caixin, L. Shaolan, L. Yuming, S. Yang, RBF neural network and ANFIS-based short-term load forecasting approach in real-time price environment, IEEE Trans. Power Syst. 23 (3) (2008) 853–858. [29] M. Milligan, A.H. Miller, F. Chapman, Estimating the economic value of wind forecasting to utilities, URL: http://www.nrel.gov/publications/, in: Proceedings of Windpower’ 95, Washington DC, USA, 1995. [30] S.J. Watson, G. Giebel, A. Joensen, The economic value of accurate wind power forecasting to utilities, in: Proceedings of European Wind Energy Conference, Nice, France, 1999. [31] T.S. Nielsen, Prediction of regional wind power, in: Proceedings of Global Wind Power Conference, Paris, France, 2002. [32] U. Focken, M. Lange, H.M. Waldl, Previento—a wind power prediction system with an innovative upscaling algorithm, in: Proceedings of European Wind Power Conference, Copenhagen, Denmark, 2001. [33] G. Giebel, R. Brownsword, G. Kariniotakis, The state-of-theart in short-term prediction of wind power. A literature review, Tech. Rep., EU Project ANEMOS. URL: http://anemos.cma.fr. [34] M. Perninge, Evaluating the uncertainties involved in net transmission capacity calculation, Licentiate Thesis, School of Electrical Engineering, KTH, Stockholm, Sweden, 2009. [35] F. Dornier, M. Queruel, Caution to the wind, Weather risk special report, Energy & power risk management, Risk Magazine, 2000. [36] N. Ahmed, T. Natarajan, K.R. Rao, Discrete cosine transform, IEEE Trans. Comput. C-23. [37] K. Conradsen, L.B. Nielsen, L.P. Prahm, Review of Weibull statistics for estimation of wind speed distributions, J. Clim. Appl. Meteorol. 23. [38] J. Matevosyan, L. Söder, Minimization of imbalance cost trading wind power on the short-term power market, IEEE Trans. Power Syst. 21(3). [39] G. Papaefthymiou, B. Klöckl, MCMC for wind power simulation, IEEE Trans. Energy Convers. 23 (1). [40] S. Ditlevsen, M. Sorensen, Inference for observations of integrated diffusion processes, Scand. J. Stat. 31. [41] U. Arvidsson, R. Murray, V. Neimane, 4000 MW wind power in Sweden: impact on regulation and reserve requirements, Tech. Rep. 05:19, Elforsk AB, 2005.
Contents lists available at ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
Modeling real-time balancing power demands in wind power systems using stochastic differential equations Magnus Olsson ∗ , Magnus Perninge, Lennart Söder Royal Institute of Technology (KTH), Electric Power Systems, Teknikringen 33, 10044 Stockholm, Sweden
a r t i c l e
i n f o
Article history: Received 7 August 2008 Received in revised form 18 February 2009 Accepted 18 February 2009 Available online 2 February 2010 Keywords: Stochastic differential equations Real-time balancing Wind power
a b s t r a c t The inclusion of wind power into power systems has a significant impact on the demand for real-time balancing power due to the stochastic nature of wind power production. The overall aim of this paper is to present probabilistic models of the impact of large-scale integration of wind power on the continuous demand in MW for real-time balancing power. This is important not only for system operators, but also for producers and consumers since they in most systems through various market solutions provide balancing power. Since there can occur situations where the wind power variations cancel out other types of deviations in the system, models on an hourly basis are not sufficient. Therefore the developed model is in continuous time and is based on stochastic differential equations (SDE). The model can be used within an analytical framework or in Monte Carlo simulations. © 2010 Elsevier B.V. All rights reserved.
1. Algebraic symbols Variables in continuous time are denoted with the time argument with parenthesis, as e.g. P(t), while discrete time is indicated with subscripts, as e.g. Pk . Stochastic variables are in uppercase Latin letters, e.g. Pk . Realizations or historical data of the same are in lowercase letters, as e.g. pk . The indicator function I¯(A) ∈ {0, 1} take the value I¯(A) = 1 if the logical expression A is true, otherwise I¯(A) = 0. Below follows a list of the variables and most important parameters. Additional parameters are defined on first occasion. D(t) P(t) L(t) Q (t) W (t) N(t) (t) ak , dk X(t) ˛, ω(t) U(t) ,
demand for balancing power [MW] base production [MW] base load [MW] controllable production [MW] uncontrollable production [MW] net base load [MW] mean function of L(t) parameters of (t) stochastic process in L(t) parameters of X(t) deterministic part of W (t) stochastic process in W (t) parameters of U(t)
∗ Corresponding author. E-mail addresses: [email protected] (M. Olsson), [email protected] (M. Perninge), [email protected] (L. Söder). 0378-7796/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.epsr.2010.01.004
(t) Y (t) rk (tk ) ˇ Z(t) (Tn , Vn )
deterministic part of Q (t) stochastic process in Q (t) segments of (t) parameter in Y (t) stochastic jump process in Y (t) stochastic pairs defining Z(t)
2. Introduction This paper addresses the problem of balance management in power systems. The aim is to describe a new probabilistic model of the time continuous balancing demand in MW in a system with large amounts of wind power. This model can, e.g. be used to estimate energy demands for primary, secondary and tertiary control. In addition to be useful within an analytical framework, the presented model can be applied in Monte Carlo simulations. Performing such analyses and simulations is of great importance not only for system operators, but also for producers (and possibly consumers) since they in most systems provide balancing power through various market solutions. The presented model provides a tool for calculation and simulation of the demand for frequency control and required ramp rates. Probabilistic load or production models are usually based on an hourly time resolution as in, e.g. [1] and [2]. However, to be able to perform precise estimations of the impact of new power sources into the power system on the continuous need for balancing power, models using an hourly time resolution are not sufficient. This is due to the fact that the production may vary significantly during the hour, as e.g. in the case of wind power.
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
The inclusion of wind power into the power system has a significant impact on the need for real-time balancing power due to the stochastic nature of the wind power production. However, the introduced unpredictable wind power production variations will also sometimes help the system in the real-time balancing process. This is in this article denoted as counter balancing. In time discrete models with, e.g. an hourly time resolution, assumptions regarding what happens during the hour become necessary to consider counter balancing. Since this can occur at any time and span a time period of arbitrary length, models in continuous time are more suitable to capture this. Concerning intra-hourly modeling, not much work has been published. In [3], an optimization model for minimization of the real-time balancing cost in systems with uncertain wind power forecasts is presented. The proposed model suggests which bids to call and when this should be performed. The model in [3] and the model presented in this paper address the same problem, but from different perspectives. The model in [3] focuses on optimizing the bids given the load and the scheduled production, while the model presented in this paper represents the continuous demand for balancing power. In [4], a simulation model of the impact of wind power production on operation of the power system is presented. However, the presented model does not have a stochastic representation. In papers [5–8] detailed models of wind turbines or farms are presented for analyses of the interaction between the power system and the wind farms. Such models are in theory possible to use for estimation of the increased demand for real-time balancing power associated with integration of wind power, but simulating a system with such accuracy for a 24 h period, for example, is not possible due to computational time. Generally, these types of simulation models concern time frames of a few seconds in length. Articles concerning the impact of wind power on real-time balancing quantities include [9] and [10]. However, these have a time resolution of one hour and thus do not include counter balancing within the hour. This paper presents stochastic models in continuous time of loads and productions in order to capture the continuous demand for balancing power. Suitable tools for modeling stochastic processes in continuous time are stochastic differential equations (SDEs) [11]. SDEs have previously been used in various applications including biology, physics, signal processing and finance [12]. 3. Background 3.1. Market structure Some of the properties that need to be considered in the modeling process are set by the market rules. These are: • The stipulated trading period length. This will affect the dispatch of power into the system since the producers tend to change the production at the borders of the trading periods where the price changes. • The possibility for actors to trade power according to, e.g. updated forecasts before the start of the actual trading period. Considering trading periods, trading is usually performed on an hourly basis as on the European Energy Exchange (EEX) [13], PJM [14], and in the Nordic system [15]. In some systems, as e.g. UK [16], half-hour periods are instead used. For convenience, this paper assumes a trading period of one hour, but the suggested model can be applied to any period length. In this paper, no specific market structure is assumed. However, one market property needs to be considered: At which time the actors no longer have the possibility to trade power ahead for a
967
specific hour. This needs to be considered in the modeling due to forecast uncertainties. In many power systems, ahead trading can be performed on the following market places: • Day-ahead spot market, where actors can trade electric energy for the various hours of the next day. Examples of day-ahead spot markets are EEX [13], Nord Pool Elspot [17] and Omel [18]. • Intra-day market, which offers possibilities to adjust the traded energy quantities from the day-ahead market according to updated forecasts or to non-accepted bids on the day-ahead market. Here it is possible to continuously trade energy for a specific hour after the day-ahead market has closed, until a few hours before the start of that hour. Examples of intra-day markets are Nord Pool Elbas [17] and APX NL [19]. Hence, in this setting, the interesting market property is when the intra-day market closes for the various hours since this represents the last possibility to trade energy ahead. Apart from the ahead markets, there can also exist a realtime balancing market where producers and consumers have the possibility to trade power with the system operator during the operational hour. This is a part of the frequency control often referred to as tertiary control. Examples of such markets are the real-time balancing market in PJM [14], regulating market in the Nordic system [15], the real-time energy market in New England [20], and frequency control ancillary services (FCAS) market in Australia [21]. A compilation of balance management in Europe can be found in [22]. The model described in this paper does not require any assumption regarding the real-time market. However, the model can, e.g. be applied in estimations or simulations of demands for real-time balancing power. 3.2. Load forecasting Load forecasting encompasses forecasts on several time scales, ranging from hourly to yearly. The models presented in this paper consider a number of hours and include load forecasting on an hourly time scale. This means that load forecasts in MWh/h are made each hour according to updated information for individual future hours. Load forecasting has been subjected to substantial amounts of research. The motivation for development of accurate forecast methods lies in the economic impact of the forecast errors [23]. Short-term load modeling and forecasting is reviewed in [24,25]. Recent work on short-term forecasting includes [26–28]. 3.3. Wind power forecasting The uncertainties associated with wind power originates from uncertainties in wind speed forecasts. This in combination with the rapid development of wind farms generates a need for better wind power production forecasting methods. The higher forecast reliability the lower the balancing costs become in the system, which in case of large-scale integration of wind power can imply substantial savings for the wind farm owners as well as better overall efficiency of the system [29,30]. Different wind power forecast techniques are presented in e.g. [31,32] and the state-of-the-art literature review [33]. 4. Definitions 4.1. Continuous demand for balancing power The overall aim of this paper is to model the continuous demand for balancing power D(t), measured in MW. To be able to achieve
968
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
4.2. Controllable and uncontrollable generation
Fig. 1. Basic system with mechanical production P(t), electrical production P e (t) and load L(t).
this, no dispatch for frequency control is included in the developed models. Hence, dispatch related to primary, secondary or tertiary control is excluded in the model of D(t). The continuous demand for balancing power is now defined as the difference between the mechanical production P(t) and the electrical load L(t): D(t) = P(t) − L(t).
(1)
Mechanical production here refers to the power from the turbines, originating from, e.g. water or steam hitting the blades. The difference between mechanical and electrical production is pointed out in the fundamental system in Fig. 1. Note that if replacing P(t) with the electrical production P e (t) in (1) will give D(t) ≡ 0 since P e (t) ≡ L(t) by the law of conservation of energy! The difference between P(t) and P e (t) is covered by changes in rotational speed of the synchronous machines in the system. The actual mechanical production, including wind power production, unforeseen shut-downs of generation units, etc. but with excluded frequency control is denoted base production throughout this paper. The corresponding electric load is denoted base load. To clarify this terminology, various schematic production and load curves are illustrated in Fig. 2. Starting with line PL in the figure, this line represents the electrical production and load including frequency control. Hence, this is the actual electrical production/consumption in the system. The other lines in the figure represent the following: The solid lines represent the planned mechanical production (black line P1) and forecasted load (grey line L1); the dashed black line P2 corresponds to the mechanical base production including e.g. unforeseen shut-downs of units and variations in mechanical wind power production; and the dashed grey line L2 is the base load. Note that in P2 and L2, neither production nor consumption changes for frequency control are included. Concerning the load, the lines L2 and PL do not coincide, implying that some of the frequency control has been performed on the load side in this figure. In terms of the lines in Fig. 2, the overall aim of this paper is to model the difference between P2 and L2, i.e. base production minus base load. The various production types in this paper consider mechanical productions, but for simplicity, the term “mechanical” is from now on omitted.
The wind power is in this paper considered as uncontrollable, implying that there are no possibilities to control the power output. This assumption can also be applied on the base load. Contrarily, except for the wind power production, the base production is controllable since it is the sum of decisions made by the producers in the power system. Thereby two different types of processes are considered: • Uncontrollable: Base load and wind power production, which are processes independent on forecasts and decisions. • Controllable: Base production with excluded wind power production. This is a result of the decisions of the producers, which are based on forecasts. Note that controllable does not imply that the production is deterministic. For example faults can occur resulting in shut-downs of production units. The uncontrollable base load and wind power production are combined into the net base load N(t), N(t) = L(t) − W (t),
(2)
where L(t) denotes the base load and W (t) the wind power production. Hence, N(t) represents the base load that is not covered by the wind power production. The corresponding for base production can be expressed as Q (t) = P(t) − W (t),
(3)
where Q (t) is the net base production and P(t) is the total base production including wind power. Q (t) thereby represents the hydrothermal base production in the system. The definitions of N(t) in (2) and Q (t) in (3) imply that N(t) is uncontrollable while Q (t) is controllable. Q (t) is dependent on decisions made by the market actors, which in turn are based on the market structure and the available forecasts of the net load. Because of the differences between (L(t), W (t)) and Q (t), the modeling is divided into two main parts: (1) Modeling L(t) and W (t): L(t) and W (t) defines N(t) according to (2). In this paper, L(t) and W (t) are considered as endogenous variables. (2) Modeling Q (t): The aim of the model of Q (t) is to simulate the sum of the actions taken by the producers, resulting in the total base production excluding wind power in the system. Q (t) is dependent on forecasts of the base load and the wind power production. The aim to model the demand of balancing power can be reached by studying the difference between P(t) and L(t) (i.e. the difference between P2 and L2 in Fig. 2). This reflects the continuous demand for balancing power, D(t), and hence D(t) = P(t) − L(t) = Q (t) − N(t).
(4)
The advantage of studying Q (t) and N(t) instead of P(t) and L(t) is the distinction between controllable and uncontrollable processes. 5. Main assumptions The proposed models in this paper assume the following:
Fig. 2. Load and production curves. The black lines correspond to productions, where P1 is the planned production and P2 is the base production P(t). The grey lines correspond to base loads where line L1 represents the forecasted base load and L2 the actual base load L(t). Line PL represents the electrical production and consumption when frequency control is included.
(i) The wind power production is independent of the base load. (ii) The wind power production does not change discontinuously. This assumption can be interpreted as that the wind power plants are 100% reliable and will never be subjected to shutdowns. Another interpretation is that each plant is of arbitrarily
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
small size. This together with the assumption that the transmission lines are 100% reliable, and that the probability of shut-downs of two units at the exact same time is zero, implies that no shut-downs will cause any discontinuities in the total wind power production. (iii) The producers (excluding wind power production) deliver a constant base production output during the hour unless the system operator orders a production change. (iv) The market structure allows for trading of electric energy in MWh/h up to a few hours before the start of the operational hour. Market actors trade power according to updated forecasts of base load and wind power production until this market closes. Hence, when the intra-day market closes the following is valid concerning the hourly quantities: • Wind power producers are in balance considering traded quantity and forecasted production. • Other producers are in balance considering traded quantity and planned production. • Consumers are in balance considering traded quantity and forecasted load. (v) Available data series are hourly base load levels and wind power production levels. Both are measured in MWh/h. 6. Modeling the net base load N(t) In this section, the models of the base load N(t) = L(t) − W (t) ∈ R,
t≥0
(5)
are presented. The ingoing parts L(t) ∈ R+ and W (t) ∈ R+ are, as stated in main assumption (i) in section 5, assumed independent and are presented separately. 6.1. Modeling L(t)
Fig. 3. Example of mean load curve. The staircase-shaped curves are the hourly load levels and the solid line represents the model of (t) given by the DCT.
Now, the process L(t) = (t) + X(t) solves the following SDE: dL(t) =
L(t) = (t) + X(t),
t ≥ 0.
(6)
d(t) dt
+ ˛ ( (t) − L(t) ) dt + dB(t).
(9)
The inclusion of d(t)/dt is necessary in order to make L(t) a meanreverting process [35]. What is now left is the mean function (t), which can be chosen in several ways. Here a model based on the discrete cosine transform (DCT) [36] is used. The available data do, as further described in Section 8, in general not consist of samples of the total base load in MW, but of energy levels during the trading periods in MWh. Thus, the time integral of the base load is modeled, which then can be differentiated to get the load itself. The integral of (t) ≥ 0 is non-decreasing in time and therefore, also a linear term is added to the model. Hence, (t) can be written
The model of the load is briefly described below. A more detailed description can be found in [34]. The base load is modeled as the sum of a function representing the mean of the load and a stochastic process:
969
t
(s)ds = 0
a0 t +
N
ak dk cos
(2k − 1)(t − 1) 2N
(10) ,
k=1
The mean (t) ∈ R+ should resemble the continuous variation of the base load between various trading hours in order to get the typical load characteristics. The stochastic process X(t) ∈ R is chosen according to the following:
where N is the number of samples, a0 is the linear coefficient, and ak and dk are given by the DCT. An example of (t) can be found in Fig. 3 where the staircase-shaped curve corresponds to historical hourly load mean values and the solid line represents the model.
• It seems reasonable to assume, in accordance with the central limit theorem, that L(t) is normally distributed since it corresponds to a sum of various loads in the system. Hence, X(t) should be normally distributed. • The square mean of L(t) cannot, e.g. increase during many successive hours, but should instead return to the mean. This implies that some sort of mean-reversion is necessary.
6.2. Modeling W (t)
Putting the above together, the deviation of the load from its mean function is modeled with an SDE of Ornstein–Uhlenbeck type [12] according to: dX(t) = −˛X(t)dt + dB(t),
˛ ∈ R+ ,
(7)
where (B(t), t ≥ 0) is a Brownian motion. The solution to this SDE can be found by applying Itô’s formula [11] and is given by
X(t) = e−˛t (X(0) +
t
e˛s dB(s)),
(8)
0
where X(0) is the initial condition for the base load deviation from the mean and the integral is taken in the Itô sense.
The following can be noted regarding wind power production: • Wind speeds do not have symmetrical distributions [37]. Thereby, the distribution of W (t) should not be symmetrical. • The square mean of W (t) cannot, e.g. increase during many successive hours. This means that some mean-reversion should be included in the model. • The wind power production W (t) do not show the same strong cyclic behavior as the load L(t) on a daily basis. • The wind power production is bounded on the interval [0, W ], ¯ is the total installed wind power production capacity. where W The wind speed is often considered to be Weibull distributed. However, in this paper the log-normal distribution is applied for W (t) due to the benefit of having normally distributed logarithms. This is an advantage from a mathematical point of view because of the possibility of using well established SDEs based on Brownian motion.
970
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
Based on the properties above and on assumption (ii) in Section 5, the wind power production is modeled as W (t) = eω(t)+U(t) ∈ R+ ,
t ≥ 0,
(11)
where ω(t) ∈ R+ is an algebraic function representing the mean of ln W (t), and U(t) ∈ R is an Ornstein–Uhlenbeck process dU(t) = −U(t)dt + dB(t),
t≥0
(12)
where B(t) is a Brownian motion. Thereby, U(t) is normally distributed and is a diffusion, which can be expressed as
t
U(t) = e−t (U(0) +
es dB(s)).
(13)
0
ω(t) should reflect cyclic patterns of the wind power production. The modeling of ω(t) depends on the wind situation in the actual power system. In this paper, ω(t) is assumed constant and will further on be denoted without time argument, i.e. ω. ¯ , The expression in (11) can result in situations where W (t) > W which is not possible in reality. However, if the wind farms are widely spatially distributed across the system, the probability that ¯ becomes very small. Therefore, this is neglected in the W (t) > W developed model. If applying the model of W (t) within an analytical framework, it might become necessary to handle the case ¯ . In Monte Carlo simulation applications, the sample W (t) > W paths can be truncated to not violate the bounds of W (t). However, ¯ is very small, this is generally not if the probability that W (t) > W a problem in such simulations.
Fig. 4. Principle behavior of (t). The time variable t ∈ [0, ∞), while tk ∈ [0, 1) with tk = 0 when t = k.
traded by the hour. This together with assumption (iii) in Section 5, results in that changes in base production tends to be performed in steps at the border of each hour. The system operators often have the possibility to order the producers to start production earlier, e.g. up to 15 min in front of and/or at the end of hour. This is to smooth out the steps in base production that otherwise will occur in the system at the start of each hour. The function (t) is therefore assumed to be flat during the middle of the hour, but will ramp up or down at the start and the end of the hour. Hence, the segments of the curve rk (tk ) looks like the following:
rk (tk ) =
7. Modeling the net base production Q (t) In this section, the models of the net base production Q (t) = P(t) − W (t) ∈ R,
t≥0
(14)
are presented. The model of W (t) is described in Section 6.2. By the definition of the net base production in (14), no wind power production is included in Q (t). This implies that Q (t) is controllable. Similarly as for the base load and wind power production, the net base production is modeled as the sum of an algebraic function and a stochastic process: Q (t) = (t) + Y (t),
t ≥ 0,
(15)
where (t) ∈ R+ is an algebraic function and Y (t) ∈ R is a stochastic process. (t) represents the planned production if all units were 100% reliable while Y (t) simulates the unforeseen events leading to shut-downs of production units. The following can be noted regarding (t): • (t) is not a mean function of Q (t). This is because the deviation of Q (t) from (t) can only be negative since Y (t) represents unforeseen production decreases. • (t) is dependent on decisions made by the producers in the system. If assuming that the producers produce according to forecasts, the hourly integral of (t) will equal the difference between the last performed forecasts of the hourly load and the hourly wind power production. This is further discussed in Section 8.3. To model (t), piecewise linear curves rk (tk ) : [0, 1) → R+ are used: (t) = rk (tk ),
t ∈ [k + 1/2, k + 3/2).
(16)
The time variable tk ∈ [0, 1) take the value 0 at the middle of hour k and asymptotically approaches the value 1 at the middle of hour k + 1. This is illustrated in Fig. 4. The motivation for choosing a piecewise linear function for modeling (t) can be found in the market structure, where electricity is
⎧ ⎪ ⎨ bk , ⎪ ⎩
tk ∈ [0, e − 1/2),
hk (tk ), bk+1 ,
tk ∈ [ e − 1/2, s + 1/2), tk
∈ [ s
(17)
+ 1/2, 1),
where hk (tk ) = k tk + mk , e denotes the time at the end of the hour that the production starts to change, and s is the corresponding for the start of the hour. Since rk (tk ) should be a continuous function, rk (0) = bk and limtk →1 rk (tk ) = bk+1 . This also gives the values for k and mk as k = (bk+1 − bk )/( e − s ) and mk = bk (1 + ( s /( e − s ))) − bk+1 ( s /( e − s )). The principles of (t) are illustrated in Fig. 4. To consider unforeseen shut-down of units which causes momentarily decreases in production, Y (t) is modeled by using a stochastic jump process [12]. This process should also include a reversion property since the production should return to (t). This is because if a power station is taken off-line, this unit will be replaced by other base production. The owner of the power plant will replace the lost generation by increasing the generation in some other units or by purchasing power from other actors. Y (t) can be described by the following SDE: dY (t) = −ˇY (t)dt − dZ(t),
(18)
where ˇ ∈ R+ is the reversion factor and Z(t) ∈ R+ denotes a pure jump process. The physical interpretation of ˇ is how fast the lost base production will be covered with other base production in relation to the size of the production level of the unit taken off-line. Z(t) is defined by the stochastic sequence (Tn ∈ R+ , Vn ∈ R+ ) as Z(t) = Z(0) +
∞
I(Tn ≤t) Vn ,
t ≥ 0.
(19)
n=1
Thus, Z(t) is the sum of all events Vn that have occurred up to time t. Thereby, dZ(t) is a series of shocks occurring at times Tn of sizes Vn . The interpretation of this in the model of Y (t) is that dZ(t) represents events where base production units are taken off line due to failures or other unpredictable events. Now, the model of the net base production in (15) solves the following SDE: dQ (t) =
d(t) dt
+ ˇ((t) − Q (t)) dt − dZ(t).
(20)
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
Concerning the process Z(t), the differences Tn+1 − Tn are called arrival times and are assumed exponentially distributed. Hence Tn+1 − Tn ∈ exp( (Q (t)). The failure rate (Q (t)) is dependent on the number of units in operation at the current time and can therefore be expressed as a function of the production level Q (t). If the number of on-line units in the system is large at all production levels,
(Q (t)) can be approximated by a constant . Vn ≥ 0 are stochastic variables with a known distribution function FV (·; Q (t)) reflecting the base production in the on-line units in the system. Hence, the distribution function is dependent of Q (t). In reality, Vn is a discrete variable since the shut-down of units can only result in a finite number of states. If the number of on-line units in the system is large and are of different sizes, FV (·; Q (t)) can be approximated with a continuous function. Further (and similarly as for the failure rate), if the number of on-line units is large for all Q (t), FV (·; Q (t)) can be approximated by FV (·).
sponding historical data consists of the sum of the mean load levels: k =
k
i ∈ R+ .
• Hourly base load levels, ln . Depending on the modeled system, this can be actual load levels (assuming that the load do not participate in the frequency control), or load levels that are adjusted to eliminate the impact of frequency control. • Hourly wind power production levels, wn . These data can be historical data from a system containing a significant amount of wind power, but can also be a simulation result. Models for generating such data are, e.g. presented in [38,39]. Let the length of the available series of hourly loads and wind power productions be denoted by N. Typically, N = 8760 corresponding to data for one year. Let K be the time horizon of the performed simulations, e.g. K = 24 hours. The historical data can then be organized in NK = N/K sequences each of length K. Now let li,j , i = 1, . . . , NK , j = 1, . . . , K and wi,j , i = 1, . . . , NK , j = 1, . . . , K denote the load and wind power production in sequence i, hour j respectively. 8.1. Estimating L(t)
(21)
i=0
What now remains are the parameters ˛ and of the stochastic process X(t) in (7). These are estimated by applying the estimation method described in [34]. Also the estimation method described in [40] can be applied. The estimation is based on studying the energy t process S(t) = 0 (s) + X(s)ds since the available data reflects this process. The variance of S(t) is given by
t
t
=
Var[S(t)]
Cov[X(u), X(v)]dudv 0
0
t
t
= 0
8. Parameter estimation In this section, possible parameter estimation methods are presented. Note that the choice of estimation method depends on the available data. Hence, in situations where other data are available, alternative estimation methods may be preferable. In the ideal situation, the available data reflects the process that one wants to estimate. However, this is not always the case, and in some situations it might even be impossible to get such data. In this paper L(t), W (t) and Q (t) are continuous stochastic processes which parameters must be estimated. However, available data consists, in accordance with assumption (v) in Section 5, of hourly integrals of these processes. This implies that information regarding the variations within the hour is lost. The estimation procedures described below therefore includes “guesstimates” where there is lack of information. In cases where data actually are available, these should be used to make a better estimate of the model parameters. In this paper, the following data series are assumed available:
971
=
0
2 −˛|u−v| dudv e 2˛
(22)
2 1 t + (e−˛t − 1) . 2˛ ˛
By estimating the variance from historical data and putting these on the left hand side of the expression above, estimates of ˛ and can be retrieved. The corresponding historical data consists of accumulated energy levels, i,k : i,k =
k
li,j ∈ R+ .
(23)
j=0
Thus, by setting 2k = Vari [ i,k ] at the left hand side of (22), estimates of ˛ and can be extracted. 8.2. Estimating W (t) Since the wind power production is assumed log-normally distributed, the logarithm of the available data can be used to estimate the parameters of ω + U(t) according to the same methods as for L(t) presented in Section 8.1. Therefore let mi,j = ln wi,j and M(t) = ln W (t) = ω + U(t) ∈ R,
t ≥ 0.
(24)
The estimate of ω is the mean of mi,j , ω=
1
NK K
NK K
mi,j .
(25)
i=1 j=1
The parameters of the Ornstein–Uhlenbeck process U(t) are estimated by the same method as for the process X(t) of the base load. t The studied energy process is R(t) = 0 ω + U(s)ds since the available data are the logarithm of the hourly productions. The variance of this process is given by Var[R(t)] =
2 1 t + (e−˛t − 1) . 2
(26)
The corresponding historical data series consists of the sum
i,k
=
k
mi,j ∈ R.
(27)
j=0
As described in Section 6.1L(t) is the mean function of the base load, (t), modeled with a sum of cosine functions according to (10). The estimation of the parameters ak and dk is performed by applying the DCT. Since (t) is the mean function, mean values of the historical data are used in the estimation. Thus, estimation data are the hourly mean load levels k ,
k = 1, . . . , K, which NK can be calculated according to k = (1/NK ) i=1 li,k . Since the expression (10) represents the accumulated energy, the corre-
By setting the left hand side of (26) to 2 = Vari [ of and can be retrieved.
k
i,k ],
estimates
8.3. Estimating Q (t) The base production model defined by (15) includes the function (t), which is modeled by linear segments as described by (17).
972
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
By applying Itô calculus, this can be expressed as Cov[S(u), S(t)] 1 1 ˛(u∧t−u∨t) 2 − e−˛(u+t) ) . = 2 (u ∧ t) − (e−˛u + e−˛t )(e˛(u∧t) − 1) + (e ˛ 2˛ ˛
(34)
Now, the hourly load for hour k is Lk = S(k + 1) − S(k) and the covariance Cov[Lk , Lu ] thereby becomes Cov[Lk , Lu ] = Cov[S(k + 1) − S(k), S(u + 1) − S(u)] Fig. 5. Forecasting the hourly load levels when n = 2. The grey area represents the previous known load levels and the dashed lines possible outcomes of the future hourly load levels. The load forecast for hour k is desired, which is conditioned on the known load levels.
= Cov[S(k + 1), S(u + 1)]
(35)
−Cov[S(k + 1), S(u)] − Cov[S(k), S(u + 1)] +Cov[S(k), S(u)].
The parameters of (t) can be estimated by studying the hourly integrals of (t), which should reflect the difference between the forecasts of the hourly load and the hourly wind power production. For convenience introduce the following integrals: Lk =
k+1
k+1
k+1 k
L(t)dt,
Wk = k W (t)dt and k = k (t)dt. Now assume, in accordance with assumption (iv) in section 5, that the production plans are set n hours before the start of the actual hour. k can then be expressed as ˆ k |W1 , . . . , Wk−n , k = Lˆ k |L1 , . . . , Lk−n − W
Hence,
2
(
the
normal
(L , . . . , Lk−n , Lk )
vector
k−n+1
k+11
have
means
(t)dt, . . . , k−n (t)dt, k (t)dt) and covariances 1 given by (34) and (35). Thereby, Lˆ kn , which is defined by (30), can be calculated by using the normal conditional distribution fLk |L1 =l1 ,...,Lk−n =lk−n (x) =
(28)
fLk ,L1 ,...,Lk−n (x, l1 , . . . , lk−n ) fL1 ,...,Lk−n (l1 , . . . , lk−n )
(36)
,
ˆ k denotes the forecasted base load and wind power where Lˆ k and W ˆ n as the forecasted production, respectively. Introducing Lˆ kn and W k hourly base load and wind power production for hour k performed n + 1 hours before hour k (hence information regarding hour n is available, but not for hour n + 1), (28) becomes
whose parameters now are known. As for the load, the forecasted wind power production is modeled as the conditional expectation of the process,
ˆ n. k = Lˆ kn − W k
Since M(t) = ln W (t) have the same model structure as L(t), the same method as for E[Lk |L1 , . . . , Lk−n ] can be applied for the calˆn= culation of E[Mk |M1 , . . . , Mk−n ]. The forecast of Wk is then W k eE[Mk |M1 ,...,Mk−n ] . ˆ n have been calculated, the When the hourly values k = Lˆ kn − W k
(29)
The load forecasting problem is depicted in Fig. 5, which illustrates an example when n = 2. The shaded area represents the known information consisting of previous hourly loads up to hour k − 2. The calculation of k require some assumption regarding the forecasts of the load and wind power production. Starting with the load, the assumption is made that the forecast equals the conditional expectation Lˆ kn = E[Lk |L1 , . . . , Lk−n ].
(30)
This can be found be studying the energy process
+
T
S(T ) =
L(s)ds = 0
0
T
+ e−˛T
(31)
(e˛T − e˛s )dB(s) . 0
Observe that s = 0 is when the measurements started. Hence, T T is very large, implying that X(0)(1 − e−˛T ) e−˛T 0 (e˛T − e˛s )dB(s). Thereby, the part originating from the initial distribution of X(0) can be neglected. Hence,
T
S(T ) ≈
(s)ds
0
+
−˛T e ˛
(32)
T
(e˛T − e˛s )dB(s). 0
(S(u), S(t)) is a normal vector with expectations E[S(t)] = Thereby t (s)ds. The covariance can be calculated as follows prerequisite 0 that the approximation in (32) holds:
Cov[S(u), S(t)] =E
−˛u e ˛
u
(e˛u − e˛s )dB(s)
0
linear segments of the curve (t) are chosen so that k , i.e.
1 1 2
rk−1 (tk−1 )dtk−1 +
1 2
rk (tk )dtk = k ,
(37)
k+1 k
(t)dt =
(38)
0
which can be expressed as
T
(s)ds
1 X(0)(1 − e−˛T ) ˛
ˆ n = E[Wk |W1 , . . . , Wk−n ]. W k
×
−˛t e ˛
t
(e˛t − e˛s )dB(s)
.
0
(33)
Abk−1 + Bbk + Cbk+1 = k ,
(39)
where A, B and C are constants depending only on s and e . Thus, by solving the system of linear equations defined by (38), values of bk can be found and (t) is defined as the piecewise curve with segments rk (tk ). The production model defined by (15) also includes the stochastic process Y (t), which in turn includes the stochastic pairs (Tn , Vn ). The distributions of the arrival times Tn+1 − Tn ∈ exp( (Q (t)) and production decreases Vn can either be estimated from historical data of the real-time operation, or by studying the actual system and the generation units within it. Getting access to information concerning all the on-line generation units might not however be possible. One might then be forced to rely on historical data, which does not necessarily reflect the total situation in the system. For example, rare events might not be represented in historical data, or recent investments in units can have changed the distributions of Tn and Vn . Finally, the mean reversion parameter ˇ also need to be set. To estimate this parameter outgoing from historical data is certainly a challenge. Instead ˇ is set heuristically so that the reversion of the deviation from the function (t) is “reasonably” quick.
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
9. Case study To show the feasibility of the developed models, a case study consisting of a Monte Carlo simulation has been performed. The case study evaluates the impact of wind power on the demand for real-time balancing power in a fictive system. To get reasonable load curves, etc., available data from the Swedish power system has been used in the estimation process. These data consists of hourly load levels ln for the whole year of 2005 and hourly wind power production levels wn from [41] for one year. Hence, N = 8760. The simulations in this case study spans over one day and thereby K = 24. This gives NK = 365. The wind power data from [41] was scaled to represent 1000 MW of installed wind power. Concerning the distributions of Tn and Vn , some information regarding the available units in the system and times between failures was available from [17]. The assumed market structure allows for trading up to one hour before the start of the operational hour. This means that information regarding the hour starting two hours before the operational hour is assumed known, while information for the hour before the operational hour is not. Thereby, n = 2 was used for the forecasts of the hourly load and wind power production levels. The developed models and estimation methods were implemented by using the software package Matlab. 9.1. Parameter estimation The model parameters were estimated according to the methods described in Section 8. The parameter values of (t) and (t) are not presented due to limitations in space but are given on request. The parameters of the stochastic process X(t) of the base load were estimated to ˛ = 0.10, = 75.
(40)
The parameters of W (t) were estimated to ω = 5.3, = 0.01,
(41)
= 0.1. For the base production, the times s and e were set to s = 0.25 and e = 0.75. This implies that the production will start to ramp up or down 15 min before the hour starts and will end to increase/decrease 15 min after the start of the hour. Since the
973
market structure allows for trading up to one hour before the operational hour starts, the hourly levels k can be expressed as ˆ 2 . In this case study, Lˆ 2 and W ˆ 2 are for simplicity calk = Lˆ k2 − W k k k ˆ 2 = E[Wk |Wk−2 ]. culated as Lˆ 2 = E[Lk |Lk−2 ] and W k
k
The intensity of Tn+1 − Tn , i.e. (Q (t)) and the distribution of Vn are assumed independent of the current production level Q (t). Thus
(Q (t)) ≡ and FV (·, Q (t)) ≡ FV (·). was set to = 24 and hence, Tn+1 − Tn ∈ Exp(24). The density function fV (·) is assumed triangular with distribution function
fV (x) =
⎧ 1 ⎪ ⎨ fV , x ∈ [100, 500] ⎪ ⎩
fV2 ,
x ∈ [500, 1000]
0,
otherwise,
(42)
where fV1 (x) = −5.56 × 10−4 + 5.56 × 10−6 x,
(43)
fV2 (x) = 4.44 × 10−3 − 4.44 × 10−6 x.
(44)
Hence, the size of the shut-downs Vn will exist on the interval [100, 1000]. Further, the reversion parameters was set to ˇ = 24. 9.2. Simulations When the model parameters of L(t), W (t) and Q (t) have been set, sample path simulations of these processes can be performed and the associated sample paths of D(t) = Q (t) − N(t) = Q (t) − (L(t) − W (t)) can be calculated. In the simulations, two different cases were run: (1) With wind power, and (2) without wind power. The simulated sample paths of the load L(t) are the same for the both cases, while W (t) only were applied in case 1. Since Q (t) is dependent of forecasts of L(t) and W (t), different simulations of Q (t) were performed for the two cases. To distinguish between these two, the denotation Q #1 (t) refers to the net base production in case 1 and Q #2 (t) for case 2. The corresponding denotation for D(t) is D#1 (t) and D#2 (t). The presented models were used to generate 100 sample paths, each covering 24 h, of L(t), W (t), Q #1 (t) and Q #2 (t). The sample paths of L(t), W (t) and Q #1 (t) are displayed in Fig. 6. Note that since the algebraic function #1 (t) is dependent on successive forecasts of Lk and Wk , the #1 (t) s are different for each sample path. Concerning the differences between the cases with and without wind power, the associated hourly energy process of D(t), defined
k+1
as Dk = k D(t)dt, was computed for the various sample paths. The histograms for case 1 and case 2 of Dk are shown in Fig. 7. The
Fig. 6. Sample paths of L(t) (left), W (t) (middle) and Q #1 (t) (right). The black lines corresponds to the algebraic functions (t), eω(t) and #1 (t). Note that for Q #1 (t) there is one algebraic function per sample path.
974
M. Olsson et al. / Electric Power Systems Research 80 (2010) 966–974
Fig. 7. Histograms of dk#1 (black bars) and dk#2 (white bars).
sample mean and variance of Dk were also calculated: E[Dk#1 ] = −15.7, E[Dk#2 ] = −18.9, Var[Dk#1 ] = 6558,
(45)
Var[Dk#2 ] = 6332. Hence, with wind power, the mean moves slightly towards zero while the variance increases. 10. Conclusion This paper presents models based on SDEs of base production, base consumption and wind power production. These models can be used to evaluate the impact of wind power on the real-time balancing of the power system. This is of great interest for system operators as well as for producers and consumers. Using continuous time models enables the possibility of considering counter balancing. Besides the description of the models, estimation methods for parameters are proposed and a case study is presented. The case study suggests that using SDEs is a possible approach for modeling imbalances in continuous time. Future research concerning this topic could consist of performing more detailed case studies where the developed models are applied. Acknowledgement Support from the Nordic power exchange Nord Pool considering market data is gratefully acknowledged. References [1] R. Weron, Modeling and Forecasting Electricity Loads and Prices—A Statistical Approach, John Wiley & Sons Inc., 2006. [2] N. Amjady, Short-term hourly load forecasting using time-series modeling with peak load estimation capability, IEEE Trans. Power Syst. 16 (3). [3] E. Lindgren, L. Söder, Minimizing regulating costs in multiarea systems with uncertain wind power forecasts, Wind Energy 11 (1) (2008) 97–108.
[4] H. Banakar, C. Luo, B.T. Ooi, Impacts of wind power minute-to-minute variations on power system operation, IEEE Trans. Power Syst. 23 (1) (2008). [5] K. Elkington, V. Knazkins, M. Ghandhari, On the stability of power systems containing doubly fed induction generator-based generation, Electric Power Syst. Res. 78 (9). [6] I. Norheim, K. Uhlen, J.O. Tande, T. Toftevang, M.P. Pálsson, Doubly fed induction generator model for power system simulation tools, in: Proceeding of the NordicWind Power Conference, Gothenburg, Sweden, 2004. [7] T. Petru, T. Thiringer, Modeling of wind turbines for power system studies, IEEE Trans. Power Syst. 17 (4). [8] V. Akhmatov, H. Knudsen, An aggregated model of a gridconnected, large-scale, offshore wind farm for power system stability investigation—importance of windmill mechanical system, Int. J. Electrical Power Energy Syst. 24 (9) (2002). [9] H. Holttinen, Impact of hourly wind power variations on the system operation in the Nordic countries, Wind Energy 8 (2). [10] Wind power plants and system operation in the hourly time domain, NREL, 2003. URL: http://www.nrel.gov. [11] B. Oksendal, Stochastic Differential Equations, Springer-Verlag, Berlin Heidelberg, 2005. [12] F.C. Klebaner, Introduction to Stochastic Calculus with Applications, Imperial Collage Press, 2005. [13] European Energy Exchange. URL: http://www.eex.com. [14] Scheduling operations, PJM Manual 11: PJM, 2007. URL: http://www.pjm.com. [15] Nordic grid code 2004. Nordel, 2004. URL: http://www.nordel.org. [16] The Grid Code, National Grid Electricity Transmission Plc., 2008. URL: http://www.nationalgrid.com/. [17] Nord Pool AS web site. URL: http://www.nordpool.com. [18] Operador del Mercado Ib’erico de Energia web site. URL: http://www.omel.es. [19] APX Group. URL: http://www.apxgroup.com. [20] Market Operations, ISO New England Inc., 2006. URL: http://www.iso-ne.com. [21] Operating Procedure: Frequency Control Ancilliry Services, National Electricity Market Management Company (NEMMCO), 2005. URL: http://www.nemmco.com.au. [22] Current State of Balance Management in Europe, ETSO, 2003. URL: http://www.etso-net.org. [23] D. Ranaweera, G.G. Karady, R.G. Farmer, Economic impact analysis of load forecasting, IEEE Trans. Power Syst. 12 (3) (1997) 1388. [24] M.A. Abu-El-Magd, N.K. Sinha, Short-term load demand modeling and forecasting: a review, IEEE Trans. Syst., Man, Cybern. SMC-12 (3) (1982) 370. [25] I. Moghram, S. Rahman, Analysis and evaluation of five shortterm load forecasting techniques, IEEE Trans. Power Syst. 4 (4). [26] S. Fan, L. Chen, Short-term load forecasting based on an adaptive hybrid method, IEEE Trans. Power Syst. 21 (1) (2006) 392–401. [27] N. Amjady, Short-term bus load forecasting of power systems by a new hybrid method, IEEE Trans. Power Syst. 22 (1) (2007) 333–341. [28] Z. Yun, Z. Quan, S. Caixin, L. Shaolan, L. Yuming, S. Yang, RBF neural network and ANFIS-based short-term load forecasting approach in real-time price environment, IEEE Trans. Power Syst. 23 (3) (2008) 853–858. [29] M. Milligan, A.H. Miller, F. Chapman, Estimating the economic value of wind forecasting to utilities, URL: http://www.nrel.gov/publications/, in: Proceedings of Windpower’ 95, Washington DC, USA, 1995. [30] S.J. Watson, G. Giebel, A. Joensen, The economic value of accurate wind power forecasting to utilities, in: Proceedings of European Wind Energy Conference, Nice, France, 1999. [31] T.S. Nielsen, Prediction of regional wind power, in: Proceedings of Global Wind Power Conference, Paris, France, 2002. [32] U. Focken, M. Lange, H.M. Waldl, Previento—a wind power prediction system with an innovative upscaling algorithm, in: Proceedings of European Wind Power Conference, Copenhagen, Denmark, 2001. [33] G. Giebel, R. Brownsword, G. Kariniotakis, The state-of-theart in short-term prediction of wind power. A literature review, Tech. Rep., EU Project ANEMOS. URL: http://anemos.cma.fr. [34] M. Perninge, Evaluating the uncertainties involved in net transmission capacity calculation, Licentiate Thesis, School of Electrical Engineering, KTH, Stockholm, Sweden, 2009. [35] F. Dornier, M. Queruel, Caution to the wind, Weather risk special report, Energy & power risk management, Risk Magazine, 2000. [36] N. Ahmed, T. Natarajan, K.R. Rao, Discrete cosine transform, IEEE Trans. Comput. C-23. [37] K. Conradsen, L.B. Nielsen, L.P. Prahm, Review of Weibull statistics for estimation of wind speed distributions, J. Clim. Appl. Meteorol. 23. [38] J. Matevosyan, L. Söder, Minimization of imbalance cost trading wind power on the short-term power market, IEEE Trans. Power Syst. 21(3). [39] G. Papaefthymiou, B. Klöckl, MCMC for wind power simulation, IEEE Trans. Energy Convers. 23 (1). [40] S. Ditlevsen, M. Sorensen, Inference for observations of integrated diffusion processes, Scand. J. Stat. 31. [41] U. Arvidsson, R. Murray, V. Neimane, 4000 MW wind power in Sweden: impact on regulation and reserve requirements, Tech. Rep. 05:19, Elforsk AB, 2005.