A semi-Markov model of the variability of power generation from renewable sources

No. 2 2013 DOI: 10.5277/ord130207

Jacek MALINOWSKI*

A SEMI-MARKOV MODEL OF THE VARIABILITY OF

POWER GENERATION FROM RENEWABLE SOURCES

The paper presents a new approach to modeling the variability of power generation from a re-newable source such as wind or flowing water. The force of the power generating agent is assumed to change according to the semi-Markov process with finite state space. For the purpose of its construc-tion, the range of possible values expressing the agent’s force is divided into a finite number of subin-tervals. It is natural to assume that the time during which the agent’s force remains within one such interval, and the probabilities of transitions to neighboring intervals depend to some extent on the agent’s earlier behavior. The model’s accuracy is determined by the number of subintervals used and the assumed degree to which the agent’s force depends on its history. This degree is expressed by the number of the most recently entered subintervals relevant to predicting the agent’s future behavior. According to the presupposed accuracy level, an appropriately complex state-space and a diagram of the inter-state transitions for the modeled process have been constructed. Subsequently, it is demon-strated how certain parameters of this process, related to forecasting power generation, can be calcu-lated by means of the calculus of the Laplace transforms.

Keywords: renewable power generation, forecasting of environmental conditions, semi-Markov

model-ing, Laplace transform

1. Introduction

There exist a considerable number of models for predicting natural phenomena in the context of electrical power generation. Due to their nature, such phenomena can only be predicted with some degree of uncertainty, hence their modeling is based on various mathematical ways of expressing non-deterministic, uncertain information. Admittedly, deterministic models are also widely used, e.g. in weather forecasting, but their implementation based mainly on systems of partial differential equations requires _________________________

*Systems Research Institute of the Polish Academy of Sciences, ul. Newelska 6, 01-447 Warsaw, Poland, e-mail: jacek.malinowski@ibspan.waw.pl

82 a on for ca lat Tr or the po pre be bil est ing bil by ran res du of the ter substantial amo nly suitable for

re, non-determi sting, since it is ted over mediu raditionally, pre statistical tools eory of possibi opular in recent esented in Fig. This paper is ehavior of renew listic modeling. tablished metho g, time series an lity of a power g y a stochastic p

nge of the force sponding to ene uring which the f transitions to n

e sequence of i rval. The mode

ount of comput short term pred inistic predictio s not much less um or long time ediction of the f s (cf. [2, 3, 5, ility, neural ne t years (cf. [4, 1. Fig. 1. An concerned with wable energy so . Short-term for ods have been nalysis). A nov generating agen rocess with a f e of the agent is ergy production agent’s force r neighboring int intervals in whi l’s accuracy is J.MALINOWS

ting power and diction – coveri on offers a plau accurate – as lo e periods, and former type is c 8, 9]). Howeve etworks or wav 6]). An overvi overview of foreca h the issue of m ources such as w recasting is not developed for t vel stochastic m nt over time. Th finite state spac s divided into a n levels) denote remains within tervals depend ich its force had

determined by SKI d memory spac ing a period of usible alternativ ong as certain a significantly le carried out with er, methods emp velet theory hav iew of the fore

asting methodology medium or long wind or flowing addressed herei this purpose (e model is presente he fluctuations o ce. In order to c finite number o ed I1, …, Im. It i

one such interv on the agent’s d been before i the total numb

ce. In addition, up to one week ve to determini average values a ess resource con h the use of prob ploying fuzzy l ve become incr ecasting method y g term forecastin g water, based o in, as a plethora .g. deterministi ed to describe t of its force are d construct this s of subintervals is assumed that val, and the prob

earlier behavio it entered the cu er of these inte they are k. There-stic fore-are calcu-nsuming. babilistic logic, the reasingly dology is ng of the on proba-a of well-ic model-the varia-described pace, the (e.g. cor-t cor-the cor-time babilities or, i.e. on urrent in-ervals (m)

an nu be rec of va the mo Z sta rea ing ac {z pe ou nd the assumed umber of the mo ehavior. A state cently entered o f the agent’s forc

As it is assum al Ix, and the nu e order of sever odeled by a se = {z1, …, zn}. ate contains in ached Ix. The nu

g the state spac curacy. In cons 1, ..., zn} exceed F when t If the degree ends on whether us interval was degree of depe ost recently ent e of the process ones, thus it inc ce. Further deta

2.

med that the tim umber of the ne

ral of the most emi-Markov pro Each interval I nformation on t umber of previo ce Z (denoted a sequence, if we ds that of {I1, ...

Fig. 2. The state sp the degree of depen

of dependence r its value has r

Ix–1 or Ix+1, wh

endence of the tered intervals r s is determined corporates infor ails are given in

The model’s

me during which ext interval (it c

recently entere ocess with the

Ix corresponds t the intervals en ously entered in as d) determines e neglect the de , Im}, so that n > ace architecture ndence is 1 and m = e is assumed to recently increas here Ix is the cu force on its hi relevant to pred by the current rmation relevan the next chapte

s basics

h the agent’s va can be either x – ed intervals, the appropriately to a number of ntered by the ntervals taken in s, along with m egenerate case m > m. = 6 be 1, the futur ed or decreased urrent interval. istory (d) define dicting the agent interval and th nt to the future e er. alue remains in – 1 or x + 1) de e agent’s change constructed sta f states in Z, wh agent’s force nto account in c m, the model’s d m = 2, the cardi re force of the a d, i.e. whether t

The state spac

ed as the t’s future he d most evolution an inter-epend on es can be ate space here each before it construct-degree of inality of agent de-the previ-ce of the

84 mo nu nu an Ix+ sp ag i.e no z10 z13 ha cre the sen ing ha dia odeling process umbered state i umbered one – nd current interv +1 and Ix, then th

ace along with If the degree gent depends no e. the last two i ow composed o

0, ..., z4m–14}, G2 3, ..., z4m–8, z4m–

as most recently eased and then en increased; in nts the consider Fig. 3. Th Let Z = {Zt, t g the agent’s va ave to be defined X = {Xk, k ≥ 0

ately after its

k-s ik-s compok-sed is entered if th if it has most r vals, then the pr he process is in the possible int e of dependence ot only on the m intervals are sig f 4m – 6 states 2 = {z2, z4, z7, z11 –6}. The states i y decreased tw decreased; in t n the fourth – if red state space a

he state space archi

t ≥ 0} be a stoc ariability. As Z d: 0} – the embed -th change in sta J.MALINOWS of n = 2m – 2 he force has m recently decreas rocess is in the s the state z2x–1. F ter-state transitio e is assumed to most recent, but

gnificant. The s, which are div

1, ..., z4m–10}, G3

in the first grou wice; in the sec

the third – if it f it has most re along with the p

itecture when the d

hastic process w

Z is assumed to

dded Markov ch ate (X0 is the ini SKI states denoted most recently i sed. Thus if Ix–1 state z2x–2; if the Figure 2 represe ons for m = 6. o be two, then

also on the pre state space of t vided into four

= {z5, z8, z12, ..

up are entered i ond – if it has t has most rece ecently increase possible intersta

degree of dependen

with the state sp be semi-Marko hain of Z, i.e. Xk

itially observed

d z1, …, z2m–2. A

increased, and

1 and Ix are the

e respective inte ents the conside the future forc eceding change the modeling p groups: G1 = { .., z4m–9, z4m–7}, if the value of s most recently ently first decre ed twice. Figure ate transitions fo

nce is 2 and m = 6

pace {z1, ..., zn}

ov, the followin

Xk is the state of d state of Z). An even-an odd-previous ervals are ered state ce of the in value, process is {z1, z3, z6, G4 = {z9, the force y first in-eased and e 3 repre-or m = 6. }, model-ng objects f Z

imme-P = [pij], i, j = 1, ..., n – the transition matrix of X, i.e. pij = Pr(Xk = zj|Xk–1 = zi). It is

assumed that P does not change with k, i.e. X is homogenous.

Tk – the moment of the k-th change in the state of Z.

Sij – the time spent by Z in the state zi given that the next state is zj. Fij – the distribution function of Sij, i.e.

( )

Pr

(

1 , 1

)

ij k k k j k i

F t = T −T₋ ≤t X =z X ₋ =z (1)

P and Fij can be obtained from the statistical analysis of the recorded values of the

agent’s force. The matrix P corresponding to the state space shown in Fig. 2 is pre-sented below. 21 24 31 34 43 46 53 56 65 68 75 78 87 8,10 97 9,10 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 p p p p p p p p p p p p p P p p p ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

Clearly, pii = 0 and pi1 + … + pin = 1 for i = 1, ..., n – these are obvious conditions

that must be fulfilled by any transition matrix.

Remark: In general, the probability of a transition from zi to zj depends on the

amount of time spent in zi, thus

(

1 1

)

(

1 1

)

Pr X_k =z X_{j k−} =z T_i, _k −T_k− ≤ ≠t Pr X_k =z X_{j k−} =z T_i, _k −T_k− ≤u (2) may hold for u ≠ t. The precise definition of pij should therefore be as follows:

(

)

1 1

Pr ,

ij k j k i k k

p = X =z X ₋ =z T −T₋ < ∞ (3)

This underscores the fact that pij is the probability of transition from zi to zj if no

J.MALINOWSKI

86

irrelevant to pij as expressed by (3), because Pr(Tk – Tk–1 < ∞) = 1. It should be noted

that this remark pertains to all semi-Markov processes.

3. Equations for basic parameters of the energy production process

The model constructed in the previous section is useful for determining many pa-rameters that characterize the process of electrical power production, especially for predicting future power output. For example, it is possible to forecast the total ex-pected energy output during a given time period, the probability that during such a period the output value will remain within certain limits, or the expected number of times it will cross these limits. It will now be shown in detail how to determine the first parameter. Let πi be the power generated when Z is in the state zi, i.e. when the

factor’s force is in the corresponding interval. Let Gi(u, t) be the expected amount of

energy produced in the time interval [u, t], given that at the moment u the process Z enters the state zi. It can be easily shown that the Gi(0, t), i = 1, ..., n satisfy the

follow-ing set of equations:

( )

( )

( )

( )

0 0, 1 , t i i ij ij ij i j ij j i j i G t πt p F t p πu G u t dF u ≠ ≠ ⎡ ⎤ ⎡ ⎤ =

∑

_⎣ − _⎦+

∑ ∫

_⎣ + _⎦ (4)

Indeed, taking into account that

( )

(

₍

₎

)

(

₍

₎

)

(

)

1 1 1 1 1 1 1 Pr , Pr , , Pr Pr , Pr , k j k i k k k j k i ij ij k i k j k i k k k j k i X z X z T T t X z X z p F t X z X z X z T T t X z X z − − − − − − − = = − ≤ = = = = = = = − ≤ = = (5)

the first component on the right-hand side of (4) is related to the amount of energy produced in the time interval [0, t], given that no change in state occurred from 0 to t. In turn, the second component is related to the amount of energy produced in that in-terval, given that the first change in state occurs at u ≤ t.

Being a semi-Markov process, Z “forgets” its history at each change in state, so that Gi(s, t) = Gi(0, t – s), i = 1, ..., n. In consequence, (4) can be transformed to:

( )

(

( )

)

( )

(

)

( )

0 0 0, 1 0, t t i i ij ij ij ij j ij j i j i G t π p t F t udF u p G t u dF u ≠ ≠ ⎡ ⎤ = _⎢ − + _⎥+ − ⎣ ⎦

∑

∫

∑

∫

(6)

Let us note that

( )

( )

( )

0 1 min , t ij ij ij t⎡_⎣ −F t ⎤_⎦+

_∫

udF u =E⎡_⎣ S t ⎤_⎦ (7)

where E is the expected value. Hence, (6) can be written in the following compact form:

( )

( )

(

)

( )

0 0, 0, t i i ij ij ij j ij j i j i G t π p H t p G t u dF u ≠ ≠ =

∑

+

∑ ∫

− (8) where

( )

min

( )

, ij ij H t _{= ⎣}E⎡ S t ⎤_{⎦ (9)}

As (8) is a system of integral equations, the calculus of Laplace transforms can be used to find its solution. Let

( )

{

( )

0,

}

i s Gi t Γ = L (10)

( )

{

( )

}

*

{

( )

}

ij s f tij F tij Φ =L =L (11)

( )

{

( )

}

ij s H tij Ψ = L (12)

where

ࣦ

and

ࣦ

* denote the Laplace and Laplace–Stieltjes transforms, respectively, and fij is the probability density function of Sij. From the basic properties of the

La-place transform, it follows that

( )

2

( )

1 _ij ij s s s Φ Ψ = ⎡⎣ − ⎤⎦ (13)

Applying

ࣦ

to both sides of (8), we obtain

( )

( )

( ) ( )

i i ij ij ij ij j j i j i s p s p s s Γ π Ψ Φ Γ ≠ ≠ =

∑

+

∑

(14)

Note that the integral in (8) is the convolution of Gj and fij, which, after applying

the Laplace transform, becomes the product of Γj and Φij. As (14) is a system of linear

algebraic equations, it can be written in the following matrix form:

( )

( )

( )

( )

( )

1 1 1 1 j j j i n n nj nj j i p s s A s s p s π Ψ Γ Γ π Ψ ≠ ≠ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢₌ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ _{⎦ ⎢ ⎥} ⎣ ⎦

∑

∑

(15)

J.MALINOWSKI

88

where the elements of the matrix A(s) are given by:

( )

( )

ij ij ij ij

a s =δ −pΦ s (16)

4. Solving the equations obtained

A closed-form solution of (15) should express Γi, i = 1, ..., n in terms of πi, pij, Φij,

and Ψij, where i, j = 1, ..., n. However, it is practically impossible to find such a

solu-tion in general. For example, in the case of the system whose state space is shown in Fig. 2, we have the following matrix A(s):

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

( )

12 21 21 24 24 31 31 34 34 43 43 53 53 65 65 75 75 46 46 56 56 68 68 78 78 87 8 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 s p s p s p s p s p s p s A s p s p s p s p s p s p s p Φ Φ Φ Φ Φ Φ Φ Φ Φ Φ Φ Φ Φ Φ ⎡ − ⎢− − ⎢ ⎢− − ⎢ − ⎢ ⎢ ₋ ⎢ = − ⎢ ⎢ ₋ ⎢ ⎢ ⎢ ⎢ ⎢⎣ − − − − −

( )

( )

( )

( )

( )

7 8,10 8, 10 97 97 9,10 9,10 10,9 1 0 0 0 1 0 0 0 1 s p s p s p s s Φ Φ Φ Φ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ − ⎥ ⎥ − − _⎥ ⎥ − _⎦

A(s) is a sparse matrix (it includes many zeros) and can be transformed into

rearrangement of the elements of [Γ1(s), ..., Γn(s)]T). Thus the system (15) can be

solved by Gaussian elimination (GE) with lower than maximal complexity, which is

O(n3_{). However, even the Thomas algorithm (a simplified form of GE designed}

espe-cially for tri-diagonal matrices), applied to (15) with A(s) as given above, results in very complicated closed formulas, which are practically impossible to derive for large n. Needless to say, these formulas become even more complicated for d > 1.

From the above considerations, (15) has to be solved numerically rather than ana-lytically. For example, it is possible to find Γi(s) for a number of discrete values of s

along a vertical line in the space of complex numbers, thus obtaining data which allow us to compute Gi(0, t) by numerical integration using the following well-known

for-mula for the reverse Laplace transform:

( )

0, 1 lim ( ) 2 x iy st i _y i i x iy G t e Γ s ds π + →∞ − =

_∫

(17)

Here, x is the real part of the complex number x + iy.

5. Finding other parameters

We will now briefly outline how other parameters, e.g. the probability that during a given time period the generated power will remain within certain limits, or the ex-pected number of times it will cross these limits, can be found in a way similar to that described in Sections 3 and 4.

Let Ni(u, t) be the expected number of times that from u to t the generated power

output crosses the lower bound of Ia from above, or the upper bound of Ib from below,

where a ≤ b, given that at the moment u the process Z enters the state zi. Let Z1 = {z∈Z: a ≤ x(z) ≤ b} and Z2 = Z \ Z1. As in the case of Gi(0, t), i = 1, ..., n (see

Sec-tion 3), it can be shown that Ni(0, t), i = 1, ..., n fulfill the following set of equations:

( )

( )

( )

( ) 0 0, , , , 1, , t i ij j ij j i N t p γ i j N u t dF u i n ≠ ⎡ ⎤ =

∑ ∫

_⎣ + _⎦ = … (18)

where γ(i, j) = 1 for i ∈ Z1 and j ∈ Z2, otherwise γ(i, j) = 0.

Let Pi(u, t) be the probability that from time u to time t the generated power output

remains above the lower bound of Ia and below the upper bound of Ib, i ∈ Z1. The

fol-lowing set of equations is obtained for Pi(0, t), i = 1, ..., n:

( )

( )

( )

( )

(

( )

)

₁ 0 0, , , 1 , t i ij j ij ij j i P t p δ i j P u t dF u F t i Z ≠ ⎡ ⎤ = _⎢ + − _⎥ ∈ ⎣ ⎦

∑

∫

(19)

J.MALINOWSKI

90

where δ(i, j) = 1 for j ∈ Z1, δ(i, j) = 0 for j ∈ Z2. Note that in (19) the number of

un-knowns is equal to the number of equations, which is equal to card(Z1). Thus a

neces-sary condition for the existence of a unique solution is fulfilled.

Both (18) and (19) can be transformed to their “convolution” versions, analogous to (8), and then solved using transform calculus, as shown in Section 4.

6. Conclusion

A new model of the stochastic variability of a renewable energy source has been presented. This model can be used to determine the number of parameters characteriz-ing the energy output over a medium or long time horizon. As the model is purely stochastic, a question may arise regarding its applicability, especially when the power generating agent is strongly dependent on some deterministic factor. For example, ocean tides, made use of in tidal power stations, depend on the positions of the Sun and the Moon in relation to the Earth. In such cases, the proposed model might de-scribe the stochastic component of the mainly deterministic variability of a given agent.

References

[1] AI B.,YANG H.,SHEN H.,LIAO X., Computer-aided design of PV/wind hybrid system, Renewable

Energy, 2003, 28, 1492–1512.

[2] ALBADI M.H.,EL-SAADANY E.F., Overview of wind power intermittency impacts on power systems,

Electric Power Systems Research, 2010, 80, 627–632.

[3] CARTA J.A.,RAMIREZ P.,VELAZQUEZ S., A review of wind speed probability distributions used in

wind energy analysis, Renewable and Sustainable Energy Reviews, 2009, 13, 933–955.

[4] JAFARIAN M.,RANJBAR A.M.,Fuzzy modeling techniques and artificial neural networks to estimate

annual energy output of a wind turbine, Renewable Energy, 2010, 35, 2008–2014.

[5] KARKI R.,PO H., BILLINTON R., A simplified wind power generation model for reliability evaluation,

IEEE Transactions on Energy Conversion, 2006, 21, 533–540.

[6] KULKARNI M.A.,PATIL S.,RAMA G.V.,SEN P.N.,Wind speed prediction using statistical regression

and neural network, Journal of Earth System Science, 2008, 117, 457–463.

[7] LIN Z.,ZHANG D.,GAO L.,KONG Z.,Using an adaptive self-tuning approach to forecast power loads,

Neurocomputing, 2008, 71, 559–563.

[8] LIU H.,TIAN H.-Q.,CHEN C.,LI Y.-F., A Hybrid Statistical Method to Predict Wind Speed and Wind

Power, Renewable Energy, 2010, 35, 1857–1861.

[9] XIAO Y.Q.,LI Q.S.,LI Z.N.,CHOW Y.W.,LI G.Q., Probability distributions of extreme wind speed